# 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");

/*
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);

/*
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 += (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>>())
;
}
```