Build a finite element method on a mesh¶
The object getfem::mesh_fem
defined in getfem/getfem_mesh_fem.h
is designed to
describe a finite element method on a whole mesh, i.e. to describe the finite
element space on which some variables will be described. This is a rather complex
object which is central in GetFEM. Basically, this structure describes the finite
element method on each element of the mesh and some additional optional
transformations. It is possible to have an arbitrary number of finite element
descriptions for a single mesh. This is particularly necessary for mixed methods,
but also to describe different data on the same mesh. One can instantiate a
getfem::mesh_fem
object as follows:
getfem::mesh_fem mf(mymesh);
where mymesh
is an already existing mesh. The structure will be linked to this
mesh and will react when modifications will be done on it.
It is possible to specify element by element the finite element method, so that
mesh_fem’s of mixed types can be treated, even if they contain elements with different
dimensions. For usual elements, the connection between two elements is done when
the two elements are compatibles (same degrees of freedom on the common face). A
numeration of the degrees of freedom is automatically done with a Cuthill Mc Kee
like algorithm. You have to keep in mind that there is absolutely no connection
between the numeration of vertices of the mesh and the numeration of the degrees
of freedom. Every getfem::mesh_fem
object has its own numeration.
There are three levels in the getfem::mesh_fem
object:
The element level: one finite element method per element. It is possible to mix the dimensions of the elements and the property to be vectorial or scalar.
The optional vectorization/tensorization (the qdim in getfem jargon, see glossary). For instance to represent a displacement or a tensor field in continuum mechanics. Scalar elements are used componentwise. Note that you can mix some intrinsic vectorial elements (RaviartThomas element for instance) which will not be vectorized and scalar elements which will be.
The optional additional linear transformation (reduction) of the degrees of freedom. It consists in defining two matrices, the reduction matrix and the extension matrix. The reduction matrix should transform the basic dofs into the reduced dofs (the number of reduced dofs should be less or equal than the number of basic dofs). The extension matrix should describe the inverse transformation. The product of the reduction matrix with the extension matrix should be the identity matrix (ensuring in particular that the two matrices are of maximal rank). This optional transformation can be used to reduce the finite element space to a certain region (tipically a boundary) or to prescribe some matching conditions between non naturally compatible fems (for instance fems with different degrees).
One has to keep in mind this construction manipulating the degrees of freedom of a
getfem::mesh_fem
object.
First level: manipulating fems on specific elements¶
To select a particular finite element method on a given element, use the method:
mf.set_finite_element(i, pf);
where i
is the index of the element and pf
is the descriptor (of type
getfem::pfem
, basically a pointer to an object which inherits from getfem::virtual_fem
) of the
finite element method. Alternative forms of this member function are:
void mesh_fem::set_finite_element(const dal::bit_vector &cvs,
getfem::pfem pf);
void mesh_fem::set_finite_element(getfem::pfem pf);
which set the finite elements for either the convexes listed in the bit_vector
cvs
, or all the convexes of the mesh. Note that the last method makes a call to
the method:
void mesh_fem::set_auto_add(pfem pf);
which defines the default finite element method which will be automatically added on new elements of the mesh (this is very useful, for instance, when a refinement of the mesh is performed).
Descriptors for finite element methods and integration methods are available thanks to the following function:
getfem::pfem pf = getfem::fem_descriptor("name of method");
where "name of method"
is to be chosen among the existing methods. A name of a
method can be retrieved thanks to the following functions:
std::string femname = getfem::name_of_fem(pf);
A non exhaustive list (see Appendix A. Finite element method list or getfem/getfem_fem.h
for
exhaustive lists) of finite element methods is given by:
"FEM_PK(n,k)"
: Classical \(P_K\) methods on simplexes of dimensionn
with degreek
polynomials."FEM_QK(n,k)"
: Classical \(Q_K\) methods on parallelepiped of dimensionn
. Tensorial product of degreek
\(P_K\) method on the segment."FEM_PK_PRISM(n,k)"
: Classical methods on prism of dimensionn
. Tensorial product of two degreek
\(P_K\) method."FEM_PRODUCT(a,b)"
: Tensorial product of the two polynomial finite element methoda
andb
."FEM_PK_DISCONTINUOUS(n,k)"
: discontinuous \(P_K\) methods on simplexes of dimensionn
with degreek
polynomials.
An alternative way to obtain a Lagrange polynomial fem suitable for a given geometric transformation is to use:
getfem::pfem getfem::classical_fem(bgeot::pgeometric_trans pg,
short_type degree);
getfem::pfem getfem::classical_discontinuous_fem(bgeot::pgeometric_trans pg,
short_type degree);
The mesh_fem can call directly these functions via:
void mesh_fem::set_classical_finite_element(const dal::bit_vector &cvs,
dim_type fem_degree);
void mesh_fem::set_classical_discontinuous_finite_element(const dal::bit_vector &cvs,
dim_type fem_degree);
void mesh_fem::set_classical_finite_element(dim_type fem_degree);
void mesh_fem::set_classical_discontinuous_finite_element(dim_type fem_degree);
Some other methods:

