My Element

We want to start with a simple linear triangular finite element for the space \(H^1\). First we create a base class, so that we can later create higher order elements as well. This base class requires elements to be able to compute values and gradients:

  class MyBaseElement : public FiniteElement
  {
  public:
    MyBaseElement(int ndof, int order) : FiniteElement(ndof, order) {}

    virtual void CalcShape(const IntegrationPoint& ip,
                           BareSliceVector<> shape) const = 0;
    virtual void CalcDShape(const IntegrationPoint& ip,
                            BareSliceMatrix<> dshape) const = 0;
  };

The linear trig implements this functionality:

  class MyLinearTrig : public MyBaseElement
  {
  public:
    MyLinearTrig () : MyBaseElement(3, 1) {}
    ELEMENT_TYPE ElementType() const override { return ET_TRIG; }

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

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

  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;
  }

Automatic Differentiation

Next we will create second order elements for our space. When creating higher order spaces the derivatives of the polynomials will get quite ugly, so we don't want to implement them by hand. NGSolve provides automatic exact differentiation by the AutoDiff class. A AutoDiff object knows its value and the value of its derivative. Combining AutoDiff objects automatically uses differentiation rules.

We implement nodal quadratic elements, note CalcShape and CalcDShape look quite similar, in CalcDShape we only use the diff value of the Autodiff variable:

  void MyQuadraticTrig :: CalcShape (const IntegrationPoint & ip,
                                     BareSliceVector<> shape) const
  {
    // now, use barycentric coordinates x, y, 1-x-y:
    double lam[3] = { ip(0), ip(1), 1-ip(0)-ip(1) };

    // vertex basis functions:
    for (int i = 0; i < 3; i++)
      shape(i) = lam[i] * (2 * lam[i] - 1);


    // edge basis functions:

    const EDGE * edges = ElementTopology::GetEdges (ET_TRIG);
    // table provides connection of edges and vertices
    // i-th edge is between vertex edges[i][0] and edges[i][1]

    for (int i = 0; i < 3; i++)
      shape(3+i) = 4 * lam[edges[i][0]] * lam[edges[i][1]];
  }



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

  {
    // Use automatic (exact !) differentiation with overloaded data-types

    AutoDiff<2> x (ip(0), 0); // value of x, gradient is 0-th unit vector (1,0)
    AutoDiff<2> y (ip(1), 1); // value of y, gradient is 1-th unit vector (0,1)


    AutoDiff<2> lam[3] = { x, y, 1-x-y };

    // vertex basis functions:
    for (int i = 0; i < 3; i++)
      {
        AutoDiff<2> shape = lam[i] * (2 * lam[i] - 1);
        dshape(i,0) = shape.DValue(0);    // x-derivative
        dshape(i,1) = shape.DValue(1);    // y-derivative
      }


    // edge basis functions:
    const EDGE * edges = ElementTopology::GetEdges (ET_TRIG);

    for (int i = 0; i < 3; i++)
      {
        AutoDiff<2> shape = 4 * lam[edges[i][0]] * lam[edges[i][1]];
        dshape(3+i,0) = shape.DValue(0);    // x-derivative
        dshape(3+i,1) = shape.DValue(1);    // y-derivative
      }
  }
}

The creation of boundary elements will be left as an exercise. Next we will create (differential) operators working with these finite elements.