The model bricks

Although it is possible to use assembly procedures on their own to make your finite element programs, it is now possible to use predefined bricks to build up very quickly a certain number of models. Most of the bricks are defined in getfem/getfem_modeling.h.

A model brick is basically an object which modifies a global tangent matrix and its associated right hand side. Typical modifications are insertion of the stiffness matrix for the problem considered (linear elasticity, laplacian, ...), handling of a set of contraints, Dirichlet condition, addition of a source term to the right hand side etc. The global tangent matrix and its right hand side are stored in a model_state structure.

The model state variable

The getfem::model_state object is an object which stores the state of the system and the tangent system with eventual constraints. There are two predefined model_state types:

getfem::standard_model_state
getfem::standard_complex_model_state

The second one is for models with complex degrees of freedom like Helmholtz problem. These two predefined model_state type are built with the following predefined sparse matrices and plain vectors:

getfem::modeling_standard_sparse_matrix (gmm::col_matrix<gmm::rsvector<double> >)
getfem::modeling_standard_plain_vector (std::vector<double>)
getfem::modeling_standard_complex_sparse_matrix (gmm::col_matrix<gmm::rsvector<std::complex<double> > >)
getfem::modeling_standard_complex_plain_vector (std::vector<std::complex<double> >)

But you can define your own model state type with arbitrary types of sparse matrices and plain vectors (see the file getfem/getfem_modeling.h)

Basic properties of a brick

A brick represents a basic problem (elasticity, Helmholtz, Poisson problems ...) or a modifier of such problems (addition of a Dirichlet or Neumann condition, source term, incompressibility term ...). Each brick will participate on the global linear system to be solved (the tangent system for non linear problem). A brick is an object which derives from getfem::mdbrick_abstract <MODEL_STATE> (which itself derive from the non-template class getfem::mdbrick_abstract_common_base) with the following virtual methods to be defined:

brick.proper_update() called each time the brick should update itself. In particular, this function is expected to assign the correct values to proper_nb_dof (the nb of new dof introduced by this brick), proper_nb_constraints and proper_mixed_variables. It may also precompute certain components (like stiffness matrices for linear problems).
brick.do_compute_tangent_matrix(MS, i0, j0) the brick has to compute its own part of the tangent and constraint matrices (i0 and j0 are optional arguments representing the shifts in the matrices defined in MS).
brick.do_compute_residual(MS, i0, j0) the brick has to compute its own part of the residual of the linear system and of the constraint system (i0 and j0 are the shifts in the residual vectors defined in MS).

Of course, each specific brick may have additional methods to build the brick, define some parameters and extract the solution from the model state variable. The brick may also do some extra efforts in order to avoid unecessary recomputations (see for example the K_uptodate flag of the getfem::mdbrick_abstract_linear_pde brick).

Brick parameters

Many bricks depend on one or more parameter fields. For example, the linear elasticity brick uses the two Lamé coefficients λ and µ. These Lamé coefficients are described as a field (i.e. and mesh_fem and a vector of dof values), in a template structure getfem::mdbrick_parameter <VECTOR_TYPE>.

Some problems require a matrix or a tensor field, instead of a scalar field. For example, the brick responsible for the Dirichlet condition is used to impose h(x)u(x) = r(x) on a region of the mesh. When the mesh_fem is a vector one (Q >= 1), h(x) is a Q × Q matrix field, and r(x) is a vector field of dimension Q. That case is also handled by the getfem::mdbrick_parameter structure.

Basically, this structure contains

This structure provides these methods:

field.mf() return the mesh_fem on which the field is defined.
field.get() return the current dof data of the field.
field.fsizes() return the field size, as a vector. For a scalar field, this is an empty vector, for a n× p matrix field, this is the vector [n,p], etc.
field.fsize() return the product of the elements of the vector fsizes().
field.set([mf, ] V) change the field value. The value V can be a scalar value (constant field), a vector of length fsize() to set a constant non-scalar field, or a large vector of length fsize()*mf().nb_dof() to set a non-constant field. The mesh_fem mf is an optional paramater, hence it is possible to change the mesh_fem associated to the parameter (typically it is a polynomial mesh_fem of degree 0).
set_diagonal(V) can be used with matrix fields, to set only the diagonal elements (V length should be fsize()[0] or fsize()[0]*mf().nb_dof()).

