High Order Spaces¶
Creating the Element¶
For higher order spaces most of the times some type of continuity of the bases functions over edges or faces must be fullfilled. Because of this the order of the basis functions, for example on an edge is important. So the edges and faces must be oriented consistently. For this NGSolve provides the interface class VertexOrientedFE<ELEMENT_TYPE> defined in fem/fe_interfaces.hpp.
This class stores the vertex numbers of the element and when asked for the vertices of an edge it returns them sorted according to the global orientation. With this information we can build basis functions which are consistent over element boundaries. The definition of the segment looks similar to the definition of low order elements:
class MyHighOrderSegm : public MyBaseElement, public VertexOrientedFE<ET_SEGM>
{
public:
using VertexOrientedFE<ET_SEGM>::SetVertexNumbers;
MyHighOrderSegm (int order);
virtual ELEMENT_TYPE ElementType() const { return ET_SEGM; }
virtual void CalcShape (const IntegrationPoint & ip,
BareSliceVector<> shape) const;
virtual void CalcDShape (const IntegrationPoint & ip,
BareSliceMatrix<> dshape) const;
private:
template <class T>
void T_CalcShape (const T& x, BareSliceVector<T> shape) const;
};
We will use the AutoDiff class from the implementation of the second order elements again to not have to compute the gradients of the shape function manually. For this we define the template function T_CalcShape which can either take an arbitrary type (we will use double and AutoDiff) and compute the value of the basis function.
For the segment the first two basis functions are the vertex ("hat") functions, then the higher order edge basis functions (integrated legendre polynomials):
template <class T>
void MyHighOrderSegm :: T_CalcShape (const T& x, BareSliceVector<T> shape) const
{
T lami[2] = { x, 1-x };
for (int i = 0; i < 2; i++)
shape[i] = lami[i];
int ii = 2;
ArrayMem<T, 20> polx(order+1);
if (order >= 2)
{
auto edge = GetVertexOrientedEdge(0);
IntegratedLegendrePolynomial (order,
lami[edge[1]]-lami[edge[0]],
polx);
for (int j = 2; j <= order; j++)
shape[ii++] = polx[j];
}
}
CalcShape and CalcDShape just use this helper function:
void MyHighOrderSegm :: CalcShape (const IntegrationPoint & ip,
BareSliceVector<> shape) const
{
double x = ip(0);
T_CalcShape (x, shape);
}
void MyHighOrderSegm :: CalcDShape (const IntegrationPoint & ip,
BareSliceMatrix<> dshape) const
{
AutoDiff<1> adx (ip(0), 0);
Vector<AutoDiff<1> > shapearray(ndof);
T_CalcShape<AutoDiff<1>> (adx, shapearray);
for (int i = 0; i < ndof; i++)
dshape(i, 0) = shapearray[i].DValue(0);
}
The implementation of trigs is similar, only inner dofs (bubbles) are added:
template <class T>
void MyHighOrderTrig :: T_CalcShape (const T & x, const T & y, BareSliceVector<T> shape) const
{
T lami[3] = { x, y, 1-x-y };
for (int i = 0; i < 3; i++)
shape[i] = lami[i];
int ii = 3;
ArrayMem<T, 20> polx(order+1), poly(order+1);
for (int i = 0; i < 3; i++)
if (order >= 2) // more general: order_edge[i]
{
auto edge = GetVertexOrientedEdge(i);
ScaledIntegratedLegendrePolynomial (order,
lami[edge[1]]-lami[edge[0]],
lami[edge[0]]+lami[edge[1]], polx);
for (int j = 2; j <= order; j++)
shape[ii++] = polx[j];
}
// inner dofs
if (order >= 3) // more general: cell order
{
T bub = x * y * (1-x-y);
ScaledLegendrePolynomial (order-2, lami[1]-lami[0], lami[1]+lami[0], polx);
LegendrePolynomial (order-1, 2*lami[2]-1, poly);
for (int i = 0; i <= order-3; i++)
for (int j = 0; j <= order-3-i; j++)
shape[ii++] = bub * polx[i] * poly[j];
}
}
}
Creating the FESpace¶
The definition of the finite element space looks similar to the low order case, we only store an array of the dof numbers of edge and cell dofs:
class MyHighOrderFESpace : public FESpace
{
int order;
Array<int> first_edge_dof;
Array<int> first_cell_dof;
public:
/*
constructor.
Arguments are the access to the mesh data structure,
and the flags from the define command in the pde-file
*/
MyHighOrderFESpace (shared_ptr<MeshAccess> ama, const Flags & flags);
// a name for our new fe-space
string GetClassName () const override
{
return "MyHighOrderFESpace";
}
void Update() override;
void GetDofNrs (ElementId ei, Array<DofId> & dnums) const override;
FiniteElement & GetFE (ElementId ei, Allocator & alloc) const override;
};
In the Update function we fill this array and compute some global variables:
void MyHighOrderFESpace :: Update()
{
// some global update:
int n_vert = ma->GetNV();
int n_edge = ma->GetNEdges();
int n_cell = ma->GetNE();
first_edge_dof.SetSize (n_edge+1);
int ii = n_vert;
for (int i = 0; i < n_edge; i++, ii+=order-1)
first_edge_dof[i] = ii;
first_edge_dof[n_edge] = ii;
first_cell_dof.SetSize (n_cell+1);
for (int i = 0; i < n_cell; i++, ii+=(order-1)*(order-2)/2)
first_cell_dof[i] = ii;
first_cell_dof[n_cell] = ii;
// cout << "first_edge_dof = " << endl << first_edge_dof << endl;
// cout << "first_cell_dof = " << endl << first_cell_dof << endl;
SetNDof (ii);
}
The GetDofNrs function collects the dof nrs into the output array:
void MyHighOrderFESpace :: GetDofNrs (ElementId ei, Array<DofId> & dnums) const
{
// returns dofs of element number elnr
dnums.SetSize(0);
auto ngel = ma->GetElement (ei);
// vertex dofs
for (auto v : ngel.Vertices())
dnums.Append(v);
// edge dofs
for (auto e : ngel.Edges())
for(auto j : Range(first_edge_dof[e], first_edge_dof[e+1]))
dnums.Append (j);
// inner dofs
if (ei.IsVolume())
for(auto j : Range(first_cell_dof[ei.Nr()], first_cell_dof[ei.Nr()+1]))
dnums.Append (j);
}
When returning the finite element we need to set the vertex numbers, so that VertexOrientedFE has the correct orientation:
FiniteElement & MyHighOrderFESpace :: GetFE (ElementId ei, Allocator & alloc) const
{
auto ngel = ma->GetElement (ei);
switch (ngel.GetType())
{
case ET_TRIG:
{
auto trig = new (alloc) MyHighOrderTrig(order);
trig->SetVertexNumbers (ngel.vertices);
return *trig;
}
case ET_SEGM:
{
auto segm = new (alloc) MyHighOrderSegm(order);
segm->SetVertexNumbers (ngel.vertices);
return *segm;
}
default:
throw Exception (string("Element type ")+ToString(ngel.GetType())+" not supported");
}
}
Finally we have to register the new space to the NGSolve register:
static RegisterFESpace<MyHighOrderFESpace> initifes ("myhofespace");