Custom Integrators

Sometimes new DifferentialOperators and CoefficientFunctions are not enough and one needs to implement special element matrices or vectors. This can be done by creating BilinearFormIntegrators and LinearFormIntegrators.

The integrators need to implement some informations about what dimension it is defined on, if it is a boundary or volume integrator and so on.

  class MyLaplaceIntegrator : public BilinearFormIntegrator
  {
    shared_ptr<CoefficientFunction> coef_lambda;
  public:
    MyLaplaceIntegrator (shared_ptr<CoefficientFunction> coef)
      : coef_lambda(coef) { ; }

    string Name () const override { return "MyLaplace"; }

    int DimElement () const override { return 2; }
    int DimSpace () const override { return 2; }
    xbool IsSymmetric () const override { return true; }
    VorB VB() const override { return VOL; }

    // Calculates the element matrix
    void CalcElementMatrix (const FiniteElement & fel,
                            const ElementTransformation & eltrans, 
                            FlatMatrix<double> elmat,
                            LocalHeap & lh) const override;
  };

For a BilinearFormIntegrator the main function to be implemented is the construction of element matrices. This is done by creating an integration rule and iterating over it using the CalcDShape function of our elements.

  void MyLaplaceIntegrator ::
  CalcElementMatrix (const FiniteElement & base_fel,
		     const ElementTransformation & eltrans, 
		     FlatMatrix<double> elmat,
		     LocalHeap & lh) const
  {
    static Timer t("MyCalcElementMatrix");
    RegionTracer reg(TaskManager::GetThreadId(), t);
    
    /*
      tell the compiler that we are expecting one of our elements.
      if not, an exception will be raised
    */
    auto & fel = dynamic_cast<const MyBaseElement &> (base_fel);

    // number of element basis functions:
    int ndof = fel.GetNDof();

    elmat = 0;

    Matrix<> dshape_ref(ndof, 2); // gradient on reference element
    Matrix<> dshape(ndof, 2);     // gradient on mapped element

    /*
      get integration rule for element geometry, 
      integration order is 2 times element order
    */
    IntegrationRule ir(fel.ElementType(), 2*fel.Order());

    // loop over integration points
    for (int i = 0 ; i < ir.GetNIP(); i++)
      {
        // calculate Jacobi matrix in the integration point
        MappedIntegrationPoint<2,2> mip(ir[i], eltrans);

        // lambda(x)
        double lam = coef_lambda -> Evaluate (mip);

        /*
          gradient on reference element
          the i-th row of the matrix is the grad of the i-th shape function
        */
        fel.CalcDShape (ir[i], dshape_ref);
        
        // transform it for the mapped element
        dshape = dshape_ref * mip.GetJacobianInverse();
        
        // integration weight and Jacobi determinant
        double fac = mip.IP().Weight() * mip.GetMeasure();

        // elmat_{i,j} += (fac*lam) * InnerProduct (grad shape_i, grad shape_j)
        elmat += (fac*lam) * dshape * Trans(dshape);
      }     
  }

  /*

Similarily a LinearFormIntegrator has to implement the construction of element vectors:

  void MySourceIntegrator ::
  CalcElementVector (const FiniteElement & base_fel,
		     const ElementTransformation & eltrans, 
		     FlatVector<double> elvec,
		     LocalHeap & lh) const
  {
    auto & fel = dynamic_cast<const MyBaseElement&> (base_fel);

    int ndof = fel.GetNDof();

    elvec = 0;

    Vector<> shape(ndof); 

    IntegrationRule ir(fel.ElementType(), 2*fel.Order());

    for (int i = 0 ; i < ir.GetNIP(); i++)
      {
        MappedIntegrationPoint<2,2> mip(ir[i], eltrans);

        double f = coef_f -> Evaluate (mip);

        // calculate shape functions 
        fel.CalcShape (ir[i], shape);
        
        // integration weight and Jacobi determinant
        double fac = mip.IP().Weight() * mip.GetMeasure();

        // elvec_{i} += (fac*f) shape_i
        elvec += (fac*f) * shape;
      }     
  }
}

Again we export the integrators to Python:

void ExportMyIntegrator(py::module m)
{
  using namespace ngfem;
  py::class_<MyLaplaceIntegrator, shared_ptr<MyLaplaceIntegrator>, BilinearFormIntegrator>
    (m, "MyLaplace")
    .def(py::init<shared_ptr<CoefficientFunction>>())
    ;
  py::class_<MySourceIntegrator, shared_ptr<MySourceIntegrator>, LinearFormIntegrator>
    (m, "MySource")
    .def(py::init<shared_ptr<CoefficientFunction>>())
    ;
}