generic elliptic brick

The generic elliptic brick is a basic brick representing a term such as
-div(k∇u) = ... where u is a scalar field and the coefficient k is a positive scalar or a symmetric positive definite order two tensor field, or a symmetric positive definite order four tensor (and u a vector field). The constructor initializes the brick for a scalar constant coefficient k:

getfem::mdbrick_generic_elliptic <MODEL_STATE> brick(mim, mf_u, k = 1.0);

where mim is a variable of type getfem::mesh_im defining the integration method used, and mf_u is a getfem::mesh_fem on the same mesh. mf_u describes the finite element method used for the unknown.

A local copy of the stiffness matrix K is stored in the brick, this obviously has a memory cost but allows not to recompute it each time when compute_tangent_matrix(...) is called.
In fact this bricks cover several situations. When k is a scalar coefficient, the brick represent a laplace operator which is componentwise if the mf_u represent a vector field (mf_u.get_qdim()>1). When k is an order two tensor coefficient, the brick represent a scalar generic elliptic operator which is componentwise if mf_u is a vector field. And finally, When k is an order four tensor coefficient, the brick represents a vectorial generic elliptic operator (for example linear elasticity with a generic Hooke tensor).

A general tensor field k can be set thanks to the two functions: brick.set_coeff_dimension(d) (with d = 0, 2 or 4) sets the tensor dimension, and brick.coeff().set(mf_data, new_k) sets the value of the tensor field (mf_data could be ommited, it is an order 0 element by default).

The following additional methods are available on this brick:

brick.coeff() gives the access to the parameter k (see section on bricks parameters).
brick.set_coeff_dimension(d) Set the tensor dimension of k. d=1 for laplacian operator, d=2 for generic scalar elliptic operator and d=4 for generic vectorial elliptic operator.
brick.get_solution(MS, V) After a solve, extract the solution of the model state variable MS and put it in the vector V.
 

Source term brick

The brick getfem::mdbrick_source_term represents either a volumic source term or a Neumann condition, i.e. a term ∫Ω f.v dx or ∫Γ f.v dx in the weak formulation , with Γ a part of∂Ω. This brick works for both real or complex terms in scalar or vectorial problems. The constructor of this brick is:

getfem::mdbrick_source_term <MODEL_STATE> brick(problem, mf_data, F, 
                                               bound=-1, num_fem=0);
getfem::mdbrick_source_term <MODEL_STATE> brick(problem, 
                                               bound=-1, num_fem=0);

where problem is the problem on which the source term will be added (a scalar elliptic brick for instance), mf_data is the finite element description for the source term f, F is a vector of type MODEL_STATE::vector_type which contains the values of the source term on each degree of freedom of mf_data, bound is an optional parameter specifying on which boundary of the main mesh fem of problem the Neumann condition is applied. If this parameter is omitted, a volumic source term will be taken into account. num_fem is an optional parameter allowing to choose a fem is the problem has several fems (for example, in a mixed problem the num_fem 0 may correspond to the mesh_fem used for the velocity, and the num_fem 1 may correspond to the mesh_fem used for the pressure).

The following additional methods are available on this brick:

brick.source_term() give the access to the parameter F defining the source term (see section on bricks parameters for how to change the value of the parameter)

Constraint brick

The constraint brick getfem::mdbrick_constraint adds constraints on the degrees of freedom of a mesh_fem. This brick is a base class for the Dirichlet condition bricks for instance. It can also be used on its own to add constraints on an unknow when this unknow is not fully determinated (for instance, when only a Neumann condition is present on the boundary of the domain). This brick deals directly with the vectorial format of the unknow, i.e. with a system representing the constraints BU=R, where B is a nc × nd matrix, U is the vector of unkown corresponding to a finite element method mf_u having nd dofs and R is a vector of size nc corresponding to the right hand side of the constraints. This correspond to add nc constraints to the system. These constraints have to be independant, which means that the matrix B has to be of maximal rank.

