These pages are intended to describe the FEniCS interface to MGIS available in the mgis.fenics Python module. This module has been developed in the spirit of providing support to the implementation of generalized nonlinear behaviours i.e. including multiple gradients and dual flux variables. It can also be used as a simple interface to standard mechanical behaviours in small or finite strain (see the plasticity demos below). We briefly describe here the general concepts underlying the module implementation.

Documented demos

Repository: a repository containing the demos sources files is available here

The provided documented demos have been designed to progressively illustrate the use of the interface and the versatility of the approach when implementing complex generalized behaviours, both on the MFront and FEniCS sides. We recommend browsing the demos in the following order:

A brief overview of the mgis.fenics module

This module has been developed based on the concepts exposed in the Elasto-plastic analysis implemented using the MFront code generator demo published on Numerical Tours of Computational Mechanics using FEniCS. We suggest reading first this demo as an introduction to this module implementation concepts. In particular, the module relies on the notion of Quadrature function spaces to express the constitutive relation at the quadrature points level.

Compared to this initial solution, the mgis.fenics module has several advantages and new features:

The MFrontNonlinearMaterial class

This class handles the loading of a MFront behaviour through MGIS. In particular, it will contain the following important attributes:

Note about finite strain behaviours

Before loading the behaviour, it is checked if the behaviour is a finite-strain one or not. In the former case, specific finite-strain options are used when calling load. Such options specify that the stress measure will be post-processed by MGIS from Cauchy to First Piola Kirchhoff (PK1) stress and that the tangent operator will be given by \(\dfrac{\partial \text{PK1}}{\partial F}\) (DPK1_DF).

A set of helper functions enables to retrieve information about MFront object names and sizes of the form


where {1} can be any of material_property, external_state_variable, internal_state_variable, gradient, flux or tangent_block and {2} can be any of names or sizes.

The MFrontNonlinearProblem class

This class handles the definition and resolution of a nonlinear problem associated with a MFrontNonlinearMaterial. Its main attributes are:

Its main methods are:

Residual and tangent bilinear form

By default, the nonlinear residual is assumed to take the following form: Find \(u\in V\) such that:

\[ \sum_{i=1}^p \int_{\Omega}\boldsymbol{\sigma}_i(u)\cdot\delta\mathbf{g}_i(v) \,\text{dx}- L(v) = 0 \quad \forall v\in V \qquad(1)\]

where the \(\boldsymbol{\sigma}_i(u)\) are a set of fluxes as defined in the MFront behaviour using @Flux and \(\mathbf{g}_i\) are the corresponding set of gradients, declared using @Gradient. \(\delta \mathbf{g}_i(v)\) denotes the directional derivative of \(\mathbf{g}_i\) along direction \(v\), a TestFunction of the function space \(V\), and \(L(v)\) is a linear form which can be expressed using standard UFL operators. From the mechanical point of view, this residual expresses the balance between internal and external forces.

This generic form is suitable for most quasi-static mechanical problems but it does not necessarily encompass all possible situations. In particular, evolution equations such as transient heat transfer cannot be expressed in this form. However, this is not a limitation since the residual can also be redefined explicitly by the user (see Transient non-linear heat equation).

Either for the default or a user-defined one, the tangent operator associated with the residual is computed automatically using either the UFL symbolic derivation ufl.derivative for simple terms (e.g. involving the unknown field \(u\)) or using tangent operator blocks defined in the MFront behaviour for the nonlinear fluxes and internal state variables. More precisely, the variation of each flux is given by:

\[ \dfrac{\partial \boldsymbol{\sigma}_i}{\partial u} = \sum_{j\in \text{blocks}(i)} \mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j}\cdot \delta \mathbf{g}_j(u) \]

where the flux \(\boldsymbol{\sigma}_i\) is assumed to depend on the gradients \(\mathbf{g}_j\) for \(j\in \text{blocks}(i)\) (which at least contains \(i\) itself in general) and \(\mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j} = \dfrac{\partial \boldsymbol{\sigma}_i}{\partial \mathbf{g}_j}\) is the tangent operator associated with the corresponding block. Variations of internal state variables are computed in the same manner.

In the default case, the tangent bilinear form therefore reads as:

\[ a_\text{tangent}(u, v) = \sum_{i=1}^p \sum_{j\in\text{blocks}(i)}\int_{\Omega}\delta\mathbf{g}_j(u)\cdot \mathbb{T}^{\boldsymbol{\sigma}_i}_{\mathbf{g}_j} \cdot\delta\mathbf{g}_i(v) \,\text{dx} \]

Flux and Gradient objects

Two helper classes have been defined to handle flux and gradient objects:

The registration concept

Inspecting the default case (1), MFront provides access to the flux names, shapes and values when performing the constitutive update and also to the corresponding gradient. The definition of the tangent operator blocks inside the MFront behaviour also gives access to the block structure \(\text{blocks}(i)\) for each flux.

The remaining information which must be provided from the FEniCS side are the unknown field \(u\) and its discretization space \(V\), the chosen integration measure \(\,\text{dx}\) (through the quadrature_degree keyword) and, finally, the expression of each declared gradients \(\mathbf{g}_i\) in terms of the unknown field \(u\).

This step is what we call registration of each gradient which will be discussed in depth in the demos. Let us just mention that the gradients can registered using the register_gradient method of the MFrontNonlinearProblem class or via an automatic procedure if the gradient name matches predefined common gradient objects e.g. "Strain", "TemperatureGradient", etc.

The MFrontOptimisationProblem class

This class is a sibing to the MFrontNonlinearProblem class. It enables to solve nonlinear optimisation problems, especially bound-constrained problems using the default PETScTAOSolver of the form:

\[ \min_{b_l \leq u \leq b_u} f(u) \] where \(b_l\) (resp. \(b_u\)) denotes a lower (resp. upper) bound on the optimisation variable \(u\in V\).

By default, the objective function \(f(u)\) corresponds to the material total energy (stored + dissipated) computed from the get_total_energy() method. The optimisation problem requires the definition of the gradient \(F(u)\) which, by default, corresponds to (1) and its jacobian which is computed as discussed before.

Current limitations

The module has been developed using FEniCS version 2019.1.0. An important reimplementation will be planned once the dolfinx project will officially release a stable version. In particular, it will aim at fixing the following current limitations: