EigenUnitTest.cpp File Reference

#include <Eigen/Sparse>
#include <vector>
#include <iostream>
#include "Helpers/Helpers.h"

Typedefs
typedef Eigen::SparseMatrix< double >	mat
typedef Eigen::Triplet< double >	T

Functions
int	main (int argc, char **argv)

Variables
std::vector< T >	tripletList

Typedef Documentation

mat

typedef Eigen::SparseMatrix<double> mat

T

typedef Eigen::Triplet<double> T

Function Documentation

main()

int main	(	int argc	,
		char **	argv
	)

                               {
    double L = 1;
    int varnum = 5;
    double dx = L/(varnum-1);
    double dt = pow(dx,2);//does not need to be this, but keeps it anyway
    double A = dt/2.0/pow(dx,2);
    double B = -(1+A*2);
 
 
    //construct coefficient sparse matrix
    //diagonal elements
    for(int i=0; i<varnum; i++) {
        tripletList.push_back(T(i, i, B));
    }
 
    //right-diagonal elements
    for(int j=1; j<varnum; j++){
        tripletList.push_back(T(j-1,j,A));
    }
    //left-diagonal elements
    for(int k=0; k<varnum-1; k++){
        tripletList.push_back(T(k+1, k, A));
    }
 
    //distribute the elements to the matrix
    mat C(varnum,varnum);
    C.setFromTriplets(tripletList.begin(), tripletList.end());
 
    //construct RHS vector
    Eigen::VectorXd T(varnum), K(varnum);
    for(int i=0; i<varnum; i++)
        T(i) = 0;
 
    for(int j=1; j<varnum-1; j++){
        K(j) = -T(j)-A*(T(j+1)-2*T(j)+T(j-1));
    }
    K(0) = -T(0)-A*(T(1)-2*T(0));
    K(varnum-1) = -T(varnum-1)-A*(1-2*T(varnum-1)+T(varnum-2))-A;
 
 
    //iterate with time
    Eigen::ConjugateGradient<mat> solver;
    for(double t=0; t<dt*100; t+=dt){
        solver.compute(C);
        T = solver.solve(K);
 
        for(int j=1; j<varnum-1; j++){
            K(j) = -T(j)-A*(T(j+1)-2*T(j)+T(j-1));
        }
        K(0) = -T(0)-A*(T(1)-2*T(0));
        K(varnum-1) = -T(varnum-1)-A*(1-2*T(varnum-1)+T(varnum-2))-A;
 
        //output the values at the grid points to the console
        std::cout << "t=" << t << "s" << "\t" << T(0) << "\t" << T((varnum-1)*1/4) << "\t" << T((varnum-1)/2) << "\t" << T((varnum-1)*3/4) << "\t" << T(varnum-1)<< std::endl;
    }
    //check last timestep
    helpers::check(T(0),0.166667,1e-5, "Gridpoint 0");
    helpers::check(T((varnum-1)*1/4),0.333333,1e-5, "Gridpoint 1/4");
    helpers::check(T((varnum-1)*2/4),0.5,1e-5, "Gridpoint 2/4");
    helpers::check(T((varnum-1)*3/4),0.666667,1e-5, "Gridpoint 3/4");
    helpers::check(T(varnum-1),0.833333,1e-5, "Gridpoint 1");
}

References helpers::check(), e(), i, j, k, PlanarWave::K, L, Eigen::ArrayBase< Derived >::pow(), solver, plotPSD::t, and tripletList.

Variable Documentation

tripletList

std::vector<T> tripletList

Referenced by main().