GetFEM++ provides a generic implementation of Nitche’s method which allows to account for Dirichlet type or contact with friction boundary conditions in a weak sense without the use of Lagrange multipliers.
The method is very attractive because it transforms a Dirichlet boundary condition into a weak term similar to a Neumann boundary condition.
However, this advantage is at the cost that the implementation of Nitche’s method is model dependent, since it requires an approximation of the corresponding Neumann term.
In order to add a boundary condition with Nitsche’s method on a variable of a model, the corresponding brick needs to have access to an approximation of the Neumann term of all partial differential terms applied to this variable.
In the following, considering a variable , we will denote by
the sum of all Neumann terms on this variable.
Note that the Neumann term will often depend on the variable
but it may also depend on other variables of the model.
This is the case for instance for mixed formulations of incompressible elasticity.
The Neumann terms depend also frequently on some parameters of the model (elasticity coefficients ...) but this is assumed to be contained in its expression.
For instance, if there is a Laplace term (), applied on the variable
, the Neumann term will be
where
is the outward unit normal on the considered boundary.
If
represents the displacements of a deformable body, the Neumann term will be
, where
is the stress tensor depending on the consitutive law.
Of course, in that case
also depends on some material parameters.
If additionally a mixed incompressibility brick is added with a variable
denoting the pressure, the Neumann term on
will depend on
in the following way:
In order to allow a generic implementation in which the brick imposing Nitsche’s method will work for every partial differential term applied to the concerned variables, each brick adding a partial differential term to a model is required to give its expression via an assembly string (weak form language).
These expressions are utilized in a special method of the model object:
expr = md.Neumann_term(variable, region)
which allows to automatically derive an expression for the sum of all Neumann terms, by scanning the expressions provided by all partial differential term bricks and performing appropriate manipulations. Of course it is required that all volumic bricks were added to the model prior to the call of this method. The derivation of the Neumann term works only for second order partial differential equations. A generic implementation for higher order pde would be more complicated.
Assume that the variable is considered and that one wants to prescribe the condition
on a part of the boundary of the considered domain.
Here
is considered equal to one in the scalar case or can be either the identity matrix in the vectorial case either a singular matrix having only 1 or 0 as eigenvalues.
This allow here to prescribe only the normal or tangent component of
.
For instance if one wants to prescribe only the normal component,
will be chosen to be equal to
where
is the outward unit normal on
.
Nitsche’s method for prescribing this Dirichlet condition consists in adding the following term to the weak formulation of the problem
where and
are two parameters of Nitsche’s method and
is the test function corresponding to
.
The parameter
can be chosen positive or negative.
corresponds to the more standard method which leads to a symmetric tangent term in standard situations,
corresponds to a non-symmetric method which has the advantage of a reduced number of terms and not requiring the second derivatives of
in the nonlinear case, and
is a kind of skew-symmetric method which ensures an inconditonal coercivity (which means independent of
) at least in standard situations.
The parameter
is a kind of penalization parameter (although the method is consistent) which is taken to be
where
is taken uniform on the mesh and
is the diameter of the element
.
Note that, in standard situations, except for
the parameter
has to be taken sufficiently small in order to ensure the convergence of Nitsche’s method.
The bricks adding a Dirichlet condition with Nitsche’s method to a model are the following:
getfem::add_Dirichlet_condition_with_Nitsche_method
(model &md, const mesh_im &mim, const std::string &varname,
const std::string &Neumannterm,
const std::string &gamma0name, size_type region,
scalar_type theta = scalar_type(1),
const std::string &dataname = std::string());
This function adds a Dirichlet condition on the variable varname and the mesh region region. This region should be a boundary. Neumannterm is the expression of the Neumann term (obtained by the Green formula) described as an expression of the weak form language. This term can be obtained with md.Neumann_term(varname, region) once all volumic bricks have been added to the model. The Dirichlet condition is prescribed with Nitsche’s method. dataname is the optional right hand side of the Dirichlet condition. It could be constant or described on a fem; scalar or vector valued, depending on the variable on which the Dirichlet condition is prescribed. gamma0name is the Nitsche’s method parameter. theta is a scalar value which can be positive or negative. theta = 1 corresponds to the standard symmetric method which is conditionnaly coercive for gamma0 small. theta = -1 corresponds to the skew-symmetric method which is inconditionnaly coercive. theta = 0 is the simplest method for which the second derivative of the Neumann term is not necessary even for nonlinear problems. Returns the brick index in the model.
getfem::add_normal_Dirichlet_condition_with_Nitsche_method
(model &md, const mesh_im &mim, const std::string &varname,
const std::string &Neumannterm,
const std::string &gamma0name, size_type region,
scalar_type theta = scalar_type(1),
const std::string &dataname = std::string());
This function adds a Dirichlet condition to the normal component of the vector (or tensor) valued variable varname and the mesh region region. This region should be a boundary. Neumannterm is the expression of the Neumann term (obtained by the Green formula) described as an expression of the weak form language. This term can be obtained with md.Neumann_term(varname, region) once all volumic bricks have been added to the model. The Dirichlet condition is prescribed with Nitsche’s method. dataname is the optional right hand side of the Dirichlet condition. It could be constant or described on a fem. gamma0name is the Nitsche’s method parameter. theta is a scalar value which can be positive or negative. theta = 1 corresponds to the standard symmetric method which is conditionnaly coercive for gamma0 small. theta = -1 corresponds to the skew-symmetric method which is inconditionnaly coercive. theta = 0 is the simplest method for which the second derivative of the Neumann term is not necessary even for nonlinear problems. Returns the brick index in the model. (This brick is not fully tested)
getfem::add_generalized_Dirichlet_condition_with_Nitsche_method
(model &md, const mesh_im &mim, const std::string &varname,
const std::string &Neumannterm,
const std::string &gamma0name, size_type region, scalar_type theta,
const std::string &dataname, const std::string &Hname);
This function adds a Dirichlet condition on the variable varname and the mesh
region region.
This version is for vector field. It prescribes a condition
where
is a matrix field. The region should be a
boundary. This region should be a boundary. Neumannterm
is the expression of the Neumann term (obtained by the Green formula)
described as an expression of the weak form language. This term can be obtained with
md.Neumann_term(varname, region) once all volumic bricks have
been added to the model. The Dirichlet
condition is prescribed with Nitsche’s method.
CAUTION : the matrix H should have all eigenvalues equal to 1 or 0.
dataname is the optional
right hand side of the Dirichlet condition. It could be constant or
described on a fem. gamma0name is the
Nitsche’s method parameter. theta is a scalar value which can be
positive or negative. theta = 1 corresponds to the standard symmetric
method which is conditionnaly coercive for gamma0 small.
theta = -1 corresponds to the skew-symmetric method which is
inconditionnaly coercive. theta = 0 is the simplest method
for which the second derivative of the Neumann term is not necessary
even for nonlinear problems. Hname is the data
corresponding to the matrix field H. It has to be a constant matrix
or described on a scalar fem. Returns the brick index in the model.
(This brick is not fully tested)
We describe here the use of Nitsch’s method to prescribe a contact with Coulomb friction condition in the small deformations framework. This corresponds to a weak integral contact condition which as some similarity with the ones which use Lagrange multipliers describe in the corresponding section, see Weak integral contact condition
In order to simplify notations, let use denote by the following map which corresponds to a couple of projections:
This application make the projection of the normal part of on
and the tangential part on the ball of center
and radius
, where
is the friction coefficient.
Using this, and considering that the sliding velocity is approximated by where the expression of
and
depend on the time integration scheme used (see Weak integral contact condition), Nitsche’s term for contact with friction reads as:
where is the contact boundary,
is the Neumann term which represents here
the stress at the contact boundary and
is the
matrix
Note that for the variant with a majority of terms vanish.
The following function adds a contact condition with or without Coulomb friction on the variable varname_u and the mesh boundary region. Neumannterm is the expression of the Neumann term (obtained by the Green formula) described as an expression of the weak form language. This term can be obtained with md.Neumann_term(varname, region) once all volumic bricks have been added to the model. The contact condition is prescribed with Nitsche’s method. The rigid obstacle should be described with the data dataname_obstacle being a signed distance to the obstacle (interpolated on a finite element method). gamma0name is the Nitsche’s method parameter. theta is a scalar value which can be positive or negative. theta = 1 corresponds to the standard symmetric method which is conditionnaly coercive for gamma0 small. theta = -1 corresponds to the skew-symmetric method which is inconditionnaly coercive. theta = 0 is the simplest method for which the second derivative of the Neumann term is not necessary. The optional parameter dataexpr_friction_coeff is the friction coefficient which could be any expression of the weak form language. Returns the brick index in the model.:
getfem::add_Nitsche_contact_with_rigid_obstacle_brick
(model &md, const mesh_im &mim, const std::string &varname_u,
const std::string &Neumannterm,
const std::string &expr_obs, const std::string &dataname_gamma0,
scalar_type theta_,
std::string dataexpr_friction_coeff,
const std::string &dataname_alpha,
const std::string &dataname_wt,
size_type region);