My FESpace

Next we want to use our created finite elements and put them into a space. While the finite element handles all the local element wise computations, the space knows the global mapping of the dofs, has access to the mesh and connects local to global computations.

Again there are only a few virtual function overrides necessary to get our space working:

  class MyFESpace : public FESpace
  {
    bool secondorder;
    size_t nvert;
    
  public:
    /*
      constructor. 
      Arguments are the access to the mesh data structure,
      and the flags from the define command in the pde-file
      or the kwargs in the Python constructor.
    */
    MyFESpace (shared_ptr<MeshAccess> ama, const Flags & flags);

    // a name for our new fe-space
    string GetClassName () const override { return "MyFESpace"; }

    static DocInfo GetDocu();

    void Update() override;
    
    void GetDofNrs (ElementId ei, Array<DofId> & dnums) const override;
    FiniteElement & GetFE (ElementId ei, Allocator & alloc) const override;

    // some new functionality our space should have in Python
    size_t GetNVert() { return nvert; }
  };

In the constructor we query the flags (if we want a second order space) and define evaluators for functions on the space, their derivatives and their boundary values (which will be left as an exercise):

  MyFESpace :: MyFESpace (shared_ptr<MeshAccess> ama, const Flags & flags)
    : FESpace (ama, flags)
  {
    cout << "Constructor of MyFESpace" << endl;
    cout << "Flags = " << flags << endl;

    // this is needed for pickling and needs to be the same as used in
    // RegisterFESpace later
    type = "myfespace";

    secondorder = flags.GetDefineFlag ("secondorder");

    if (!secondorder)
      cout << "You have chosen first order elements" << endl;
    else
      cout << "You have chosen second order elements" << endl;

    // needed for symbolic integrators and to draw solution
    evaluator[VOL] = make_shared<T_DifferentialOperator<MyDiffOpId>>();
    flux_evaluator[VOL] = make_shared<T_DifferentialOperator<MyDiffOpGradient>>();
  }

In the Update function, the number of degrees of freedom are set:

  void MyFESpace :: Update()
  {
    // some global update:
    cout << "Update MyFESpace, #vert = " << ma->GetNV()
         << ", #edge = " << ma->GetNEdges() << endl;

    // number of vertices
    nvert = ma->GetNV();

    // number of dofs:
    size_t ndof = nvert;
    if (secondorder)
      ndof += ma->GetNEdges();  // num vertics + num edges

    SetNDof (ndof);
  }

Further we must define the mapping from local dofs for each ElementId to the global dof numbers:

  void MyFESpace :: GetDofNrs (ElementId ei, Array<DofId> & dnums) const
  {
    // returns dofs of element ei
    // may be a volume triangle or boundary segment

    dnums.SetSize(0);

    // first dofs are vertex numbers:
    for (auto v : ma->GetElVertices(ei))
      dnums.Append (v);

    if (secondorder)
      {
        // more dofs on edges:
        for (auto e : ma->GetElEdges(ei))
          dnums.Append (nvert+e);
      }
  }

as well as a function that returns our new element:

  FiniteElement & MyFESpace :: GetFE (ElementId ei, Allocator & alloc) const
  {
    if (ei.IsVolume())
      {
        if (!secondorder)
          return * new (alloc) MyLinearTrig;
        else
          return * new (alloc) MyQuadraticTrig;
      }
    else
      throw Exception("Boundary elements not implemented yet!");
    // else
    //   {
    //     if (!secondorder)
    //       return * new (alloc) MyLinearSegm;
    //     else
    //       return * new (alloc) MyQuadraticSegm;
    //   }
  }

NGSolve has an internal register of finite element spaces, where we need to register our new space:

  static RegisterFESpace<MyFESpace> initifes ("myfespace");

To register FESpace for Python we provide the helper function ExportFESpace in python_comp.hpp. This automatically creates the constructor and pickling support for the space. The function returns the Python class and additional functionality can be appended.

  ExportFESpace<MyFESpace>(m, "MyFESpace")
    .def("GetNVert", &MyFESpace::GetNVert)
    ;

If we want to use the secondorder flag now, we would get a warning, that we use an undocumented flag. This is a safety guide to prevent typos for keyword arguments, what we need to do is to provide documentation for our flag:

  DocInfo MyFESpace :: GetDocu()
  {
    auto docu = FESpace::GetDocu();
    docu.Arg("secondorder") = "bool = False\n"
      "  Use second order basis functions";
    return docu;
  }

Next we will show how to implement custom integrators.