Some Standard assembly procedures (low-level generic assembly)

Procedures defined in the file getfem/getfem_assembling.h allow the assembly of stiffness matrices, mass matrices and boundary conditions for a few amount of classical partial differential equation problems. All the procedures have vectors and matrices template parameters in order to be used with any matrix library.

CAUTION: The assembly procedures do not clean the matrix/vector at the begining of the assembly in order to keep the possibility to perform several assembly operations on the same matrix/vector. Consequently, one has to clean the matrix/vector before the first assembly operation.

Laplacian (Poisson) problem

An assembling procedure is defined to solve the problem:

\[\begin{split}-\mbox{div}(a(x)\cdot\mbox{grad}(u(x))) &= f(x)\ \mbox{ in }\Omega, \\ u(x) & = U(x)\ \mbox{ on }\Gamma_D, \\ \frac{\partial u}{\partial\eta}(x) & = F(x)\ \mbox{ on }\Gamma_N,\end{split}\]

where \(\Omega\) is an open domain of arbitrary dimension, \(\Gamma_{D}\) and \(\Gamma_{N}\) are parts of the boundary of \(\Omega\), \(u(x)\) is the unknown, \(a(x)\) is a given coefficient, \(f(x)\) is a given source term, \(U(x)\) the prescribed value of \(u(x)\) on \(\Gamma_{D}\) and \(F(x)\) is the prescribed normal derivative of \(u(x)\) on \(\Gamma_{N}\). The function to be called to assemble the stiffness matrix is:

getfem::asm_stiffness_matrix_for_laplacian(SM, mim, mfu, mfd, A);

where

• SM is a matrix of any type having the right dimension (i.e. mfu.nb_dof()),

• mim is a variable of type getfem::mesh_im defining the integration method used,

• mfu is a variable of type getfem::mesh_fem and should define the finite element method for the solution,

• mfd is a variable of type getfem::mesh_fem (possibly equal to mfu) describing the finite element method on which the coefficient \(a(x)\) is defined,

• A is the (real or complex) vector of the values of this coefficient on each degree of freedom of mfd.

Both mesh_fem should use the same mesh (i.e. &mfu.linked_mesh() == &mfd.linked_mesh()).

It is important to pay attention to the fact that the integration methods stored in mim, used to compute the elementary matrices, have to be chosen of sufficient order. The order has to be determined considering the polynomial degrees of element in mfu, in mfd and the geometric transformations for non-linear cases. For example, with linear geometric transformations, if mfu is a \(P_{K}\) FEM, and mfd is a \(P_{L}\) FEM, the integration will have to be chosen of order \(\geq 2(K-1) + L\), since the elementary integrals computed during the assembly of SM are \(\int\nabla\varphi_i\nabla\varphi_j\psi_k\) (with \(\varphi_i\) the basis functions for mfu and \(\psi_i\) the basis functions for mfd).

To assemble the source term, the function to be called is:

getfem::asm_source_term(B, mim, mfu, mfd, V);

where B is a vector of any type having the correct dimension (still mfu.nb_dof()), mim is a variable of type getfem::mesh_im defining the integration method used, mfd is a variable of type getfem::mesh_fem (possibly equal to mfu) describing the finite element method on which \(f(x)\) is defined, and V is the vector of the values of \(f(x)\) on each degree of freedom of mfd.

The function asm_source_term also has an optional argument, which is a reference to a getfem::mesh_region (or just an integer i, in which case mim.linked_mesh().region(i) will be considered). Hence for the Neumann condition on \(\Gamma_{N}\), the same function:

getfem::asm_source_term(B, mim, mfu, mfd, V, nbound);

is used again, with nbound is the index of the boundary \(\Gamma_{N}\) in the linked mesh of mim, mfu and mfd.

There is two manner (well not really, since it is also possible to use Lagrange multipliers, or to use penalization) to take into account the Dirichlet condition on \(\Gamma_{D}\), changing the linear system or explicitly reduce to the kernel of the Dirichlet condition. For the first manner, the following function is defined:

getfem::assembling_Dirichlet_condition(SM, B, mfu, nbound, R);

where nbound is the index of the boundary \(\Gamma_D\) where the Dirichlet condition is applied, R is the vector of the values of \(R(x)\) on each degree of freedom of mfu. This operation should be the last one because it transforms the stiffness matrix SM. It works only for Lagrange elements. At the end, one obtains the discrete system:

\[[SM] U = B,\]

where \(U\) is the discrete unknown.

For the second manner, one should use the more general:

getfem::asm_dirichlet_constraints(H, R, mim, mf_u, mf_mult,
                                  mf_r, r, nbound).

See the Dirichlet condition as a general linear constraint that must satisfy the solution \(u\). This function does the assembly of Dirichlet conditions of type \(\int_{\Gamma} u(x)v(x) = \int_{\Gamma}r(x)v(x)\) for all \(v\) in the space of multiplier defined by mf_mult. The fem mf_mult could be often chosen equal to mf_u except when mf_u is too “complex”.

This function just assemble these constraints into a new linear system \(H u=R\), doing some additional simplification in order to obtain a “simple” constraints matrix.

Then, one should call:

ncols = getfem::Dirichlet_nullspace(H, N, R, Ud);

which will return a vector \(U_d\) which satisfies the Dirichlet condition, and an orthogonal basis \(N\) of the kernel of \(H\). Hence, the discrete system that must be solved is:

\[(N'[SM]N) U_{int}=N'(B-[SM]U_d),\]

and the solution is $U=N U_{int}+U_d$. The output matrix \(N\) should be a \(nbdof \times nbdof\) (sparse) matrix but should be resized to ncols columns. The output vector \(U_d\) should be a \(nbdof\) vector. A big advantage of this approach is to be generic, and do not prescribed for the finite element method mf_u to be of Lagrange type. If mf_u and mf_d are different, there is implicitly a projection (with respect to the \(L^2\) norm) of the data on the finite element mf_u.

If you want to treat the more general scalar elliptic equation \(\mbox{div}(A(x)\nabla u)\), where \(A(x)\) is square matrix, you should use:

getfem::asm_stiffness_matrix_for_scalar_elliptic(M, mim, mfu,
                                                 mfdata, A);

The matrix data A should be defined on mfdata. It is expected as a vector representing a \(n \times n \times nbdof\) tensor (in Fortran order), where \(n\) is the mesh dimension of mfu, and \(nbdof\) is the number of dof of mfdata.

Linear Elasticity problem

The following function assembles the stiffness matrix for linear elasticity:

getfem::asm_stiffness_matrix_for_linear_elasticity(SM, mim, mfu,
                                                   mfd, LAMBDA, MU);

where SM is a matrix of any type having the right dimension (i.e. here mfu.nb_dof()), mim is a variable of type getfem::mesh_im defining the integration method used, mfu is a variable of type getfem::mesh_fem and should define the finite element method for the solution, mfd is a variable of type getfem::mesh_fem (possibly equal to mfu) describing the finite element method on which the Lamé coefficient are defined, LAMBDA and MU are vectors of the values of Lamé coefficients on each degree of freedom of mfd.

Caution

Linear elasticity problem is a vectorial problem, so the target dimension of mfu (see mf.set_qdim(Q)) should be the same as the dimension of the mesh.

In order to assemble source term, Neumann and Dirichlet conditions, same functions as in previous section can be used.

Stokes Problem with mixed finite element method

The assembly of the mixed term \(B = - \int p\nabla.v\) is done with:

getfem::asm_stokes_B(MATRIX &B, const mesh_im &mim,
                     const mesh_fem &mf_u, const mesh_fem &mf_p);

Assembling a mass matrix

Assembly of a mass matrix between two finite elements:

getfem::asm_mass_matrix(M, mim, mf1, mf2);

It is also possible to obtain mass matrix on a boundary with the same function:

getfem::asm_mass_matrix(M, mim, mf1, mf2, nbound);

where nbound is the region index in mim.linked_mesh(), or a mesh_region object.