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";
}
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");
}
void ExportMyFESpace(py::module m)
{
using namespace ngcomp;
/*
We just export the class here and use the FESpace constructor to create our space.
This has the advantage, that we do not need to specify all the flags to parse (like
dirichlet, definedon,...), but we can still append new functions only for that space.
*/
auto myfes = py::class_<MyFESpace, shared_ptr<MyFESpace>, FESpace>
(m, "MyFESpace", "FESpace with first order and second order trigs on 2d mesh");
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);
auto pyfes = py::cast(fes);
pyfes.attr("__initialize__")(**kwa);
return fes;
}), py::arg("mesh"))
/*
this is, so that we do not get an 'undocumented flag' warning
*/
.def_static("__flags_doc__", [] ()
{
auto doc = py::cast<py::dict>(py::module::import("ngsolve").
attr("FESpace").
attr("__flags_doc__")());
doc["secondorder"] = "bool = False \n"
" Use second order elements.";
return doc;
})
.def("GetNVert", &MyFESpace::GetNVert)
;
}
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", "FESpace with first order and second order trigs on 2d mesh");
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);
auto pyfes = py::cast(fes);
pyfes.attr("__initialize__")(**kwa);
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__", [] ()
{
auto doc = py::cast<py::dict>(py::module::import("ngsolve").
attr("FESpace").
attr("__flags_doc__")());
doc["secondorder"] = "bool = False \n"
" Use second order elements.";
return doc;
})
Here we query the flags_doc from the base class and append our secondorder flag to the dictionary.