My Element

We want to start with a simple linear triangular finite element for the space \(H^1\). Therefore we can use the ScalarFiniteElement as a base class:

  class MyLinearTrig : public ScalarFiniteElement<2>
  {
  public:
    // constructor
    MyLinearTrig ();

    virtual ELEMENT_TYPE ElementType() const { return ET_TRIG; }
    
    /*
      Calculate the vector of shape functions in the point ip.
      ip is given in the reference element.
     */
    virtual void CalcShape (const IntegrationPoint & ip, 
                            BareSliceVector<> shape) const;
  
    /*
      Calculate the matrix of derivatives of the shape functions in the point ip.
      dshape is an 3 by 2 matrix in our case.
     */
    virtual void CalcDShape (const IntegrationPoint & ip, 
                             BareSliceMatrix<> dshape) const;

    // there are some more functions to bring in ...
    using ScalarFiniteElement<2>::CalcShape;    
    using ScalarFiniteElement<2>::CalcDShape;
  };

So we need to define 3 functions: The constructor, CalcShape and CalcDShape to compute the shape functions and their derivatives.

  MyLinearTrig :: MyLinearTrig ()
  /*
    Call constructor for base class:
    geometry is ET_TRIG, number of dofs is 3, maximal order is 1
   */
    : ScalarFiniteElement<2> (3, 1)
  { ; }


  void MyLinearTrig :: CalcShape (const IntegrationPoint & ip,
                                  BareSliceVector<> shape) const
  {
    // coordinates in reference elements
    double x = ip(0);
    double y = ip(1);

    /*
      Vertex coordinates have been defined to be (1,0), (0,1), (0,0)
      see ElementTopology::GetVertices(ET_TRIG)
     */

    // define shape functions
    shape(0) = x;
    shape(1) = y;
    shape(2) = 1-x-y;
  }

  void MyLinearTrig :: CalcDShape (const IntegrationPoint & ip,
                                   BareSliceMatrix<> dshape) const

  {
    // matrix of derivatives:

    dshape(0,0) = 1;
    dshape(0,1) = 0;
    dshape(1,0) = 0;
    dshape(1,1) = 1;
    dshape(2,0) = -1;
    dshape(2,1) = -1;
  }

Finally we can export it to Python, as we have done it with the CoefficientFunction. Note that this is not necessary if we just want to build a finite element space with it.

void ExportMyElement(py::module m)
{
  using namespace ngfem;
  /*
    Our Trig is derived
    from the classes ScalarFiniteElement<2> -> BaseScalarFiniteElement -> FiniteElement.
    Only BaseScalarFiniteElement and FiniteElement are exported to python
    (see ngsolve/fem/python_fem.cpp), so we derive from BaseScalarFiniteElement (which derives from
    FiniteElement).
    If we only want to use it in our FESpace we do not need the Python hierarchy, but it's nice for
    debugging :)
  */
  py::class_<MyLinearTrig, shared_ptr<MyLinearTrig>, BaseScalarFiniteElement>
    (m, "MyLinearTrig", "My new linear Trig")
    .def(py::init<>())
    ;
}

Todo

Explain AutoDiff with second order elements

Next we want to build this finite element into a space, so we can use it.