The brick offers three different manners to take the constraints into account. This is representied by an enum in getfem/getfem_modeling.h:

The constructor of this brick is

getfem::mdbrick_constraint <MODEL_STATE> brick(problem, num_fem);

where problem is the problem on which the constraints will be added (a generic elliptic brick for instance) and num_fem is an optional parameter allowing to choose a fem is the problem has several fems (0 is the default).

The specification of the constraints is done using the method

brick.set_constraints(B, R);

The following additional methods are available on this brick:

brick.set_constraints_rhs(R) changes only the right hand side of the constraints.
brick.set_constraints_type(c) set the method to take into account the constraints. The parameter c is either getfem::AUGMENTED_CONSTRAINTS, getfem::PENALIZED_CONSTRAINTS, or getfem::ELIMINATED_CONSTRAINTS.
brick.set_penalization_parameter(eps) set the penalization parameter for the method getfem::PENALIZED_CONSTRAINTS.

Dirichlet condition brick

The getfem::mdbrick_Dirichlet brick allows to define a Dirichlet condition on a part Γ of the boundary of the domain (i.e. set the value of the unknown on this part of the boundary u = r). This brick is derived from the constraint brick and automatically set the constraints system BU=R. As a consequence, the methods of the constraint brick are available (brick.set_constraints_type(c) and brick.set_penalization_parameter(eps)).
In order to be able to treat arbitraries finite element methods, the Dirichlet condition is considered in the weak form ∫Γ u(x)v(x) dΓ= ∫Γ r(x)v(x) dΓ for all v taken in a space of convenient multipliers. This allows to describe the data r(x) on a different fem than the unknown u(x) and allows also to have a more stable condition when the unknow u(x) is describe on a complex fem like an Xfem using a standard lagrangian fem for the multipliers.

getfem::mdbrick_Dirichlet <MODEL_STATE> brick(problem, bound, mf_mult, num_fem);

where problem is the problem on which the Dirichlet condition will be added, bound specifies on which boundary of the main mesh of problem the Dirichlet condition is applied, mf_mult is an optional parameter representing the fem for the multipliers (the default value is to take the same fem as the unknow) and num_fem is an optional parameter allowing to choose a fem is the problem has several fems (0 is the default).

The fem mf_mult has to be chosen in order to satify the Babuska-Brezzi inf-sup condition (i.e. the fact that the matrix B, which represent here a mass matrix on the boundary Γ, is of maximal rank). It is satisfied when mf_mult is the same fem as the one for the unkown and generally when mf_mult is "less rich" than this fem.

By default, the prescribed value is zero. For non-homogenous Dirichlet condition u = r, the parameter rhs of the brick has to be set by the command:

brick.rhs().set(mf_data, R);

where mf_data is the finite element description for the data and F is a vector of type MODEL_STATE::vector_type which contains the values of the data on each degree of freedom of mf_data.

The following additional methods are available on this brick:

brick.rhs() gives the access to the parameter representing the value of the Dirichlet condition.
brick.set_constraints_type(c) set the method to take into account the constraints. The parameter c is either getfem::AUGMENTED_CONSTRAINTS, getfem::PENALIZED_CONSTRAINTS, or getfem::ELIMINATED_CONSTRAINTS.
brick.set_penalization_parameter(eps) set the penalization parameter for the method getfem::PENALIZED_CONSTRAINTS.

Remark: except for the getfem::AUGMENTED_CONSTRAINTS option, an algorithm of simplification tries when it is possible to have a matrix B with only one element per line. This is possible when the mass matrix on Γ of mf_u the fem for the unknown and mf_mult is invertible.

Example of a complete Poisson problem

If mf_u and mf_data are valid finite element description on a mesh representing the domain on which the problem is defined, the following sequence will define a Poisson (laplacian) problem with a Dirichlet condition on boundary 5 of mf_u and a Neumann condition on boundary 7 of mf_u:

typedef getfem::modeling_standard_plain_vector  plain_vector;

plain_vector U(mf_u.nb_dof()), F(mf_data.nb_dof());

getfem::mdbrick_generic_elliptic <> laplacian(mim, mf_u);

plain_vector F(mf_data.nb_dof());
for (int i = 0; i < mf_data.nb_dof(); ++i) F(i) = ...;
getfem::mdbrick_source_term <> volumic_source_term(laplacian, mf_data, F);

for (int i = 0; i < mf_data.nb_dof(); ++i) F(i) = ...;
getfem::mdbrick_source_term <> neumann_condition(volumic_source_term, mf_data, F, 7);

for (int i = 0; i < mf_data.nb_dof(); ++i) F(i) = ...;
getfem::mdbrick_Dirichlet <> final_model(neumann_condition, 5);
final_model.rhs().set(mf_data, F);

getfem::standard_model_state MS(final_model);
gmm::iteration iter(residual, 1, 40000);

getfem::standard_solve(MS, final_model, iter);

laplacian.get_solution(MS, U);

Remark how the bricks are linked, each condition is applied to the brick defined with the previous condition. The order of the conditions is of course arbitrary, you can define the Dirichlet condition before the source term for instance.

Predefined solvers

Of course, for many problems, it will be more convenient to make a specific solver. Even so, one generic solver is at the moment available to test your models quickly. It can also be taken as a model to build your own solvers. It is defined in getfem/getfem_model_solvers.h and the call is

getfem::standard_solve(MS, problem, iter);

where MS is a model state variable, problem is the brick that represent your global problem and iter is an iteration object from Gmm++. See also the previous section for an example of use.
Note that SuperLu is used by default on "small" problems. You can also link MUMPS with Getfem (see section *) and used the parallele version.

Isotropic linearized elasticity brick

The getfem::mdbrick_isotropic_linearized_elasticity is a basic brick representing a term such as
-div(σ) = ...;

with
σ= λtr(ε(u))I + 2µε(u); ε(u) = (∇u + ∇u^T)/2.

ε(u) is the small strain tensor, σ is the stress tensor, λ and µ are the Lamé coefficients. This represents the system of linearized isotropic elasticity. It can also be used with λ=0 together with the linear incompressible brick to build the Stokes problem.

The constructors build the brick for constant Lamé coefficients:

getfem::mdbrick_isotropic_linearized_elasticity <MODEL_STATE>
   brick(mim, mf_u);

where mim is a variable of type getfem::mesh_im defining the integration method used, mf_u is a valid fem descriptor. mf_u describe the finite element method used for the unknown.

