fem2d Namespace Reference

Functions
void	assemble_heat_problem (const numcxx::SimpleGrid &Grid, const numcxx::DArray1 &BCfac, const numcxx::DArray1 &BCval, const numcxx::DArray1 &Source, const numcxx::DArray1 &Kappa, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Rhs)
void	assemble_transient_heat_problem (const numcxx::SimpleGrid &grid, const numcxx::DArray1 &bcfac, const numcxx::DArray1 &bcval, const numcxx::DArray1 &source, const numcxx::DArray1 &kappa, double tau, double theta, bool lump, numcxx::DArray1 &OldSol, numcxx::DSparseMatrix &SGlobal, numcxx::DArray1 &Rhs)
double	l2norm (const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
double	h1norm (const numcxx::SimpleGrid &grid, const numcxx::DArray1 &u)
void	compute_cell_volume (const int icell, const numcxx::DArray2 &points, const numcxx::IArray2 &cells, double &vol)
void	compute_local_stiffness_matrix (const int icell, const numcxx::DArray2 &points, const numcxx::IArray2 &cells, numcxx::DArray2 &SLocal, double &vol)

Variables
const double	Dirichlet =1.0e30
	BC value marking Dirichlet boundary condition. More...

Function Documentation

void fem2d::assemble_heat_problem	(	const numcxx::SimpleGrid &	Grid,
		const numcxx::DArray1 &	BCfac,
		const numcxx::DArray1 &	BCval,
		const numcxx::DArray1 &	Source,
		const numcxx::DArray1 &	Kappa,
		numcxx::DSparseMatrix &	SGlobal,
		numcxx::DArray1 &	Rhs
	)

Examples:

40-stationary-heat-fe.cxx, and 42-convtest-fem-sin.cxx.

Definition at line 74 of file fem2d.cxx.

  {
    
    auto ndim=grid.spacedim();      // space dimension
    assert(ndim==2);
    auto points=grid.get_points(); // Array of global nodes
    auto cells=grid.get_cells();   // Local-global dof map

    int npoints=grid.npoints();
    int ncells=grid.ncells();
    double vol=0.0;

    // Local stiffness matrix, 
    numcxx::DArray2 SLocal{{0,0,0},{0,0,0},{0,0,0}};
    numcxx::DArray2 MLocal0{{2,1,1},{1,2,1},{1,1,2}};
    MLocal0*=1.0/12.0;
    Rhs=0.0;
    SGlobal(0,0)=0;
    SGlobal.clear();
    double third=1.0/3.0;

    
    // Loop over all elements (cells) of the triangulation
    for (int icell=0; icell<ncells; icell++)
    {
      compute_local_stiffness_matrix(icell, points,cells, SLocal, vol);
      
      // Assemble into global stiffness matrix
      double klocal=(kappa(cells(icell,0))+kappa(cells(icell,1))+kappa(cells(icell,2)))/3.0;
      
      for (int i=0;i<=ndim;i++)
        for (int j=0;j<=ndim;j++)
        {
          Rhs(cells(icell,i))+=vol*MLocal0(i,j)*source(cells(icell,j));
          SGlobal(cells(icell,i),cells(icell,j))+=klocal*SLocal(i,j);
        }
    }    
    

    // Assemble boundary conditions (Dirichlet penalty method)
    int nbfaces=grid.nbfaces();
    auto bfaces=grid.get_bfaces();
    auto bfaceregions=grid.get_bfaceregions();
    
    for (int ibface=0; ibface<nbfaces; ibface++)
    {
      // Obtain number of boundary condition
      int ireg=bfaceregions(ibface);
      int i0=bfaces(ibface,0);
      int i1=bfaces(ibface,1);
      
      // Check if it is "Dirichlet"
      if (bcfac(ireg)>=Dirichlet)
      {
        // Assemble penalty values
        SGlobal(i0,i0)+=bcfac(ireg);
        Rhs(i0)+=bcfac(ireg)*bcval(ireg);
        
        SGlobal(i1,i1)+=bcfac(ireg);
        Rhs(i1)+=bcfac(ireg)*bcval(ireg);
        
      }
      else if (bcfac(ireg)>0.0)
      {

        double dx= points(bfaces(ibface,1),0)-points(bfaces(ibface,0),0);
        double dy= points(bfaces(ibface,1),1)-points(bfaces(ibface,0),1);
        double h=sqrt(dx*dx+dy*dy);
        
        SGlobal(i0,i0)+=h*bcfac(ireg)/3.0;
        SGlobal(i0,i1)+=h*bcfac(ireg)/6.0;
        
        SGlobal(i1,i0)+=h*bcfac(ireg)/6.0;
        SGlobal(i1,i1)+=h*bcfac(ireg)/3.0;
        
        Rhs(i0)+=0.5*h*bcval(ireg);
        Rhs(i1)+=0.5*h*bcval(ireg);
      }
    }
    SGlobal.flush();
  }

Here is the call graph for this function:

void fem2d::assemble_transient_heat_problem	(	const numcxx::SimpleGrid &	grid,
		const numcxx::DArray1 &	bcfac,
		const numcxx::DArray1 &	bcval,
		const numcxx::DArray1 &	source,
		const numcxx::DArray1 &	kappa,
		double	tau,
		double	theta,
		bool	lump,
		numcxx::DArray1 &	OldSol,
		numcxx::DSparseMatrix &	SGlobal,
		numcxx::DArray1 &	Rhs
	)

Examples:

44-transient-heat-fe.cxx.

Definition at line 165 of file fem2d.cxx.

  {
    
    auto ndim=grid.spacedim();
    auto points=grid.get_points(); // Array of global nodes
    auto cells=grid.get_cells();   // Local-global dof map
    int npoints=grid.npoints();
    int ncells=grid.ncells();
    double tauinv=1.0/tau;
    
    // Local stiffness matrix
    numcxx::DArray2 SLocal{{0,0,0},{0,0,0},{0,0,0}};
    numcxx::DArray2 MFull{{2,1,1},{1,2,1},{1,1,2}}; // from Stroud quadrature...
    numcxx::DArray2 MLumped{{4,0,0},{0,4,0},{0,0,4}}; // mass lumping
    numcxx::DArray2 MLocal(3,3);  

    if (lump)
      MLocal=MLumped/12.0;
    else
      MLocal=MFull/12.0;

    // Loop over all elements (cells) of the triangulation
    Rhs=0.0;
    SGlobal(0,0)=0;
    SGlobal.clear();
    double vol;
    for (int icell=0; icell<ncells; icell++)
    {
      compute_local_stiffness_matrix(icell, points,cells, SLocal, vol);
      double klocal=(kappa(cells(icell,0))+kappa(cells(icell,1))+kappa(cells(icell,2)))/3.0;
      
      // Quadrature for heat transfer coefficient
      double KLocal=(kappa(cells(icell,0))+kappa(cells(icell,1))+kappa(cells(icell,2)))/3.0;
      
      
      int i;
      for (int i=0;i<=ndim;i++)
        for (int j=0;j<=ndim;j++)
        {
          SGlobal(cells(icell,i),cells(icell,j))+=theta*KLocal*SLocal(i,j)+ tauinv*vol*MLocal(i,j);
          
          Rhs(cells(icell,i))+=((theta-1)*KLocal*SLocal(i,j)+tauinv*vol*MLocal(i,j))*
            OldSol(cells(icell,j))+vol*MLocal(i,j)*source(cells(icell,j));
        }    
    }

    // Assemble boundary conditions (Dirichlet penalty method)
    int nbfaces=grid.nbfaces();
    auto bfaces=grid.get_bfaces();
    auto bfaceregions=grid.get_bfaceregions();
    
    for (int ibface=0; ibface<nbfaces; ibface++)
    {
      // Obtain number of boundary condition
      int ireg=bfaceregions(ibface);
      int i0=bfaces(ibface,0);
      int i1=bfaces(ibface,1);
      
      // Check if it is "Dirichlet"
      if (bcfac(ireg)>=Dirichlet)
      {
        // Assemble penalty values
        SGlobal(i0,i0)+=bcfac(ireg);
        Rhs(i0)+=bcfac(ireg)*bcval(ireg);
        
        SGlobal(i1,i1)+=bcfac(ireg);
        Rhs(i1)+=bcfac(ireg)*bcval(ireg);
      }
      else if (bcfac(ireg)>0.0)
      {

        double dx= points(bfaces(ibface,1),0)-points(bfaces(ibface,0),0);
        double dy= points(bfaces(ibface,1),1)-points(bfaces(ibface,0),1);
        double h=sqrt(dx*dx+dy*dy);
        // first oder quadrature, needs to be replaced by second order
        int i0=bfaces(ibface,0);
        int i1=bfaces(ibface,1);
        
        SGlobal(i0,i0)+=theta*h*bcfac(ireg)/3.0;
        SGlobal(i0,i1)+=theta*h*bcfac(ireg)/6.0;
        
        SGlobal(i1,i0)+=theta*h*bcfac(ireg)/6.0;
        SGlobal(i1,i1)+=theta*h*bcfac(ireg)/3.0;
        
        Rhs(i0)+=OldSol(i0)*(1.0-theta)*h*bcfac(ireg)/3.0;
        Rhs(i0)+=OldSol(i1)*(1.0-theta)*h*bcfac(ireg)/6.0;
        Rhs(i1)+=OldSol(i0)*(1.0-theta)*h*bcfac(ireg)/6.0;
        Rhs(i1)+=OldSol(i1)*(1.0-theta)*h*bcfac(ireg)/3.0;
        
        Rhs(i0)+=0.5*h*bcval(ireg);
        Rhs(i0)+=0.5*h*bcval(ireg);
      }
    }
    SGlobal.flush();
  }

Here is the call graph for this function:

double fem2d::l2norm	(	const numcxx::SimpleGrid &	grid,
		const numcxx::DArray1 &	u
	)

Examples:

42-convtest-fem-sin.cxx.

Definition at line 275 of file fem2d.cxx.

  {
    auto ndim=grid.spacedim();
    auto points=grid.get_points(); // Array of global nodes
    auto cells=grid.get_cells();   // Local-global dof map
    int npoints=grid.npoints();
    int ncells=grid.ncells();
    
    // Local mass matrix
    numcxx::DArray2 MLocal{{2,1,1},{1,2,1},{1,1,2}}; // from Stroud quadrature...
    MLocal*=1.0/12.0;
    
    double norm=0.0;
    for (int icell=0; icell<ncells; icell++)
    {
      double vol;
      compute_cell_volume(icell,points,cells,vol);
      
      for (int i=0;i<=ndim;i++)
        for (int j=0;j<=ndim;j++)
          norm+=u(cells(icell,i))*u(cells(icell,j))*vol*MLocal(i,j);
    }
    return sqrt(norm);
  }

Here is the call graph for this function:

double fem2d::h1norm	(	const numcxx::SimpleGrid &	grid,
		const numcxx::DArray1 &	u
	)

Examples:

42-convtest-fem-sin.cxx.

Definition at line 302 of file fem2d.cxx.

  {
    auto ndim=grid.spacedim();
    auto points=grid.get_points(); // Array of global nodes
    auto cells=grid.get_cells();   // Local-global dof map
    int npoints=grid.npoints();
    int ncells=grid.ncells();


    numcxx::DArray2 SLocal(ndim+1,ndim+1);

    double norm=0.0;
    for (int icell=0; icell<ncells; icell++)
    {
      double vol;
      compute_local_stiffness_matrix(icell, points,cells, SLocal, vol);
      for (int i=0;i<=ndim;i++)
        for (int j=0;j<=ndim;j++)
          norm+=u(cells(icell,i))*u(cells(icell,j))*SLocal(i,j);
    }
    return sqrt(norm);
  }

Here is the call graph for this function:

void fem2d::compute_cell_volume ( const int  icell, const numcxx::DArray2 &  points, const numcxx::IArray2 &  cells, double &  vol  )

inline

Definition at line 11 of file fem2d.cxx.

  {
    int i0=cells(icell,0);
    int i1=cells(icell,1);
    int i2=cells(icell,2);
    
    // Fill matrix V
    double V00= points(i1,0)- points(i0,0);
    double V01= points(i2,0)- points(i0,0);
    
    double V10= points(i1,1)- points(i0,1);
    double V11= points(i2,1)- points(i0,1);
    
    // Compute determinant
    double det=V00*V11 - V01*V10;
    vol=0.5*det;
  }

void fem2d::compute_local_stiffness_matrix ( const int  icell, const numcxx::DArray2 &  points, const numcxx::IArray2 &  cells, numcxx::DArray2 &  SLocal, double &  vol  )

inline

Definition at line 33 of file fem2d.cxx.

  {
    
    int i0=cells(icell,0);
    int i1=cells(icell,1);
    int i2=cells(icell,2);
    
    // Fill matrix V
    double V00= points(i1,0)- points(i0,0);
    double V01= points(i2,0)- points(i0,0);
    
    double V10= points(i1,1)- points(i0,1);
    double V11= points(i2,1)- points(i0,1);
    
    // Compute determinant
    double det=V00*V11 - V01*V10;
    double fac = 0.5/det;
    
    // Compute entries of local stiffness matrix
    SLocal(0,0)= fac * (  ( V10-V11 )*( V10-V11 )+( V01-V00 )*( V01-V00 ) );
    SLocal(0,1)= fac * (  ( V10-V11 )* V11          - ( V01-V00 )*V01 );
    SLocal(0,2)= fac * ( -( V10-V11 )* V10          + ( V01-V00 )*V00 );
    
    SLocal(1,1)=  fac * (  V11*V11 + V01*V01 );
    SLocal(1,2)=  fac * ( -V11*V10 - V01*V00 );
    
    SLocal(2,2)=  fac * ( V10*V10+ V00*V00 );
    
    SLocal(1,0)=SLocal(0,1);
    SLocal(2,0)=SLocal(0,2);
    SLocal(2,1)=SLocal(1,2);
    
    vol=0.5*det;
  }

Variable Documentation

const double fem2d::Dirichlet =1.0e30

BC value marking Dirichlet boundary condition.

Examples:

40-stationary-heat-fe.cxx, 42-convtest-fem-sin.cxx, and 44-transient-heat-fe.cxx.

Definition at line 11 of file fem2d.hxx.