Approximation spaces and unknowns#

XLiFE++ allows to solve PDE with the Finite Elements Method, the Boundary Element Method, Discontinuous Galerkin Methods and the Spectral Elements Method or a coupling of them.

Given a basis of \(n\) functions \(\varphi_i(x, y, z)\), an approximation space of finite dimension is defined as:

\[V_h=\left\{w=\sum_{i=0}^n w_i \varphi_i(x, y, z)\right\}\]

XLiFE++ deals with two types of spaces : finite element spaces and spectral spaces. Finite element and spectral shape functions cannot be mixed in a same space. All spaces are managed by the Space top class which hides the specific classes FeSpace and SpSpace.

Attention

The essential conditions are taken into account algebraically in XLiFE++. They are not attached to the space, unlike in PDE variational formulations.

Finite elements spaces#

With the Finite Elements Method, the basis function \(\varphi_i\) is related to the \(i^{th}\) degree of freedom (DoF) localized in few elements, and constructed from local definition on these elements. It is, generally, named shape function. A finite element space is defined from

  • a mesh domain (Domain), generally the domain where the problem is set,

  • a finite element interpolation, such as \(P_k, k\in \mathbb{N}\) for instance.

Here are some examples of FE space construction, considering omega as a Domain:

Space V1(_domain=omega,_interpolation=P1,_name="V1");
Space V2(_domain=omega,_FE_type=Lagrange,_order=1,_name="V2"); // same as V1
Space V3(_domain=omega,_FE_type=Lagrange,_order=1,_Sobolev_type=H1,_name="V3"); // same as V1

Space V4(_domain=omega,_FE_type=RaviartThomas,_order=3,_Sobolev_type=Hdiv,_name="V4");
Space V5(_domain=omega,_interpolation=RT3,_name="V5"); // same as V4

Space V6(_domain=omega,_FE_type=Nedelec,_FE_subtype=firstFamily,_order=4,_Sobolev_type=Hcurl,_name="V6");
Space V7(_domain=omega,_interpolation=N1_4,_name="V7"); // same as V6

Space V8(_domain=omega,_FE_type=NedelecFace,_FE_subtype=firstFamily,_order=1,_soboloev_type=Hcurl,_name="V8");
Space V9(_domain=omega,_interpolation=NF1_1,_name="V9"); // same as V8

The definition of a finite element space is done with a key-value system:

  • _domain: the geometrical support of the space

  • _FE_type: the type of finite elements (Lagrange, Hermite, CrouzeixRaviart, Nedelec, RaviartThomas, NedelecEdge, NedelecFace, Morley)

  • _FE_subtype: the subtype of finite elements (standard, GaussLobattoPoints, firstFamily, secondFamily)

  • _Sobolev_type: the Sobolev space (L2, H1, Hdiv, Hcurl, H2, Hinf, Linf)

  • _order: the interpolation order

  • _interpolation: an other way to define the interpolation with shortcuts (P0 up to P10, P1Bubble, Q0 up to Q10, RT1 up to RT5, NF1_1 up to NF1_5, BDM_1 up to BDM_5, NF2_1 up to NF2_5, N1_1 up to N1_5, N2_1 up to N2_5). This is a shortcut to the set of keys {_FE_type, _FE_subtype, _Sobolev_type, _order}

  • _name: a name for print purpose (not optional)

See Finite elements for a description of finite elements.

When using the set of keys {_FE_type, _FE_subtype, _Sobolev_type, _order}, _FE_subtype and _Sobolev_type keys are optional; they take default values:

  • standard, H1 for Lagrange interpolation,

  • Hdiv, firstFamily for RaviartThomas and NedelecFace interpolation,

  • firstFamily, Hrot for NedelecEdge interpolation,

  • standard, H2 for Hermite and Morley interpolation.

These are also the default values when using the shortcuts _interpolation.

Note

To deal with FE discontinuous approximation, specify L2 as _Sobolev_type. In that case, dofs located on sides of elements are no longer shared.

