# My DiffOpΒΆ

DiffOps are the static polymorphism classes implementing Operators for GridFunctions and Test- and TrialFunctions. They only have static member functions and are never created but only used as a template argument for DifferentialOperators. They use CRTP for "static virtual" methods. DifferentialOperator is the runtime polymorphism class wrapping them. Use T_DifferentialOperator<DiffOp> as a coupling mechanism.

We will need the identity and the gradient operator. These operators have to define the following constexpr values:

 DIM Dimension of input DIM_SPACE Dimension of FESpace DIM_ELEMENT Dimension of Element DIM_DMAT Dimension of range DIFFORDER Order of differentiation

The order of differentiation is used to choose the correct integration rules for integrating the polynomials. Further we need to implement a function that generates the result of the Diffop for one mapped integration point. This function gets a matrix because the result could be matrix valued. A minimal implementation of the identity would be this:

```  class MyDiffOpId : public DiffOp<MyDiffOpId>
{
public:
// ******************** Necessary stuff *************************

static constexpr int DIM = 1; // dimension of the input
static constexpr int DIM_SPACE = 2; // dimension of the space
static constexpr int DIM_ELEMENT = 2; // dimension of the element
static constexpr int DIM_DMAT = 1; // dimension of the output
static constexpr int DIFFORDER = 0; // order of differentiation

template<typename MIP, typename MAT>
static void GenerateMatrix(const FiniteElement& fel, const MIP& mip,
MAT& mat, LocalHeap& lh)
{
HeapReset hr(lh);
Cast(fel).CalcShape(mip.IP(), mat.Row(0));
}

// ******************** Helper function ***************************
static const MyBaseElement& Cast(const FiniteElement& fel)
{ return static_cast<const MyBaseElement&> (fel); }

};
}
```

Again some more functions can be implemented to improve performance, for example evaluation in the whole integration rule and SIMD implementations.

```  class MyDiffOpGradient : public DiffOp<MyDiffOpGradient>
{
public:
// ******************** Necessary stuff *************************

static constexpr int DIM = 1; // dimension of the input
static constexpr int DIM_SPACE = 2; // dimension of the space
static constexpr int DIM_ELEMENT = 2; // dimension of the element
static constexpr int DIM_DMAT = 2; // dimension of the output
static constexpr int DIFFORDER = 1; // order of differentiation

// so that you can call grad(u)
static string Name() { return "grad"; }

template<typename MIP, typename MAT>
static void GenerateMatrix(const FiniteElement& fel, const MIP& mip,
MAT& mat, LocalHeap& lh)
{
// We need to compute the gradient on the reference element and
// then transform it to the using the jacobian of the mapped
// integration point
HeapReset hr(lh);
// create a temporary matrix on the local heap
FlatMatrixFixWidth<2> dshape(fel.GetNDof(), lh);
Cast(fel).CalcDShape(mip.IP(), dshape);
mat = Trans(dshape * mip.GetJacobianInverse());
}

// ******************** Helper function ***************************
static const MyBaseElement& Cast(const FiniteElement& fel)
{ return static_cast<const MyBaseElement&> (fel); }

};
}
```

Again the implementation of a boundary evaluation operator will be left as an exercise.

To make the DiffOp available at runtime we have to put it into a T_DifferentialOperator class as template argument. This happens next in the FESpace implementation.