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;
    int ndof, 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
    virtual string GetClassName () const
    {
      return "MyFESpace";
    }

    static DocInfo GetDocu();

    virtual void Update(LocalHeap & lh);
    virtual size_t GetNDof () const { return ndof; }
    
    virtual void GetDofNrs (ElementId ei, Array<DofId> & dnums) const;
    virtual FiniteElement & GetFE (ElementId ei, Allocator & alloc) const;

    // some new functionality our space should have in Python
    int 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:

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

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

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

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

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

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

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

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
      {
        if (!secondorder)
          return * new (alloc) FE_Segm1;
        else
          return * new (alloc) FE_Segm2;
      }
  }

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

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

Registering you FESpace for Python is quite similar as before, only the init (optional) function is different:

  auto myfes = py::class_<MyFESpace, shared_ptr<MyFESpace>, FESpace>
    (m, "MyFESpace", (docu.short_docu + "\n\n" + docu.long_docu).c_str());
  myfes
    /*
       this is optional, if you don't write an init function, you can create your fespace
       with FESpace("myfes",mesh,...), but it's nicer to write MyFESpace(mesh,...) ;)
    */
    .def(py::init([myfes] (shared_ptr<MeshAccess> ma, py::kwargs kwa)
                  {
                    py::list info;
                    info.append(ma);
                    auto flags = CreateFlagsFromKwArgs(myfes, kwa, info);
                    auto fes = make_shared<MyFESpace>(ma, flags);
                    LocalHeap glh(100000000, "init-fes-lh");
                    fes->Update(glh);
                    fes->FinalizeUpdate(glh);
                    return fes;
                  }), py::arg("mesh"))

The init function takes the mesh as argument and all the flags as keyword arguments, then we call the CreateFlagsFromKwArgs function with the mesh as additional information (you can just copy paste this from the base class constructor). This is necessary, because the base FESpace may have some flags we don't know anything about (and do not have to care about) in the derived class, i.e. dirichlet boundaries or it is defined on only a part of the domain. Then we create our new space and call the __initialize__ function from the base class, which calls amongst other things the Update function of our space. Instead of this init function you can create the registered FESpace in Python as well just with

FESpace("myfespace", mesh, ...)

with the name we registered it.

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, so we want to register our flag for the space:

    .def_static("__flags_doc__", [docu] ()
                {
                  auto doc = py::cast<py::dict>(py::module::import("ngsolve").
                                                attr("FESpace").
                                                attr("__flags_doc__")());
                  for (auto & flagdoc : docu.arguments)
                    doc[get<0>(flagdoc).c_str()] = get<1>(flagdoc);
                  return doc;
                })

Here we query the flags_doc from the base class and append our secondorder flag to the dictionary.