The brick has two parameters, lambda() and mu() for the usual Lamé coefficients. As they are getfem::mdbrick_parameter, it is possible to use either a constant value (with for example brick.lambda().set(100.), or a non constant value (for example brick.lambda().set(mf_lambda, lambdav) with lambdav of type MODEL_STATE::vector_type).

The stiffness matrix is "cached" in the brick (and available with brick.get_K()), it has a memory cost but avoids unnecessary recomputations each time compute_tangent_matrix(...) is called.
The following additional methods are available on this brick:

brick.lambda() gives access to the brick parameter lambda.
brick.mu() gives access to the brick parameter mu.
brick.get_solution(MS, V) After a solve, extract the solution of the model state variable MS and put it in the vector V.
brick.compute_Von_Mises_or_Tresca(MS, mf_vm, VM, tresca_flag) Compute the Von Mises criterion (or Tresca is tresca_flag is set to true) on the displacement field stored in MS. The stress is evaluated on the mesh_fem mf_vm and stored in the vector VM.

The program elastostatic.cc in the tests directory of Getfem++ distribution can be taken as a model of use of this brick.

Qu term brick

The getfem::mdbrick_QU_term brick can be used to add boundary conditions of Fourier-Robin type like
(∂u)/(∂n) = Qu

for scalar problems, or
σn = Qu

for linearized elasticity problems. Q is a scalar field in the scalar case or a matrix field in the vectorial case. This brick works for both real or complex terms in scalar or vectorial problems.

The constructor is the following :

getfem::mdbrick_QU_term <MODEL_STATE> brick(problem, Q_diag=0.0, bound=-1, numfem=0);

where problem is the problem on which the condition will be added. Q_diag is a real for the homogeneous and diagonal case (Q = Qdia * I). bound is the number of the boundary of main_mesh_fem() on which the condition will be applied. num_fem is an optional parameter allowing to choose a fem is the problem has several fems.

The following additional methods are available on this brick:

brick.Q() gives access to the brick parameter Q.

Helmholtz brick

This brick represent the complex or real Helmholtz problem
Δu + k^2 u = ... where k the wave number is a real or complex value. For a complex version, a complex model state variable has to be used (for example getfem::standard_complex_model_state, see helmholtz.cc in the tests directory)

The constructor is

getfem::mdbrick_Helmholtz<MODEL_STATE> brick(mim, mf_u, k);

where mim is a variable of type getfem::mesh_im defining the integration method used, mf_u describe the finite element method used for the unknown. k is the (homogeneous) wave_number. It can be changed to a non-homogeneous wave number afterwards, with brick.wave_number().set(mf_k, k) etc.

The following additional methods are available on this brick:

brick.wave_number() gives access to the brick parameter k.
brick.get_solution(MS, V) After a solve, extract the solution of the model state variable MS and put it in the vector V.

linear incompressibility (or nearly incompressibility) brick

The getfem::mdbrick_linear_incomp brick adds a linear incompressibility condition (or a nearly incompressible condition) in a problem of type
div(u) = 0, (or div(u) = εp)

This constraint is enforced with Lagrange multipliers representing the pressure, introduced in a mixed formulation.

The constructor for the incompressibility condition is:

getfem::mdbrick_linear_incomp <MODEL_STATE> brick(problem, mf_p, numfem);

where problem is the problem on which the incompressibility condition is applied, mf_p is the finite element description for the pressure (be aware that the LBB inf-sup condition has to be satisfied between the main_mesh_fem() and mf_p. num_fem is an optional parameter allowing to choose a fem is the problem has several fems.

The nearly incompressibility condition is used when it is switched on with brick.set_penalized(true). The penalization parameter ε is accessed with brick.penalization_coeff().
In nearly incompressible homogeneous linearized elasticity, one has ε= 1 / λ where λ is one of the Lamé coefficient.
For instance, the following program defines a Stokes problem with a source term and an homogeneous Dirichlet condition on boundary 0. mf_u, mf_data and mf_p have to be valid finite element description on the same mesh.

typedef getfem::modeling_standard_plain_vector  plain_vector;

plain_vector U(mf_u.nb_dof()), F(mf_data.nb_dof());

double mu = 1.0;
getfem::mdbrick_isotropic_linearized_elasticity<>
  stokes(mim, mf_u, 0.0, mu);

getfem::mdbrick_linear_incomp<> incomp(stokes, mf_p);

plain_vector F(mf_data.nb_dof());
for (int i = 0; i < mf_data.nb_dof(); ++i) F(i) = ...;
getfem::mdbrick_source_term<> volumic_source_term(incomp, mf_data, F);

gmm::clear(F);
getfem::mdbrick_Dirichlet<> final_model(volumic_source_term, 0);
final_model.rhs().set(mf_data, F);

getfem::standard_model_state MS(final_model);
gmm::iteration iter(residual, 1, 40000);
getfem::standard_solve(MS, final_model, iter);

stokes.get_solution(MS, U);

A example for a nearly incompressibility condition can be found in the program tests/elastostatic.cc.

Small displacement plasticity brick

The getfem::mdbrick_plasticity brick modelizes small-displacement quasi-static plasticity problems. It is defined in getfem/getfem_plasticity.h Plasticity happens when you stress an object too much, so that even when you remove the charge, constraints remain 'trapped' into the object. When a stress is applied to an object, if the stress is small enough, the displacements stays elastic. If that stress gets greater than a constant value called stress threshold (intrinsic to the object), then the displacements becomes plastic.

Quasi-static means that we do not take inertia into account, however the algorithm used needs some kind of a time representation. It is not possible to put the charge on the model at once, as it will render false results. Instead, we have to put the charge only a bit at a time, and calculate the deformation each time, hence the 'quasi-static' name. For instance, if you wish to put a 100N/m charge on your object, you should put it by small steps (20,40,60,80,100 is ok). We have to use that method to keep the consistency of the problem.

The constructor is

getfem::mdbrick_plasticity (mesh_im &mim_, mesh_fem &mf_u_, 
                            value_type lambdai, value_type mui,
                            value_type stress_threshold, 
                            const abstract_constraints_project &t_proj_);

The stress_threshold is the 'elasticity limit' : if constraints are below that limit, the problem is an elasticity problem. Else, it is a plasticity one. t_proj_ is an instance of a constraints projection objects, such as getfem::VM_projection.

Using that constructor, you can build the problem like with any other bricks, adding the Volumic, Neumann and Dirichlet bricks as usual (look at elastostatic.cc and plasticity.cc in tests directory for examples), and call the solver with getfem::standard_solve(MS, final_model, iter);.

Once it's done, you need to know what the results are. You can know the displacement using brick.get_solution(MS, U); where MS is defined by getfem::standard_model_state MS(final_model);, and U the displacement vector.

The Von Mises constraints can be obtained with brick.compute_Von_Mises_or_Tresca(mf_vm,VM,tresca_flag) (see the description in the linearized isotropic elasticity brick).

Contact and friction conditions brick

(to be documented, see the test program tests/dynamic_friction.cc, and the source file getfem/getfem_Coulomb_friction.h)

Linearized plate brick

(to be documented, see the test program tests/plate.cc)

Many specialized bricks are defined in getfem/getfem_linearized_plates.h:

Large strain elasticity brick

The getfem::mdbrick_nonlinear_elasticity brick represents a large strain elasticity problem. It is defined in getfem/getfem_nonlinear_elasticity.h The constructor is:

getfem::mdbrick_nonlinear_elasticity <MODEL_STATE> brick
     (Hyperelastic_Law, mim, mf_u, lawparams);

where Hyperelastic_Law is an object of type getfem::abstract_hyperelastic_law representing the considered hyperelastic law. It has to be chosen between:

getfem::SaintVenant_Kirchhoff_hyperelastic_law Hyperelastic_Law;
getfem::Ciarlet_Geymonat_hyperelastic_law Hyperelastic_Law;
getfem::Mooney_Rivlin_hyperelastic_law Hyperelastic_Law;

The Saint-Venant Kirchhoff law is a linearized law defined with the two Lamé coefficients, Ciarlet Geymonat law is defined with the two Lamé coefficients and an additional coefficient and the Mooney-Rivlin law is defined with two coefficients and is to be used with the large strain incompressibility condition.
mf_u describe the finite element method used for the unknown and mf_data the finite element method used for the parameters of the selected hyperelastic law. The parameters of the hyperelastic law are supplied in the vector lawparams (by default they are constant over the mesh, but of you can change that later with brick.params().set(mf, V)).

The program nonlinear_elastostatic.cc in tests directory is an example of use of this brick with or without an incompressibility condition.

Large strain incompressibility brick

The getfem::mdbrick_nonlinear_incomp brick adds an incompressibility condition in a large strain problem of type
det(I+∇u) = 1,

For this, Lagrange multipliers representing the pressure are introduced in a mixed formulation.

The constructor is:

getfem::mdbrick_nonlinear_incomp <MODEL_STATE> brick(problem, mf_p, numfem);

where problem is the problem on which the incompressibility condition is applied, mf_p is the finite element description for the pressure (be aware that the LBB (Ladyzhenskaja-Babuska-Brezzi) inf-sup condition has to be satisfied between the main_mesh_fem() and mf_p. num_fem is an optional parameter allowing to choose a fem is the problem has several fems.

The program nonlinear_elastostatic.cc in tests directory is an example of use of this brick.