getfem::mesh_fem::convex_index()¶
Set of indexes (a
dal::bit_vector
) on which a finite element method is defined.

getfem::mesh_fem::linked_mesh()¶
gives a reference to the linked mesh.

getfem::mesh_fem::fem_of_element(i)¶
gives a descriptor on the finite element method defined on element of index
i
(does not take into account the qdim nor the optional reduction).

getfem::mesh_fem::clear()¶
Clears the structure, no finite element method is still defined.
Examples¶
For instance if one needs to have a description of a \(P_1\) finite element method on a triangle, the way to set it is:
mf.set_finite_element(i, getfem::fem_descriptor("FEM_PK(2, 1)"));
where i
is still the index of the triangle. It is also possible to select a
particular method directly on a set of element, passing to
mf.set_finite_element
a dal::bit_vector
instead of a single index. For instance:
mf.set_finite_element(mymesh.convex_index(),
getfem::fem_descriptor("FEM_PK(2, 1)"));
selects the method on all the elements of the mesh.
Second level: optional “vectorization/tensorization”¶
If the finite element represents an unknown which is a vector field, the method mf.set_qdim(Q)
allows set the target dimension for the definition of the
target dimension \(Q\).
If the target dimension \(Q\) is set to a value different of \(1\), the
scalar FEMs (such as \(P_k\) fems etc.) are automatically “vectorized” from
the mesh_fem object point of view, i.e. each scalar degree of freedom appears \(Q\)
times in order to represent the \(Q\) components of the vector field. If an
intrinsically vectorial element is used, the target dimension of the fem
and
the one of the mesh_fem object have to match. To sum it up,
if the fem of the \(ith\) element is intrinsically a vector FEM, then:
mf.get_qdim() == mf.fem_of_element(i)>target_dim() && mf.nb_dof_of_element(i) == mf.fem_of_element(i).nb_dof()
if the fem has a
target_dim
equal to \(1\), then:mf.nb_dof_of_element(i) == mf.get_qdim()*mf.fem_of_element(i).nb_dof()
Additionally, if the field to be represented is a tensor field instead of a vector field (for instance the stress or strain tensor field in elasticity), it is possible to specify the tensor dimensions with the methods:
mf.set_qdim(dim_type M, dim_type N)
mf.set_qdim(dim_type M, dim_type N, dim_type O, dim_type P)
mf.set_qdim(const bgeot::multi_index &mii)
respectively for a tensor field of order two, four and arbitrary (but limited to 6). For most of the operations, this is equivalent to declare a vector field of the size the product of the dimensions. However, the declared tensor dimensions are taken into account into the high level generic assembly. Remember that the components inside a tensor are stored in Fortran order.
At this level are defined the basic degrees of freedom. Some methods of the
getfem::mesh_fem
allows to obtain information on the basic dofs:

getfem::mesh_fem::nb_basic_dof_of_element(i)¶
gives the number of basic degrees of freedom on the element of index
i
.

getfem::mesh_fem::ind_basic_dof_of_element(i)¶
gives a container (an array) with all the global indexes of the basic degrees of freedom of element of index
i
.

getfem::mesh_fem::point_of_basic_dof(i, j)¶
gives a
bgeot::base_node
which represents the point associated with the basic dof of local indexj
on element of indexi
.

getfem::mesh_fem::point_of_basic_dof(j)¶
gives a
bgeot::base_node
which represents the point associated with the basic dof of global indexj
.

getfem::mesh_fem::reference_point_of_basic_dof(i, j)¶
gives a
bgeot::base_node
which represents the point associated with the basic dof of local indexj
on element of indexi
in the coordinates of the reference element.

