NGSolve  5.3
Some representative functions

This is an overview of some representative functions of NGSolve.

For instruction purposes, the functions are simplified.

Element matrix assembling for BDB integrators:

Many bilinear-forms are of the structure

$\int D \partial u \cdot \partial v \, dx $

where $\partial$ is some differential operator, i.e., the gradient, and $D$ is the material tensor. Then, the finite elemen matrix $A_T$ is computed as

\[ A_T = \sum_{Int.Pts.} \omega |J| B^T D B. \]

The elements of the matrix $B$ are the values of the differential operator applied to the shape functions in the integration point, i.e., $B_{i,j} = (\partial \varphi_j)_i$.

The BDBIntegrator is a derived class from BilinearFormIntegrator. Here, the differential operator and the material tensor are defined by template arguments.

  template <class DIFFOP, class DMATOP>
  void BDBIntegrator<DIFFOP,DMATOP> :: 
  AssembleElementMatrix (const FiniteElement & fel, 
                         const ElementTransformation & eltrans, 
                         Matrix & elmat,
                         LocalHeap & locheap) const
  {
    enum { DIM_SPACE   = DIFFOP::DIM_SPACE };
    enum { DIM_ELEMENT = DIFFOP::DIM_ELEMENT };
    enum { DIM_DMAT    = DIFFOP::DIM_DMAT };
    enum { DIM         = DIFFOP::DIM };

    int ndof = fel.GetNDof();
   
    elmat.AssignMemory (ndof*DIM, ndof*DIM, locheap);
    elmat = 0;
        
    MatrixFixHeight<DIM_DMAT> bmat (ndof * DIM, locheap);
    MatrixFixHeight<DIM_DMAT> dbmat (ndof * DIM, locheap);
    Mat<DIM_DMAT,DIM_DMAT> dmat;

    const IntegrationRule & ir = GetIntegrationRule (fel);

    // optimized memory management
    void * heapp = locheap.GetPointer();

    for (int i = 0; i < ir.GetNIP(); i++)
      {
        // the mapped point, including Jacobian
        SpecificIntegrationPoint<DIM_ELEMENT,DIM_SPACE> 
          sip(ir, i, eltrans, locheap);

        DIFFOP::GenerateMatrix (fel, sip, bmat, locheap);
        dmatop.GenerateMatrix (fel, sip, dmat, locheap);

        double fac = fabs (sip.GetJacobiDet()) * sip.IP().Weight();

        dbmat = dmat * bmat;
        elmat += fac * (Trans (bmat) * dbmat);

        locheap.CleanUp (heapp);
      } 
    }

Global matrix assembling The bilinear-form consists of a couple of integrators acting on the domain or on the boundary. It has a reference to the finite element space it is defined on. The sparse matrices are stored in the bilinear-form object.

template <typename SCAL = double>
void BilinearForm :: Assemble ()
{
  Array<int> dnums;
  Matrix<SCAL> Matrix elmat;
  ElementTransformation eltrans;
      
  int ndof = fespace.GetNDof();
      
  BaseMatrix & mat = GetMatrix();
  mat = 0.0;
      
  LocalHeap locheap (1000000);

  // assembling of volume terms
  int ne = ma.GetNE();
  for (int i = 0; i < ne; i++)
    {
      locheap.CleanUp();
      if (!fespace.DefinedOn (ma.GetElIndex (i))) continue;

      ma.GetElementTransformation (i, eltrans);
      const FiniteElement & fel = fespace.GetFE (i);
      fespace.GetDofNrs (i, dnums);

      for (int j = 0; j < parts.Size(); j++)
        {
          const BilinearFormIntegrator & bfi = *parts.Get(j);
          if (bfi.BoundaryForm()) continue;
              
          bfi.AssembleElementMatrix (fel, eltrans, elmat, &locheap);
          fespace.TransformMatrix (i, elmat);
          mat->AddElementMatrix (dnums, elmat);
        }
    }
       
  // assembling of surface terms
  int nse = ma.GetNSE();
  for (int i = 0; i < nse; i++)
    {
      locheap.CleanUp();
      if (!fespace.DefinedOnBoundary (ma.GetSElIndex (i))) continue;

      ma.GetSurfaceElementTransformation (i, eltrans);
      const FiniteElement & fel = fespace.GetSFE (i);
      fespace.GetSDofNrs (i, dnums);

      for (int j = 1; j <= parts.Size(); j++)
        {
          const BilinearFormIntegrator & bfi = *parts.Get(j);
          if (!bfi.BoundaryForm()) continue;
              
          bfi.AssembleElementMatrix (fel, eltrans, elmat, &locheap);
          fespace.TransformSurfMatrix (i, elmat);
          mat->AddElementMatrix (dnums, elmat);
        }
    }
}