When creating a space, all dofs are created according to chosen interpolation, and then comes the dof numbering. The defaut behaviour is that dof numbering is optimized to reduce the bandwidth of matrices. This can be controlled by using _optimizeNumbering or _notOptimizeNumbering keys.

When evaluate a function of the FE space at a point, XLiFE++ first searches for the location on this point in the mesh and then uses local interpolation. To make the location faster, you can activate precomputation of the search tree by using _withLocateData key. By default, it is not constructed (key _withoutLocateData), but will be automatically constructed at the first evaluation calling.

Warning

The second family edge or face elements are not yet available.

When dealing with problem with vector unknown where each component is approximated in the same space (for instance \(P1 \times P1 \times P1\) for the displacement field in elasticity problem), the product space is not built, but a ‘vector’ unknown on a ‘scalar’ space must be invoked; see Unknowns and test functions.

Some subspaces or trace spaces are automatically created by XLiFE++. For instance, when an integral on a boundary (\(\Sigma\)) is involved in a bilinear or linear form, the trace space \(V_{|\Sigma}=\{v_{|\Sigma},\ v\in V\}\) is created. This subspace can be got using the restriction operator:

Space& Vs=V|Sigma;  // do not forget the & (reference)

Spectral spaces#

Spectral spaces are spaces defined from basis functions defined on a mesh domain (Domain). Contrary to finite element basis function given by a local definition on elements, the spectral basis functions are given as global functions either by their analytic forms or by a set of interpolated functions (say vectors related to an other space).

The definition of spectral spaces is done with a key-value system:

  • _basis: the spectral basis as a Function (for analytical basis) or a TermVectors (for interpolated basis).

  • _dim: the spectral space dimension (for analytical basis).

  • _basis_dim: the basis function dimension: 1 for scalar function, 2 or 3 for vector function. Default value is 1.

  • _name: a name for print purpose (not optional).

Analytic spectral space#

Analytic spectral spaces are constructed from a standard C++ function taking a Point and a Parameters as argument and returning function values with the right type (Real, Complex,…) depending on the problem. It managed the basis index and , if necessary, the derivation operator. Let’s look at an example with the basis \(\varphi_n(x, y)=\sqrt{\frac{2}{h}}\sin \frac{n\pi}{h},\ n>0\):

Real sin_n(const Point& P, Parameters& pa = defaultParameters)
{
  Real h = pa("h");              //get the parameter h (user definition)
  Number n = getBasisIndex();    //get the index of function (start from 1)
  DiffOpType d = getDerivative();//get the derivative index
  switch(d)
  { case _id : return sqrt(2. / h) * sin(n *pi_ * P.x() / h);
    case _dx : return sqrt(2. / h) * (n * pi_/h)*cos(n *pi_ * P.x() / h);
    default : error("free_error","unsupported dif operator in sin_n()");
  }
  return 0.;
}

int main (int argc, char** argv)
{
  Mesh m(...);
  Domain omega=m.domain("Omega");
  Parameters ps(1.,"h");
  Number n=10;
  Space sp(_domain=omega, _basis=Function(sin_n, "sin_n", ps), _dim=n, _name="sinus basis");
  ...
}

Note that the C++ function is passed to the Space constructor using a Function object, taking the name of the function, a string and the Parameters object.

Note

For advanced usage, furtherer running data are available in user functions, depending on the FE computation context:

funcion

action

xlifepp::getN()

return current normal vector (reference)

xlifepp::getT()

return current tangent vector (reference)

xlifepp::getB()

return current bitangent vector (reference)

xlifepp::getElement()

return current element (reference)

xlifepp::getDof()

return current dof (reference)

xlifepp::getDomain()

return current domain (reference)

xlifepp::BasisIndex()

return current basis index (Number)

xlifepp::Derivative()

return current derivative index (DiffOpType)

To avoid type errors, use auto type to get references !

Interpolated spectral space#