getfem::mesh_fem::first_convex_of_basic_dof(j)¶
gives the index of the first element on which the basic degree of freedom of global index
j
is defined.

getfem::mesh_fem::nb_basic_dof()¶
gives the total number of different basic degrees of freedom.

getfem::mesh_fem::get_qdim()¶
gives the target dimension
Q
.

getfem::mesh_fem::basic_dof_on_region(i)¶
Return a
dal::bit_vector
which represents the indices of basic dof which are in the set of convexes or the set of faces of indexi
(see thegetfem::mesh
object).

getfem::mesh_fem::dof_on_region(i)¶
Return a
dal::bit_vector
which represents the indices of dof which are in the set of convexes or the set of faces of indexi
(see thegetfem::mesh
object). For a reduced mesh_fem, a dof is lying on a region if its potential corresponding shape function is nonzero on this region. The extension matrix is used to make the correspondence between basic and reduced dofs.
Third level: optional linear transformation (or reduction)¶
As described above, it is possible to provide two matrices, a reduction matrix \(R\) and an extension matrix \(E\) which will describe a linear transformation of the degrees of freedom. If \(V\) is the vector of basic degrees of freedom, then \(U=RV\) will be the vector of reduced degrees of freedom. Contrarily, given a vector \(U\) of reduced dof, \(V=EU\) will correspond to a vector of basic dof. In simple cases, \(E\) will be simply the transpose of \(R\). NOTE that every line of the extension matrix should be sparse. Otherwise, each assembled matrix will be plain !
A natural condition is that \(RE = I\) where \(I\) is the identity matrix.

getfem::mesh_fem::nb_dof()¶
gives the total number of different degrees of freedom. If the optional reduction is used, this will be the number of columns of the reduction matrix. Otherwise it will return the number of basic degrees of freedom.

getfem::mesh_fem::is_reduced()¶
return a boolean. True if the reduction is used.

getfem::mesh_fem::reduction_matrix()¶
return a const reference to the reduction matrix \(R\).

getfem::mesh_fem::extension_matrix()¶
return a const reference to the extension matrix \(E\).

getfem::mesh_fem::set_reduction_matrices(R, E)¶
Set the reduction and extension matrices to
R
andE
and validate their use.

getfem::mesh_fem::set_reduction(b)¶
Where \(b\) is a boolean. Cancel the reduction if \(b\) is false and validate it if
b
is true. Ifb
is true, the extension and reduction matrices have to be set previously.

getfem::mesh_fem::reduce_to_basic_dof(idof)¶
Set the reduction and extension matrices corresponding to keep only the basic dofs present in
idof
. The parameteridof
is either adal::bit_vector
or astd::set<size_type>
. This is equivalent to the use of agetfem::partial_mesh_fem
object.
Obtaining generic mesh_fem’s¶
It is possible to use the function:
const mesh_fem &getfem::classical_mesh_fem(const getfem::mesh &mymesh, dim_type K);
to get a classical polynomial mesh_fem of order \(K\) on the given mymesh
.
The returned mesh_fem will be destroyed automatically when its linked mesh is
destroyed. All the mesh_fem built by this function are stored in a cache, which means
that calling this function twice with the same arguments will return the same mesh_fem
object. A consequence is that you should NEVER modify this mesh_fem!
The partial_mesh_fem object¶
The getfem::partial_mesh_fem
object defined in the file
getfem_partial_mesh_fem.h
allows to reduce a getfem::mesh_fem
object to a set of dofs.
The interest is this is not a complete description of a finite element method, it
refers to the original getfem::mesh_fem
and just add reduction and extension matrices. For
instance, you can reduce a mesh_fem obtained by the function
getfem::classical_mesh_fem(mesh, K)
to obtain a finite element method on a
mesh region (which can be a boundary). The getfem::partial_mesh_fem
is in
particular used to obtain multiplier description to prescribed boundary
conditions.
The declaration of a getfem::partial_mesh_fem
object is the following:
getfem::partial_mesh_fem partial_mf(mf);
Then, one has to call the adapt method as follows:
partial_mf.adapt(kept_dof, rejected_elt = dal::bit_vector());
where kept_dof
and rejected_elt
are some dal::bit_vector
. kept_dof
is the
list of dof indices of the original mesh_fem mf
to be kept. rejected_elt
is an
optional parameter that contains a list of element indices on which the
getfem::partial_mesh_fem
states that there is no finite element method. This
is to avoid unnecessary computations during assembly procedures.