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.

The gradient operator looks similar:

  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.