An interpolated spectral space is defined from a set of interpolated functions, say vectors of another space (TermVectors, see TermMatrix and TermVector in details chapter). The following example shows how it works:

int main (int argc, char** argv)
{
  Mesh m(...);
  Domain omega=m.domain("Omega");
  Space V(_domain=omega,_interpolation=P1,_name="V",_optimiseNumbering); // FE space
  Unknown u(V, _name="u"); // FE unknown
  Real h=1;
  Parameters params(h,"h");
  Function sinBasis(sin_n, params);

  TermVectors sinInt(N); // interpolated functions
  for (Number n=0; n<N; n++)
  {
    sinBasis.parameter("basis index")=n;
    sinInt(n+1)=TermVector(u, omega, sinBasis, "c"+tostring(n));
  }
  Space S(_basis=sinInt, _name="V interpolated sin(n*pi*y)");

The Unknown object is described in the following section.

Advanced usage

It is possible to manipulate spectral basis by instantiate such objects:

SpectralBasisFun sbFun(omega, sinBasis, N, 1);
SpectralBasisInt sbInt(sinInt);

Thus, it is possible to evaluate basis functions at a point:

Real r;
Point P(1.,0.);
sbFun.function(n, P, r);
sbInt.function(n, P, r);

Be cautious, the type of returned argument r has to be consistent with the type of basis functions.

Unknowns and test functions#

Once the space is defined, the next step is to define unknowns (Unknown objects) and test functions (TestFunction objects) on this space. Unknowns and test functions are symbolic objects used to define linear or bilinear forms or to relate the vector representation of elements of the space.

Space sp(...);
Unknown u(sp, _name="u");           // scalar unknown
Unknown u2(sp, _name="u2", _dim=2); // vector unknown
TestFunction v(u, _name="v");       // scalar test function, dual of u
TestFunction v2(u2, _name="v2");    // vector test function, dual of u2

According to the problem, scalar or vector unknowns or test functions may be defined. The key _dim is dedicated to specify the dimension of a vector unknown (1 by default).

In the case of a problem with several unknowns, the order of the unknowns determines the organization of the sub-blocks in the FE global matrix, and can therefore have an effect on solving times. By default, they are sorted by the construction order, using the rank property of unknown:

Unknown u(V, _name="u", _dim=2);
Unknown p(V, _name="p");
TestFunction v(u, _name="v");
TestFunction q(p, _name="q");
cout << u.rank() << " " << v.rank() << " " << p.rank() << " " << q.rank();

This example gives 1 3 2 4.

It is possible to assign the rank of an unknown at the construction:

Unknown u(V, _name="u", _dim=2, _rank=2);  //rank 2
Unknown p(V, _name="p", _rank=1);  //rank 1
TestFunction v(u, _name="v", _rank=4); //rank 4
TestFunction q(p, _name="q", _rank=3); //rank 3
cout << u.rank() << " " << v.rank() << " " << p.rank() << " " << q.rank();

This example gives 2 4 1 3.

Be cautious, rank has to be unique! It is not mandatory that ranks are consecutives.

The setRanks function may be used to change the ranks of a collection of unknowns:

setRanks(u, 1, p, 2, v, 11, q, 12);

Warning

Do not confuse Unknown that represents a symbolic element (dummy) of approximation space, used to define forms, with a real element of the approximation space, represented by a TermVector. Generally, defining only one Unknown is sufficient.

Dealing with collections of unknowns/test functions#

When there is a lot of unknowns it may be more friendly to work with indexed collections. This is the purpose of the classes Unknowns and TestFunctions. As an example, suppose there are 4 domains to deal with:

Strings sn("Gamma1", "Gamma2", "Gamma3", "Gamma4");
Mesh mesh2d(Rectangle(_xmin=0, _xmax=0.5, _ymin=0, _ymax=1, _nnodes=Numbers(3, 6), _domain_name="Omega", _side_names=sn), _generator=structured);
// get the mesh domains in a Domains object
Domains doms(4);
for (Number i=1; i<=4; i++) doms(i)=mesh2d.domain(i-1);
// create one space by domain (say P1)
Spaces Vs(4);
for (Number i=1; i<=4; i++) Vs(i)=Space(_domain=doms(i), _interpolation=P1, _name="V_"+tostring(i));
// create unknowns
Unknowns us(4);
for (Number i=1; i<=4; i++) us(i)=Unknown(Vs(i), _name="u"+tostring(i));
// create TestFunctions
TestFunctions vs=dualOf(us);

Other syntaxes are available:

Unknown u1(V1, _name="u1"), u2(V2, _name="u2"), u3(V3, _name="u3");
Unknowns us1(u1, u2, u3);
Unknowns us2; us2 << u1 << u2 << u3;
Unknowns usi={u1, u2, u3};

Spaces and projections#

Sometimes, it may be useful to project an element of one FE space, say \(V\) to an other FE space, say \(W\). XLiFE++ deals with projection based on a bilinear form, say \(a(.,.)\), defined on both projection spaces. The projection \(w\in W\) of \(v \in V\) is the solution of:

\[a(w, \tilde{w}) = a(v, \tilde{w}) \ \ \forall \tilde{w} \in W\]

Let \((w_i)_{i=1,n}\) a basis of \(W\) and \((v_i)_{i=1,m}\) a basis of \(V\), the above problem is equivalent to the matrix problem:

\[\mathbb A W = \mathbb B V\]

where \(\mathbb A_{ij} = a(w_j, w_i), \mathbb B_{ij} = a(v_j, w_i), v=\sum_{i=,m} V_iv_i\) and \(w=\sum_{i=1,n} W_iw_i\).

The bilinear form should be symmetric and positive on the space \(W\) for the matrix \(\mathbb A\) to be invertible. The most simple example is the \(L^2\) projection related to the bilinear form:

\[a\left( w, \tilde{w} \right)= \int_{\Omega} w\,\tilde{w}\, d\Omega.\]

Projection operations are related to the Projector class. Let us see how to construct such Projector object :

Mesh mesh2d(Rectangle(_xmin=0, _xmax=1, _ymin=0, _ymax=1, _nnodes=6, _side_names="Gamma"), _generator=structured);
Domain omega=mesh2d.domain("Omega"), gamma = mesh2d.domain("Gamma");
Space V1(_domain=omega, _interpolation=P1, _name="V1", _notOptimizeNumbering);
Unknown u1(V1, _name="u1"); TestFunction v1(u1, _name="v1");
Space V2(_domain=omega, _interpolation=P2, _name="V2", _notOptimizeNumbering);
Unknown u2(V2, _name="u2");
Space N1(_domain=omega, _interpolation=NE1_1, _name="N1", _notOptimizeNumbering);
Unknown n1(N1, _name="n1"); // Nedelec FE
// create some L2 projectors
Projector P2toP1(V2, V1, L2Projector, "P2toP1");
Projector N1toP1(N1, 1, V1, 2, L2Projector, "N1toP1");

Note the particular construction of the \(N1\) to \(V1 \times V1\) projector. Because the space \(V1 \times V1\) is never constructed, the size of the vector unknowns (1 for N1, 2 for P1) has to be specified. Note, that such projectors may be useful when exporting solution of problems in order to vizualize it; most of vizualization softwares supporting only Lagrange representation.

The projector types available are L2Projector, H1Projector, H10Projector associated respectively to the \(L^2\), \(H^1\) and \(H^1_0\) inner product.

Warning

Be cautious, \(H^1\) and \(H^1_0\) projectors are not consistent with some FE spaces (e.g. P0)!

It is also possible to use any bilinear form as projector:

BilinearForm myblf=intg(gamma, u1*v1)+intg(omega, u1*v1);
Projector P1toP2(V1, V2, myblf, "P1toP2");

The projection may be restricted to a subdomain:

Projector P2toP1_Gamma(V2, V2, gamma, L2Projector, "P2toP1_Gamma");

To get more detailis on projectors, see Projectors.