getfem-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Getfem-commits] (no subject)


From: Tetsuo Koyama
Subject: [Getfem-commits] (no subject)
Date: Sun, 15 Dec 2019 09:07:04 -0500 (EST)

branch: fixmisspell
commit d36accf844c37a1194aae9d1c9204f17f58010fa
Author: Tetsuo Koyama <address@hidden>
Date:   Sun Dec 15 23:01:08 2019 +0900

    :art: imgmath :arrow_right: mathjax
---
 doc/sphinx/source/conf.py                          |  11 +-
 doc/sphinx/source/matlab/examples.rst              |   4 +-
 doc/sphinx/source/matlab/pre.rst                   |   4 +-
 doc/sphinx/source/project/femdesc.rst              |   6 +-
 doc/sphinx/source/python/examples.rst              |   4 +-
 doc/sphinx/source/python/pre.rst                   |   4 +-
 doc/sphinx/source/tutorial/thermo_coupling.rst     |  16 +--
 doc/sphinx/source/tutorial/wheel.rst               |   6 +-
 doc/sphinx/source/userdoc/appendixA.rst            |  70 ++++++------
 doc/sphinx/source/userdoc/appendixB.rst            |  24 ++---
 doc/sphinx/source/userdoc/gasm_high.rst            |   2 +-
 doc/sphinx/source/userdoc/model_ALE_rotating.rst   |  36 +++----
 doc/sphinx/source/userdoc/model_Mindlin_plate.rst  |   6 +-
 doc/sphinx/source/userdoc/model_Nitsche.rst        |   8 +-
 .../model_contact_friction_large_sliding.rst       |  20 ++--
 .../source/userdoc/model_linear_elasticity.rst     |   4 +-
 .../source/userdoc/model_nonlinear_elasticity.rst  |   2 +-
 .../userdoc/model_plasticity_small_strain.rst      | 120 ++++++++++-----------
 doc/sphinx/source/userdoc/model_solvers.rst        |   2 +-
 19 files changed, 170 insertions(+), 179 deletions(-)

diff --git a/doc/sphinx/source/conf.py b/doc/sphinx/source/conf.py
index 1cb9c9f..7452e4f 100644
--- a/doc/sphinx/source/conf.py
+++ b/doc/sphinx/source/conf.py
@@ -27,18 +27,9 @@ user_preamble = '''
 \\usepackage{mathrsfs}
 \\usepackage{amsmath}
 \\usepackage{amssymb}
-\\newcommand\\Reel{\\rm I\\hspace{-0.15em}R}
-\\newcommand\\R{\\rm I\\hspace{-0.15em}R}
-\\newcommand{\\ds}{\\displaystyle}
-\\newcommand{\\Frac}[2]{{\\ds \\frac{\\ds #1}{\\ds #2}}}
-\\usepackage[draft]{minted}\\fvset{breaklines=true}
 % end user_preamble
 '''
 
-imgmath_use_preview = True
-imgmath_dvipng_args = ['-gamma', '1.5', '-D', '110', '-bg', 'Transparent']
-imgmath_latex_preamble = user_preamble
-
 autoclass_content = "both"
 
 _stdauthor = getfem_env('authors')
@@ -49,7 +40,7 @@ _stdauthor = getfem_env('authors')
 
 # Add any Sphinx extension module names here, as strings. They can be 
extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.imgmath','sphinx.ext.autodoc',
+extensions = ['sphinx.ext.mathjax','sphinx.ext.autodoc',
               'sphinx.ext.coverage',
               'sphinx.ext.doctest']
 
diff --git a/doc/sphinx/source/matlab/examples.rst 
b/doc/sphinx/source/matlab/examples.rst
index 55de4f4..4a045c4 100644
--- a/doc/sphinx/source/matlab/examples.rst
+++ b/doc/sphinx/source/matlab/examples.rst
@@ -59,8 +59,8 @@ The first instruction builds a new |mlab_mf| object, the 
second argument specifi
 that this object will be used to interpolate scalar fields (since the unknown
 :math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM
 to every convex (each basis function is a polynomial of degree 4, remember that
-:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while
-:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a
+:math:`P^k\rm I\hspace{-0.15em}Rightarrow` polynomials of degree :math:`k`, 
while
+:math:`Q^k\rm I\hspace{-0.15em}Rightarrow` polynomials of degree :math:`2k`). 
As :math:`Q^2` is a
 polynomial FEM, you can view the expression of its basis functions on the
 reference convex::
 
diff --git a/doc/sphinx/source/matlab/pre.rst b/doc/sphinx/source/matlab/pre.rst
index 7f50df0..6474ce9 100644
--- a/doc/sphinx/source/matlab/pre.rst
+++ b/doc/sphinx/source/matlab/pre.rst
@@ -67,8 +67,8 @@ numerical integrations on these fields. Most of the finite 
elements implemented
 |gf| are scalar (however, :math:`TR_0` and edges elements are also available). 
Of
 course, these scalar FEMs can be used to approximate each component of a vector
 field. This is done by setting the :math:`Qdim` of the |mf| to the dimension of
-the vector field (i.e. :math:`Qdim=1` :math:`\Rightarrow` scalar field,
-:math:`Qdim=2` :math:`\Rightarrow` 2D vector field etc.).
+the vector field (i.e. :math:`Qdim=1` :math:`\rm I\hspace{-0.15em}Rightarrow` 
scalar field,
+:math:`Qdim=2` :math:`\rm I\hspace{-0.15em}Rightarrow` 2D vector field etc.).
 
 When solving a PDE, one often has to use more than one FEM. The most important 
one
 will be of course the one on which is defined the solution of the PDE. But most
diff --git a/doc/sphinx/source/project/femdesc.rst 
b/doc/sphinx/source/project/femdesc.rst
index 21e7ea6..6b855f8 100644
--- a/doc/sphinx/source/project/femdesc.rst
+++ b/doc/sphinx/source/project/femdesc.rst
@@ -151,7 +151,7 @@ A geometric transformation is a polynomial application:
 
 .. math::
 
-   \tau : \widehat{T} \subset \Reel^P \longrightarrow T \subset \Reel^N,
+   \tau : \widehat{T} \subset \rm I\hspace{-0.15em}R^P \longrightarrow T 
\subset \rm I\hspace{-0.15em}R^N,
 
 which maps the reference element :math:`\widehat{T}` to the real element 
:math:`T`. The
 geometric nodes are denoted:
@@ -240,12 +240,12 @@ Finite element methods description
 ----------------------------------
 
 A finite element method is defined on a reference element
-:math:`\widehat{T}\subset\Reel^P` by a set of :math:`n_d` nodes :math:`a^i` and
+:math:`\widehat{T}\subset\rm I\hspace{-0.15em}R^P` by a set of :math:`n_d` 
nodes :math:`a^i` and
 corresponding base functions
 
 .. math::
 
-   (\widehat{\varphi})^i : \widehat{T}\subset\Reel^P \longrightarrow \Reel^Q
+   (\widehat{\varphi})^i : \widehat{T}\subset\rm I\hspace{-0.15em}R^P 
\longrightarrow \rm I\hspace{-0.15em}R^Q
 
 Denoting
 
diff --git a/doc/sphinx/source/python/examples.rst 
b/doc/sphinx/source/python/examples.rst
index f9afddd..966acc3 100644
--- a/doc/sphinx/source/python/examples.rst
+++ b/doc/sphinx/source/python/examples.rst
@@ -42,8 +42,8 @@ The first instruction builds a new |py_mf| object, the second 
argument specifies
 that this object will be used to interpolate scalar fields (since the unknown
 :math:`u` is a scalar field). The second instruction assigns the :math:`Q^2` 
FEM
 to every convex (each basis function is a polynomial of degree 4, remember that
-:math:`P^k\Rightarrow` polynomials of degree :math:`k`, while
-:math:`Q^k\Rightarrow` polynomials of degree :math:`2k`). As :math:`Q^2` is a
+:math:`P^k\rm I\hspace{-0.15em}Rightarrow` polynomials of degree :math:`k`, 
while
+:math:`Q^k\rm I\hspace{-0.15em}Rightarrow` polynomials of degree :math:`2k`). 
As :math:`Q^2` is a
 polynomial FEM, you can view the expression of its basis functions on the
 reference convex:
 
diff --git a/doc/sphinx/source/python/pre.rst b/doc/sphinx/source/python/pre.rst
index d8fddf3..4bc0b5b 100644
--- a/doc/sphinx/source/python/pre.rst
+++ b/doc/sphinx/source/python/pre.rst
@@ -67,8 +67,8 @@ numerical integrations on these fields. Most of the finite 
elements implemented
 |gf| are scalar (however, :math:`TR_0` and edges elements are also available). 
Of
 course, these scalar FEMs can be used to approximate each component of a vector
 field. This is done by setting the :math:`Qdim` of the |mf| to the dimension of
-the vector field (i.e. :math:`Qdim=1` :math:`\Rightarrow` scalar field,
-:math:`Qdim=2` :math:`\Rightarrow` 2D vector field etc.).
+the vector field (i.e. :math:`Qdim=1` :math:`\rm I\hspace{-0.15em}Rightarrow` 
scalar field,
+:math:`Qdim=2` :math:`\rm I\hspace{-0.15em}Rightarrow` 2D vector field etc.).
 
 When solving a PDE, one often has to use more than one FEM. The most important 
one
 will be of course the one on which is defined the solution of the PDE. But most
diff --git a/doc/sphinx/source/tutorial/thermo_coupling.rst 
b/doc/sphinx/source/tutorial/thermo_coupling.rst
index 8267ec1..1c0d8a2 100644
--- a/doc/sphinx/source/tutorial/thermo_coupling.rst
+++ b/doc/sphinx/source/tutorial/thermo_coupling.rst
@@ -12,7 +12,7 @@ This example aims to present a simple example of a 
multiphysics problem with a n
 The problem setting
 -------------------
 
-Let :math:`\Omega \subset \R^2` be the reference configuration of a 2D plate 
(see the geometry :ref:`here <tut-fig-meshthermo>`) of thickness 
:math:`\varepsilon` submitted to external forces, electric potential and 
heating. We will denote by  :math:`\theta : \Omega \rightarrow \R` the 
temperature field (in |degreC|),  :math:`V : \Omega \rightarrow \R` the 
electric potential field and :math:`u : \Omega \rightarrow \R^2` the membrane 
displacement field.
+Let :math:`\Omega \subset \rm I\hspace{-0.15em}R^2` be the reference 
configuration of a 2D plate (see the geometry :ref:`here <tut-fig-meshthermo>`) 
of thickness :math:`\varepsilon` submitted to external forces, electric 
potential and heating. We will denote by  :math:`\theta : \Omega \rightarrow 
\rm I\hspace{-0.15em}R` the temperature field (in |degreC|),  :math:`V : \Omega 
\rightarrow \rm I\hspace{-0.15em}R` the electric potential field and :math:`u : 
\Omega \rightarrow \rm I\hspace{-0 [...]
 
 Thermal problem
 ***************
@@ -51,7 +51,7 @@ where :math:`\sigma` is still the electrical conductivity. 
Moreover, we consider
 
 .. math::
 
-  \sigma = \Frac{1}{\rho_0(1+\alpha(\theta - T_0))},
+  \sigma = \dfrac{1}{\rho_0(1+\alpha(\theta - T_0))},
 
 where :math:`T_0` is a reference temperature (air temperature here), 
:math:`\rho_0` the resistance temperature coefficient at :math:`T_0` and 
:math:`\alpha` a second resistance temperature coefficient.
 
@@ -81,9 +81,9 @@ where :math:`F` is the force density applied on the right 
lateral boundary and :
 
 .. math::
 
-  &\lambda = \Frac{E\nu}{(1+\nu)(1-2\nu)}, \\
-  &\mu = \Frac{E}{2(1+\nu)}, \\
-  &\lambda^* = \Frac{2\lambda\mu}{\lambda+2*\mu},
+  &\lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, \\
+  &\mu = \dfrac{E}{2(1+\nu)}, \\
+  &\lambda^* = \dfrac{2\lambda\mu}{\lambda+2*\mu},
 
 from :math:`E` the Young modulus and :math:`\nu` the Poisson ratio of the 
material.
 
@@ -100,9 +100,9 @@ Weak formulation of each partial differential equation is 
obtained by multiplyin
 .. math::
 
   &\mbox{Find } \theta, V, u \mbox{ with } V = 0.1, u = 0 \mbox{ on the left 
face}, V = 0 \mbox{ on the right face}, \\
-  &\ds \int_{\Omega} \varepsilon\kappa\nabla\theta\cdot\nabla\delta_{\theta} + 
2D\theta\delta_{\theta}d\Omega = \int_{\Omega} (2DT_0 + 
\varepsilon\sigma|\nabla V|^2)\delta_{\theta} d\Omega ~~~\mbox{ for all } 
\delta_{\theta}, \\
-  &\ds \int_{\Omega} \varepsilon\sigma\nabla V\cdot\nabla\delta_V = 0 d\Omega 
~~~ \mbox{ for all } \delta_V \mbox{ satisfying } \delta_V = 0 \mbox{ on the 
left and right faces}, \\
-  &\ds \int_{\Omega} \bar{\sigma}(u):\bar{\varepsilon}(\delta_u)d\Omega = 
\int_{\Gamma_N} F\cdot \delta_u d\Gamma ~~~ \mbox{ for all } \delta_{u} \mbox{ 
satisfying } \delta_u = 0 \mbox{ on the left face},
+  & \int_{\Omega} \varepsilon\kappa\nabla\theta\cdot\nabla\delta_{\theta} + 
2D\theta\delta_{\theta}d\Omega = \int_{\Omega} (2DT_0 + 
\varepsilon\sigma|\nabla V|^2)\delta_{\theta} d\Omega ~~~\mbox{ for all } 
\delta_{\theta}, \\
+  & \int_{\Omega} \varepsilon\sigma\nabla V\cdot\nabla\delta_V = 0 d\Omega ~~~ 
\mbox{ for all } \delta_V \mbox{ satisfying } \delta_V = 0 \mbox{ on the left 
and right faces}, \\
+  & \int_{\Omega} \bar{\sigma}(u):\bar{\varepsilon}(\delta_u)d\Omega = 
\int_{\Gamma_N} F\cdot \delta_u d\Gamma ~~~ \mbox{ for all } \delta_{u} \mbox{ 
satisfying } \delta_u = 0 \mbox{ on the left face},
 
 
 where :math:`\delta_{\theta}, \delta_V, \delta_u` are the test functions 
corresponding to :math:`\theta, V, u`, respectively, :math:`\Gamma_N` denotes 
the right boundary where the density of force :math:`F` is applied and 
:math:`\bar{\sigma}:\bar{\varepsilon}` is the Frobenius scalar product between 
second order tensors.
diff --git a/doc/sphinx/source/tutorial/wheel.rst 
b/doc/sphinx/source/tutorial/wheel.rst
index 09d5a72..fc5eac5 100644
--- a/doc/sphinx/source/tutorial/wheel.rst
+++ b/doc/sphinx/source/tutorial/wheel.rst
@@ -12,7 +12,7 @@ In this example of a deformable ''wheel'' enters in contact 
with a deformable fo
 The problem setting
 -------------------
 
-Let :math:`\Omega^1 \subset \R^2` be the reference configuration of a 2D wheel 
and :math:`\Omega^2 \subset \R^2` the reference configuration of a deformable 
foundation. We consider small deformation of these two bodies (linearized 
elasticity) and the contact between them. We also consider that the rim of the 
wheel is rigid and apply a vertical force on the wheel.
+Let :math:`\Omega^1 \subset \rm I\hspace{-0.15em}R^2` be the reference 
configuration of a 2D wheel and :math:`\Omega^2 \subset \rm 
I\hspace{-0.15em}R^2` the reference configuration of a deformable foundation. 
We consider small deformation of these two bodies (linearized elasticity) and 
the contact between them. We also consider that the rim of the wheel is rigid 
and apply a vertical force on the wheel.
 
 
 Building the program
@@ -163,8 +163,8 @@ Using the defined transformation, we can write an integral 
contact condition usi
 
 .. math::
 
-  & \cdots + \ds \int_{\Gamma_c} \lambda_N(X) 
(\delta_{u^1}(X)-\delta_{u^2}(\Pi(X)))\cdot n d\Gamma \\
-  & -  \ds \int_{\Gamma_c} \left(\lambda_N(X) + \left(\lambda_N(X) + 
\Frac{1}{h_T\gamma_0}((X + u^1(X))\cdot n - (\Pi(X) - u^2(\Pi(X)))\cdot 
n\right)_-\right)\delta_{\lambda_N}(X) d\Gamma = 0 ~~~~ \forall 
\delta_{\lambda_N}, \forall \delta_{u^1}, \forall \delta_{u^2},
+  & \cdots +  \int_{\Gamma_c} \lambda_N(X) 
(\delta_{u^1}(X)-\delta_{u^2}(\Pi(X)))\cdot n d\Gamma \\
+  & -   \int_{\Gamma_c} \left(\lambda_N(X) + \left(\lambda_N(X) + 
\dfrac{1}{h_T\gamma_0}((X + u^1(X))\cdot n - (\Pi(X) - u^2(\Pi(X)))\cdot 
n\right)_-\right)\delta_{\lambda_N}(X) d\Gamma = 0 ~~~~ \forall 
\delta_{\lambda_N}, \forall \delta_{u^1}, \forall \delta_{u^2},
 
 where :math:`\Gamma_c` is the slave contact boundary, :math:`\lambda_N` is the 
contact multiplier (contact pressure), :math:`h_T` is the radius of the 
element, :math:`\Pi` is the transformation, `n` is the outward normal vector to 
the master contact boundary (here :math:`n = (0,1)`), :math:`\gamma_0` is an 
augmentation parameter, :math:`(\cdot)_-:I\hspace{-0.2em}R\rightarrow 
I\hspace{-0.2em}R_+` is the negative part and :math:`\delta_{\lambda_N}, 
\delta_{u^1}, \delta_{u^2}` are the test  [...]
 
diff --git a/doc/sphinx/source/userdoc/appendixA.rst 
b/doc/sphinx/source/userdoc/appendixA.rst
index 7d40678..b0d1321 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -125,9 +125,9 @@ node is the so-called Lagrange grid. Figures 
:ref:`ud-fig-segmentpk`.
        - :math:`P_6` element, 28 d.o.f., :math:`C^0`
 
 The number of degrees of freedom for a classical :math:`P_K` Lagrange element 
of
-dimension :math:`P` and degree :math:`K` is :math:`\Frac{(P+K)!}{P!K!}`. For
-instance, in dimension 2 :math:`(P = 2)`, this value is :math:`\Frac{(K+1)
-(K+2)}{2}` and in dimension 3 :math:`(P = 3)`, it is :math:`\Frac{(K+1) (K+2)
+dimension :math:`P` and degree :math:`K` is :math:`\dfrac{(P+K)!}{P!K!}`. For
+instance, in dimension 2 :math:`(P = 2)`, this value is :math:`\dfrac{(K+1)
+(K+2)}{2}` and in dimension 3 :math:`(P = 3)`, it is :math:`\dfrac{(K+1) (K+2)
 (K+3)}{6}`.
 
   .. _ud-fig-tetrahedronpk:
@@ -169,16 +169,16 @@ Then, the coordinate of a node can be computed as
 
 .. math::
 
-   a_{i_0, i_1, ... i_P} = \sum_{n = 0}^{P} \Frac{i_n}{K}S_n, \ \ \mbox{ for } 
K \neq 0,
+   a_{i_0, i_1, ... i_P} = \sum_{n = 0}^{P} \dfrac{i_n}{K}S_n, \ \ \mbox{ for 
} K \neq 0,
 
 where :math:`S_0, S_1, ... S_N` are the vertices of the simplex (for :math:`K 
= 0`
-the particular choice :math:`a_{0, 0, ... 0} = \ds \sum_{n = 0}^{P}
-\Frac{1}{P+1}S_n` has been chosen). Then each base function, corresponding of 
each
+the particular choice :math:`a_{0, 0, ... 0} =  \sum_{n = 0}^{P}
+\dfrac{1}{P+1}S_n` has been chosen). Then each base function, corresponding of 
each
 node :math:`a_{i_0, i_1, ... i_P}` is defined by
 
 .. math::
 
-  \phi_{i_0, i_1, ... i_P} = \prod_{n = 0}^{P} \prod_{j=0}^{i_n-1} 
\left(\Frac{K \lambda_n - j}{j+1}\right).
+  \phi_{i_0, i_1, ... i_P} = \prod_{n = 0}^{P} \prod_{j=0}^{i_n-1} 
\left(\dfrac{K \lambda_n - j}{j+1}\right).
 
 where :math:`\lambda_n` are the barycentric coordinates, i.e. the polynomials 
of
 degree 1 whose value is :math:`1` on the vertex :math:`S_n` and whose value is
@@ -218,7 +218,7 @@ the classical :math:`P_K` Lagrange element.
 
      * - :math:`K`, :math:`0 \leq K \leq 255`
        - :math:`P`, :math:`~ 1 \leq P \leq 255`
-       - :math:`\Frac{(K+P)!}{K! P!}`
+       - :math:`\dfrac{(K+P)!}{K! P!}`
        - :math:`C^0`
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -240,7 +240,7 @@ the classical :math:`P_K` Lagrange element.
 
      * - :math:`K`, :math:`0 \leq K \leq 255`
        - :math:`P`, :math:`~ 1 \leq P \leq 255`
-       - :math:`\Frac{(K+P)!}{K! P!}`
+       - :math:`\dfrac{(K+P)!}{K! P!}`
        - discontinuous
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -262,7 +262,7 @@ the classical :math:`P_K` Lagrange element.
 
      * - :math:`K`, :math:`0 \leq K \leq 255`
        - :math:`P`, :math:`~ 1 \leq P \leq 255`
-       - :math:`\Frac{(K+P)!}{K! P!}`
+       - :math:`\dfrac{(K+P)!}{K! P!}`
        - discontinuous
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -284,7 +284,7 @@ functions of the tensorial product (on the reference 
element) as
 
 .. math::
 
-  \widehat{\varphi}_{ij}(x,y) = \widehat{\varphi}^1_i(x) 
\widehat{\varphi}^2_j(y), ~~ x \in \Reel^{P^1}, y \in  \Reel^{P^2},
+  \widehat{\varphi}_{ij}(x,y) = \widehat{\varphi}^1_i(x) 
\widehat{\varphi}^2_j(y), ~~ x \in \rm I\hspace{-0.15em}R^{P^1}, y \in  \rm 
I\hspace{-0.15em}R^{P^2},
 
 where :math:`\widehat{\varphi}^1_i` and :math:`\widehat{\varphi}^2_i` are 
respectively the base functions
 of the first and second element.
@@ -388,7 +388,7 @@ not to have the same degree on each dimension. An example 
is shown on figure
 
      * - :math:`2K`, :math:`0 \leq K \leq 255`
        - :math:`P`, :math:`~ 2 \leq P \leq 255`
-       - :math:`(K+1)` :math:`\times~\Frac{(K+P-1)!}{K! (P-1)!}`
+       - :math:`(K+1)` :math:`\times~\dfrac{(K+P-1)!}{K! (P-1)!}`
        - :math:`C^0`
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -410,7 +410,7 @@ not to have the same degree on each dimension. An example 
is shown on figure
 
      * - :math:`K_1+K_2`, :math:`0 \leq K_1,K_2 \leq 255`
        - :math:`P`, :math:`~ 2 \leq P \leq 255`
-       - :math:`(K_2+1)` :math:`\times~\Frac{(K_1+P-1)!}{K_1! (P-1)!}`
+       - :math:`(K_2+1)` :math:`\times~\dfrac{(K_1+P-1)!}{K_1! (P-1)!}`
        - :math:`C^0`
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -486,7 +486,7 @@ Hierarchical elements with respect to the degree
 
      * - :math:`K`, :math:`0 \leq K\leq 255`
        - :math:`P`, :math:`~ 1 \leq P \leq 255`
-       - :math:`\Frac{(K+P)!}{K! P!}`
+       - :math:`\dfrac{(K+P)!}{K! P!}`
        - :math:`C^0`
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -530,7 +530,7 @@ Hierarchical elements with respect to the degree
 
      * - :math:`K`, :math:`0 \leq K\leq 255`
        - :math:`P`, :math:`~ 2 \leq P \leq 255`
-       - :math:`(K+1)` :math:`\times~\Frac{(K+P-1)!}{K! (P-1)!}`
+       - :math:`(K+1)` :math:`\times~\dfrac{(K+P-1)!}{K! (P-1)!}`
        - :math:`C^0`
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -609,7 +609,7 @@ Hierarchical composite elements
 
      * - :math:`K`
        - :math:`P`
-       - :math:`\Frac{(SK+P)!}{(SK)! P!}`
+       - :math:`\dfrac{(SK+P)!}{(SK)! P!}`
        - variable
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -631,7 +631,7 @@ Hierarchical composite elements
 
      * - :math:`K`
        - :math:`P`
-       - :math:`\Frac{(SK+P)!}{(SK)! P!}`
+       - :math:`\dfrac{(SK+P)!}{(SK)! P!}`
        - variable
        - No :math:`(Q = 1)`
        - Yes :math:`(M = Id)`
@@ -1360,10 +1360,10 @@ The shape functions are not polynomial ones but 
rational fractions. For the firs
 .. math::
 
    \begin{array}{l}
-   \widehat{\varphi}_{0}(x,y,z) =  
\frac{1}{4}\left(1-x-y-z+\Frac{xy}{1-z}\right), \\
-   \widehat{\varphi}_{1}(x,y,z) =  
\frac{1}{4}\left(1+x-y-z-\Frac{xy}{1-z}\right), \\
-   \widehat{\varphi}_{2}(x,y,z) =  
\frac{1}{4}\left(1-x+y-z-\Frac{xy}{1-z}\right), \\
-   \widehat{\varphi}_{3}(x,y,z) =  
\frac{1}{4}\left(1+x+y-z+\Frac{xy}{1-z}\right), \\
+   \widehat{\varphi}_{0}(x,y,z) =  
\frac{1}{4}\left(1-x-y-z+\dfrac{xy}{1-z}\right), \\
+   \widehat{\varphi}_{1}(x,y,z) =  
\frac{1}{4}\left(1+x-y-z-\dfrac{xy}{1-z}\right), \\
+   \widehat{\varphi}_{2}(x,y,z) =  
\frac{1}{4}\left(1-x+y-z-\dfrac{xy}{1-z}\right), \\
+   \widehat{\varphi}_{3}(x,y,z) =  
\frac{1}{4}\left(1+x+y-z+\dfrac{xy}{1-z}\right), \\
    \widehat{\varphi}_{4}(x,y,z) =  z.\\
    \end{array}
 
@@ -1371,26 +1371,26 @@ For the second degree, setting
 
 .. math::
 
-   \xi_0 = \Frac{1-z-x}{2}, ~~~\xi_1 = \Frac{1-z-y}{2}, ~~~\xi_2 = 
\Frac{1-z+x}{2}, ~~~\xi_3 = \Frac{1-z+y}{2}, ~~~\xi_4 = z,
+   \xi_0 = \dfrac{1-z-x}{2}, ~~~\xi_1 = \dfrac{1-z-y}{2}, ~~~\xi_2 = 
\dfrac{1-z+x}{2}, ~~~\xi_3 = \dfrac{1-z+y}{2}, ~~~\xi_4 = z,
 
 the shape functions read:
 
 .. math::
 
    \begin{array}{l}
-   \widehat{\varphi}_{0}(x,y,z) = \Frac{\xi_0 
\xi_1}{(1-\xi_4)^2}((1-\xi_4-2\xi_0)(1-\xi_4-2\xi_1) -\xi_4(1-\xi_4)), \\
-   \widehat{\varphi}_{1}(x,y,z) = 
4\Frac{\xi_0\xi_1\xi_2}{(1-\xi_4)^2}(2\xi_1-(1-\xi_4)), \\
-   \widehat{\varphi}_{2}(x,y,\xi_4) = \Frac{\xi_1 
\xi_2}{(1-\xi_4)^2}((1-\xi_4-2\xi_1)(1-\xi_4-2\xi_2) -\xi_4(1-\xi_4)), \\
-   \widehat{\varphi}_{3}(x,y,z) = 
4\Frac{\xi_3\xi_0\xi_1}{(1-\xi_4)^2}(2\xi_0-(1-\xi_4)), \\
-   \widehat{\varphi}_{4}(x,y,z) = 16\Frac{\xi_0\xi_1\xi_2\xi_3}{(1-\xi_4)^2}, 
\\
-   \widehat{\varphi}_{5}(x,y,z) = 
4\Frac{\xi_1\xi_2\xi_3}{(1-\xi_4)^2}(2\xi_2-(1-\xi_4)), \\
-   \widehat{\varphi}_{6}(x,y,z) = \Frac{\xi_3 
\xi_0}{(1-\xi_4)^2}((1-\xi_4-2\xi_3)(1-\xi_4-2\xi_0) -\xi_4(1-\xi_4)), \\
-   \widehat{\varphi}_{7}(x,y,z) = 
4\Frac{\xi_2\xi_3\xi_0}{(1-\xi_4)^2}(2\xi_3-(1-\xi_4)), \\
-   \widehat{\varphi}_{8}(x,y,z) = \Frac{\xi_2 
\xi_3}{(1-\xi_4)^2}((1-\xi_4-2\xi_2)(1-\xi_4-2\xi_3) -\xi_4(1-\xi_4)), \\
-   \widehat{\varphi}_{9}(x,y,z) = 4\Frac{\xi_4}{1-\xi_4}\xi_0\xi_1, \\
-   \widehat{\varphi}_{10}(x,y,z) = 4\Frac{\xi_4}{1-\xi_4}\xi_1\xi_2,  \\
-   \widehat{\varphi}_{11}(x,y,z) = 4\Frac{\xi_4}{1-\xi_4}\xi_3\xi_0,  \\
-   \widehat{\varphi}_{12}(x,y,z) = 4\Frac{\xi_4}{1-\xi_4}\xi_2\xi_3,  \\
+   \widehat{\varphi}_{0}(x,y,z) = \dfrac{\xi_0 
\xi_1}{(1-\xi_4)^2}((1-\xi_4-2\xi_0)(1-\xi_4-2\xi_1) -\xi_4(1-\xi_4)), \\
+   \widehat{\varphi}_{1}(x,y,z) = 
4\dfrac{\xi_0\xi_1\xi_2}{(1-\xi_4)^2}(2\xi_1-(1-\xi_4)), \\
+   \widehat{\varphi}_{2}(x,y,\xi_4) = \dfrac{\xi_1 
\xi_2}{(1-\xi_4)^2}((1-\xi_4-2\xi_1)(1-\xi_4-2\xi_2) -\xi_4(1-\xi_4)), \\
+   \widehat{\varphi}_{3}(x,y,z) = 
4\dfrac{\xi_3\xi_0\xi_1}{(1-\xi_4)^2}(2\xi_0-(1-\xi_4)), \\
+   \widehat{\varphi}_{4}(x,y,z) = 16\dfrac{\xi_0\xi_1\xi_2\xi_3}{(1-\xi_4)^2}, 
\\
+   \widehat{\varphi}_{5}(x,y,z) = 
4\dfrac{\xi_1\xi_2\xi_3}{(1-\xi_4)^2}(2\xi_2-(1-\xi_4)), \\
+   \widehat{\varphi}_{6}(x,y,z) = \dfrac{\xi_3 
\xi_0}{(1-\xi_4)^2}((1-\xi_4-2\xi_3)(1-\xi_4-2\xi_0) -\xi_4(1-\xi_4)), \\
+   \widehat{\varphi}_{7}(x,y,z) = 
4\dfrac{\xi_2\xi_3\xi_0}{(1-\xi_4)^2}(2\xi_3-(1-\xi_4)), \\
+   \widehat{\varphi}_{8}(x,y,z) = \dfrac{\xi_2 
\xi_3}{(1-\xi_4)^2}((1-\xi_4-2\xi_2)(1-\xi_4-2\xi_3) -\xi_4(1-\xi_4)), \\
+   \widehat{\varphi}_{9}(x,y,z) = 4\dfrac{\xi_4}{1-\xi_4}\xi_0\xi_1, \\
+   \widehat{\varphi}_{10}(x,y,z) = 4\dfrac{\xi_4}{1-\xi_4}\xi_1\xi_2,  \\
+   \widehat{\varphi}_{11}(x,y,z) = 4\dfrac{\xi_4}{1-\xi_4}\xi_3\xi_0,  \\
+   \widehat{\varphi}_{12}(x,y,z) = 4\dfrac{\xi_4}{1-\xi_4}\xi_2\xi_3,  \\
    \widehat{\varphi}_{13}(x,y,z) = \xi_4(2\xi_4-1). \\
    \end{array}
 
diff --git a/doc/sphinx/source/userdoc/appendixB.rst 
b/doc/sphinx/source/userdoc/appendixB.rst
index 9d1bfa0..5c9997d 100644
--- a/doc/sphinx/source/userdoc/appendixB.rst
+++ b/doc/sphinx/source/userdoc/appendixB.rst
@@ -209,9 +209,9 @@ Gauss Integration methods on dimension 2
 
          7 points, order 5
 
-         :math:`a = \Frac{6+\sqrt{15}}{21}`
+         :math:`a = \dfrac{6+\sqrt{15}}{21}`
          :math:`b = 4/7 - a`
-         :math:`c = \Frac{155+\sqrt{15}}{2400}`
+         :math:`c = \dfrac{155+\sqrt{15}}{2400}`
          :math:`d = 31/240 - c`
 
      * - .. image:: images/getfemlistintmethodtriangle6.png
@@ -500,9 +500,9 @@ Gauss Integration methods on dimension 3
 
          4 points, order 2
 
-         :math:`a = \Frac{5 - \sqrt{5}}{20}`
+         :math:`a = \dfrac{5 - \sqrt{5}}{20}`
 
-         :math:`b = \Frac{5 + 3\sqrt{5}}{20}`
+         :math:`b = \dfrac{5 + 3\sqrt{5}}{20}`
 
   .. list-table:: Integration methods on dimension 3
      :widths: 40 20 10 20
@@ -601,21 +601,21 @@ Gauss Integration methods on dimension 3
 
          15 points, order 5
 
-         :math:`a = \Frac{7 + \sqrt{15}}{34}`
+         :math:`a = \dfrac{7 + \sqrt{15}}{34}`
 
-         :math:`b = \Frac{7 - \sqrt{15}}{34}`
+         :math:`b = \dfrac{7 - \sqrt{15}}{34}`
 
-         :math:`c = \Frac{13 + 3\sqrt{15}}{34}`
+         :math:`c = \dfrac{13 + 3\sqrt{15}}{34}`
 
-         :math:`d = \Frac{13 - 3\sqrt{15}}{34}`
+         :math:`d = \dfrac{13 - 3\sqrt{15}}{34}`
 
-         :math:`e = \Frac{5 - \sqrt{15}}{20}`
+         :math:`e = \dfrac{5 - \sqrt{15}}{20}`
 
-         :math:`f = \Frac{5 + \sqrt{15}}{20}`
+         :math:`f = \dfrac{5 + \sqrt{15}}{20}`
 
-         :math:`h = \Frac{2665 - 14\sqrt{15}}{226800}`
+         :math:`h = \dfrac{2665 - 14\sqrt{15}}{226800}`
 
-         :math:`i = \Frac{2665 + 14\sqrt{15}}{226800}`
+         :math:`i = \dfrac{2665 + 14\sqrt{15}}{226800}`
 
 Others methods are:
 
diff --git a/doc/sphinx/source/userdoc/gasm_high.rst 
b/doc/sphinx/source/userdoc/gasm_high.rst
index 3305c9f..8dc50df 100644
--- a/doc/sphinx/source/userdoc/gasm_high.rst
+++ b/doc/sphinx/source/userdoc/gasm_high.rst
@@ -317,7 +317,7 @@ where ``region`` is the mesh region number.
 As another example, let us describe a simple nonlinear elasticity problem. 
Assume that we consider a Saint-Venant Kirchhoff constitutive law which means 
that we consider the following elastic energy on a body of reference 
configuration :math:`\Omega`:
 
 .. math::
-  \int_{\Omega} \Frac{\lambda}{2} (\mbox{tr}(E))^2 + \mu \mbox{tr}(E^2) dx
+  \int_{\Omega} \dfrac{\lambda}{2} (\mbox{tr}(E))^2 + \mu \mbox{tr}(E^2) dx
 
 where :math:`\lambda, \mu` are the |Lame| coefficients and  :math:`E` is the 
strain tensor given by :math:`E = (\nabla u + (\nabla u)^T + (\nabla u)^T\nabla 
u)/2`.
 
diff --git a/doc/sphinx/source/userdoc/model_ALE_rotating.rst 
b/doc/sphinx/source/userdoc/model_ALE_rotating.rst
index 833ce78..4561f01 100644
--- a/doc/sphinx/source/userdoc/model_ALE_rotating.rst
+++ b/doc/sphinx/source/userdoc/model_ALE_rotating.rst
@@ -93,24 +93,24 @@ Due to the objectivity properties of standard constitutive 
laws, the expression
 
 .. math::
 
-  \Frac{\partial \varphi}{\partial t} = \Frac{\partial \bar{\varphi}}{\partial 
t} + \dot{\theta} \nabla \bar{\varphi} A \bar{X} + \dot{Z}(t),
+  \dfrac{\partial \varphi}{\partial t} = \dfrac{\partial 
\bar{\varphi}}{\partial t} + \dot{\theta} \nabla \bar{\varphi} A \bar{X} + 
\dot{Z}(t),
 
-  \Frac{\partial^2 \varphi}{\partial t^2} = \Frac{\partial^2 
\bar{\varphi}}{\partial t^2} + 2\dot{\theta} \nabla\Frac{\partial 
\bar{\varphi}}{\partial t}A \bar{X} + 
\dot{\theta}^2\mbox{div}((\nabla\bar{\varphi}A \bar{X}) \otimes (A \bar{X}) )  
+ \ddot{\theta}\nabla\bar{\varphi}A \bar{X} + \ddot{Z}(t).
+  \dfrac{\partial^2 \varphi}{\partial t^2} = \dfrac{\partial^2 
\bar{\varphi}}{\partial t^2} + 2\dot{\theta} \nabla\dfrac{\partial 
\bar{\varphi}}{\partial t}A \bar{X} + 
\dot{\theta}^2\mbox{div}((\nabla\bar{\varphi}A \bar{X}) \otimes (A \bar{X}) )  
+ \ddot{\theta}\nabla\bar{\varphi}A \bar{X} + \ddot{Z}(t).
 
 Note that the term :math:`\dot{\theta} A \bar{X} = 
\left(\hspace*{-0.5em}\begin{array}{c} \dot{\theta}\bar{X}_2 \\ 
-\dot{\theta}\bar{X}_1 \\ 0 \end{array}\hspace*{-0.5em}\right)` is the rigid 
motion velocity vector. Now, If :math:`\Theta(t,X)` is a quantity attached to 
the material points (for instance the temperature), then, with 
:math:`\bar{\Theta}(t,\bar{X}) = \Theta(t,X)` , one simply has
 
 .. math::
 
-  \Frac{\partial \Theta}{\partial t} = \Frac{\partial \bar{\Theta}}{\partial 
t} + \dot{\theta} \nabla \bar{\Theta} A \bar{X}
+  \dfrac{\partial \Theta}{\partial t} = \dfrac{\partial \bar{\Theta}}{\partial 
t} + \dot{\theta} \nabla \bar{\Theta} A \bar{X}
 
 This should not be forgotten that a correction has to be provided for each 
evolving variable for which the time derivative intervene in the considered 
model (think for instance to platic flow for plasticity). So that certain model 
bricks canot be used directly (plastic bricks for instance).
 
 |gf| bricks for structural mechanics are mainly considering the displacement 
as the amin unknown. The expression for the displacement is the following:
 
 .. math::
-  \Frac{\partial u}{\partial t} = \Frac{\partial \bar{u}}{\partial t} + 
\dot{\theta} (I_d + \nabla \bar{u}) A \bar{X} + \dot{Z}(t),
+  \dfrac{\partial u}{\partial t} = \dfrac{\partial \bar{u}}{\partial t} + 
\dot{\theta} (I_d + \nabla \bar{u}) A \bar{X} + \dot{Z}(t),
 
-  \Frac{\partial^2 u}{\partial t^2} = \Frac{\partial^2 \bar{u}}{\partial t^2} 
+ 2\dot{\theta} \nabla\Frac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t).
+  \dfrac{\partial^2 u}{\partial t^2} = \dfrac{\partial^2 \bar{u}}{\partial 
t^2} + 2\dot{\theta} \nabla\dfrac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t).
 
 
 
@@ -122,16 +122,16 @@ Weak formulation of the transient terms
 Assuming :math:`\rho^0` the density in the reference configuration having a 
rotation symmetry, the term corresponding to acceleration in the weak 
formulation reads (with :math:`v(X) = \bar{v}(\bar{X})` a test function):
 
 .. math::
-   \int_{\Omega^0} \rho^0 \Frac{\partial^2 u}{\partial t^2}\cdot vdX =
+   \int_{\Omega^0} \rho^0 \dfrac{\partial^2 u}{\partial t^2}\cdot vdX =
 
-   \int_{\bar{\Omega}^0} \rho^0 \left[\Frac{\partial^2 \bar{u}}{\partial t^2} 
+ 2\dot{\theta} \nabla\Frac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t) \right] \cdot 
\bar{v} d\bar{X}.
+   \int_{\bar{\Omega}^0} \rho^0 \left[\dfrac{\partial^2 \bar{u}}{\partial t^2} 
+ 2\dot{\theta} \nabla\dfrac{\partial \bar{u}}{\partial t}A \bar{X} +  
\dot{\theta}^2\mbox{div}(((I_d + \nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) 
 + \ddot{\theta} (I_d + \nabla\bar{u}) A \bar{X}  + \ddot{Z}(t) \right] \cdot 
\bar{v} d\bar{X}.
 
 The third term in the right hand side can be integrated by part as follows:
 
 .. math::
    \begin{array}{rcl}
-   \ds \int_{\bar{\Omega}^0} \rho^0 \dot{\theta}^2\mbox{div}(((I_d + 
\nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) \cdot \bar{v} d\bar{X} &=& - \ds 
\int_{\bar{\Omega}^0} (\dot{\theta}^2 (I_d + \nabla\bar{u})A \bar{X})) \cdot 
(\nabla (\rho^0 \bar{v}) A \bar{X}) d\bar{X} \\
-   &&\ds + \int_{\partial \bar{\Omega}^0} \rho^0 \dot{\theta}^2 (((I_d + 
\nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) \bar{N} \cdot \bar{v} 
d\bar{\Gamma}.
+    \int_{\bar{\Omega}^0} \rho^0 \dot{\theta}^2\mbox{div}(((I_d + 
\nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) \cdot \bar{v} d\bar{X} &=& -  
\int_{\bar{\Omega}^0} (\dot{\theta}^2 (I_d + \nabla\bar{u})A \bar{X})) \cdot 
(\nabla (\rho^0 \bar{v}) A \bar{X}) d\bar{X} \\
+   && + \int_{\partial \bar{\Omega}^0} \rho^0 \dot{\theta}^2 (((I_d + 
\nabla\bar{u})A \bar{X}) \otimes (A \bar{X}) ) \bar{N} \cdot \bar{v} 
d\bar{\Gamma}.
   \end{array}
 
 Since :math:`\bar{N}` the outward unit normal vector on :math:`\partial 
\bar{\Omega}^0` is orthogonal to :math:`A \bar{X}` the boundary term is zero 
and :math:`\nabla (\rho^0 \bar{v}) = \bar{v} \otimes \nabla \rho^0   + \rho^0 
\nabla \bar{v}` and since :math:`\nabla \rho^0.(A\bar{X}) = 0` because of the 
assumption on :math:`\rho^0` to have a rotation symmetry, we have
@@ -143,9 +143,9 @@ Thus globally
 
 .. math::
    \begin{array}{rcl}
-   \ds \int_{\Omega^0} \rho^0 \Frac{\partial^2 u}{\partial t^2}\cdot vdX &=&
-   \ds \int_{\bar{\Omega}^0} \rho^0 \left[\Frac{\partial^2 \bar{u}}{\partial 
t^2} + 2\dot{\theta} \nabla\Frac{\partial \bar{u}}{\partial t}A \bar{X} + 
\ddot{\theta} \nabla\bar{u} A \bar{X}   \right] \cdot \bar{v} d\bar{X}\\
-   &&\ds - \int_{\bar{\Omega}^0} \rho^0 \dot{\theta}^2(\nabla\bar{u}A \bar{X}) 
\cdot (\nabla \bar{v} A \bar{X}) d\bar{X} - \int_{\bar{\Omega}^0} \rho^0 
(\dot{\theta}^2 A^2 \bar{X} + \ddot{\theta} A\bar{X} + \ddot{Z}(t))\cdot 
\bar{v} d\bar{X}.
+    \int_{\Omega^0} \rho^0 \dfrac{\partial^2 u}{\partial t^2}\cdot vdX &=&
+    \int_{\bar{\Omega}^0} \rho^0 \left[\dfrac{\partial^2 \bar{u}}{\partial 
t^2} + 2\dot{\theta} \nabla\dfrac{\partial \bar{u}}{\partial t}A \bar{X} + 
\ddot{\theta} \nabla\bar{u} A \bar{X}   \right] \cdot \bar{v} d\bar{X}\\
+   && - \int_{\bar{\Omega}^0} \rho^0 \dot{\theta}^2(\nabla\bar{u}A \bar{X}) 
\cdot (\nabla \bar{v} A \bar{X}) d\bar{X} - \int_{\bar{\Omega}^0} \rho^0 
(\dot{\theta}^2 A^2 \bar{X} + \ddot{\theta} A\bar{X} + \ddot{Z}(t))\cdot 
\bar{v} d\bar{X}.
    \end{array}
 
 Note that two terms can deteriorate the coercivity of the problem and thus its 
well posedness and the stability of time integration schemes: the second one 
(convection term) and the fifth one. This may oblige to use additional 
stabilization techniques for large values of the angular velocity 
:math:`\dot{\theta}`.
@@ -172,7 +172,7 @@ This section present a set of bricks facilitating the use 
of an ALE formulation
 Theoretical background
 ++++++++++++++++++++++
 
-Let us consider an object whose reference configuration :math:`\Omega^0 \in 
\R^{d}` is infinite in the direction :math:`E_1`, i.e. :math:`\Omega^0 = \R 
\times \omega^0` where :math:`\omega^0 \in \R^{d-1}`. At a time :math:`t`, only 
a "windows" of this object is considered
+Let us consider an object whose reference configuration :math:`\Omega^0 \in 
\rm I\hspace{-0.15em}R^{d}` is infinite in the direction :math:`E_1`, i.e. 
:math:`\Omega^0 = \rm I\hspace{-0.15em}R \times \omega^0` where :math:`\omega^0 
\in \rm I\hspace{-0.15em}R^{d-1}`. At a time :math:`t`, only a "windows" of 
this object is considered
 
 .. math::
    \Omega^{0t} = (\alpha + z(t), \beta + z(t)) \times \omega^0
@@ -184,7 +184,7 @@ If :math:`\varphi(t, X)` is the deformation of the body 
which maps the reference
 .. math::
    \bar{\Omega}^{0} = (\alpha, \beta) \times \omega^0
 
-and :math:`\bar{\varphi}(t, X) : \R_+ \times \bar{\Omega}^{0} \rightarrow 
\R^d` defined by
+and :math:`\bar{\varphi}(t, X) : \rm I\hspace{-0.15em}R_+ \times 
\bar{\Omega}^{0} \rightarrow \rm I\hspace{-0.15em}R^d` defined by
 
 .. math::
   \bar{\varphi}(t,\bar{X}) = \varphi(t,X), ~~~\mbox{ with } \bar{X} = X - Z(t),
@@ -208,9 +208,9 @@ If we denote
 the displacement on the intermediary configuration, then it is easy to check 
that
 
 .. math::
-   \Frac{\partial \varphi}{\partial t} = \Frac{\partial \bar{u}}{\partial t} - 
\nabla \bar{u} \dot{Z}
+   \dfrac{\partial \varphi}{\partial t} = \dfrac{\partial \bar{u}}{\partial t} 
- \nabla \bar{u} \dot{Z}
 
-   \Frac{\partial^2 \varphi}{\partial t^2} = \Frac{\partial^2 
\bar{u}}{\partial t^2} - \nabla\Frac{\partial \bar{u}}{\partial t}\dot{Z} + 
\Frac{\partial^2 \bar{u}}{\partial \dot{Z}^2} - \nabla\bar{u}\ddot{Z}.
+   \dfrac{\partial^2 \varphi}{\partial t^2} = \dfrac{\partial^2 
\bar{u}}{\partial t^2} - \nabla\dfrac{\partial \bar{u}}{\partial t}\dot{Z} + 
\dfrac{\partial^2 \bar{u}}{\partial \dot{Z}^2} - \nabla\bar{u}\ddot{Z}.
 
 
 
@@ -221,9 +221,9 @@ Weak formulation of the transient terms
 Assuming :math:`\rho^0` the density in the reference being invariant with the 
considered translation, the term corresponding to acceleration in the weak 
formulation reads (with :math:`v(X) = \bar{v}(\bar{X})` a test function and 
after integration by part):
 
 .. math::
-   \int_{\Omega^0} \rho^0 \Frac{\partial^2 u}{\partial t^2}\cdot vdX =
+   \int_{\Omega^0} \rho^0 \dfrac{\partial^2 u}{\partial t^2}\cdot vdX =
 
-   \int_{\bar{\Omega}^{0}} \rho^0 \left[\Frac{\partial^2 \bar{u}}{\partial 
t^2} - 2\nabla\Frac{\partial \bar{u}}{\partial t}\dot{Z} - 
\nabla\bar{u}\ddot{Z}\right]\cdot \bar{v}  - \rho^0 
(\nabla\bar{u}\dot{Z}).(\nabla\bar{v}\dot{Z}) d\bar{X} + \int_{\partial 
\bar{\Omega}^0} \rho^0 (\nabla\bar{u}\dot{Z}).\bar{v}(\dot{Z}.\bar{N}) 
d\bar{\Gamma},
+   \int_{\bar{\Omega}^{0}} \rho^0 \left[\dfrac{\partial^2 \bar{u}}{\partial 
t^2} - 2\nabla\dfrac{\partial \bar{u}}{\partial t}\dot{Z} - 
\nabla\bar{u}\ddot{Z}\right]\cdot \bar{v}  - \rho^0 
(\nabla\bar{u}\dot{Z}).(\nabla\bar{v}\dot{Z}) d\bar{X} + \int_{\partial 
\bar{\Omega}^0} \rho^0 (\nabla\bar{u}\dot{Z}).\bar{v}(\dot{Z}.\bar{N}) 
d\bar{\Gamma},
 
 where :math:`\bar{N}` is the outward unit normal vector on :math:`\partial 
\bar{\Omega}^0`. Note that the last term vanishes on :math:`(\alpha, \beta) 
\times \partial \omega^0` but not necessarily on :math:`\{\alpha\} \times 
\omega^0` and :math:`\{\beta\} \times \omega^0`.
 
diff --git a/doc/sphinx/source/userdoc/model_Mindlin_plate.rst 
b/doc/sphinx/source/userdoc/model_Mindlin_plate.rst
index 35aac80..35f6703 100644
--- a/doc/sphinx/source/userdoc/model_Mindlin_plate.rst
+++ b/doc/sphinx/source/userdoc/model_Mindlin_plate.rst
@@ -16,7 +16,7 @@ This brick implements the classical Mindlin-Reissner bending 
model for isotropic
 The Mindlin-Reissner plate model
 ++++++++++++++++++++++++++++++++
 
-Let :math:`\Omega \subset \R^2` be the reference configuration of the 
mid-plane of a plate of thickness :math:`\epsilon`.
+Let :math:`\Omega \subset \rm I\hspace{-0.15em}R^2` be the reference 
configuration of the mid-plane of a plate of thickness :math:`\epsilon`.
 
 The weak formulation of the Mindlin-Reissner model for isotropic material can 
be written as follows for :math:`u_3` the transverse displacement and 
:math:`\theta` the rotation of fibers normal to the mid-plane:
 
@@ -25,10 +25,10 @@ The weak formulation of the Mindlin-Reissner model for 
isotropic material can be
   & \int_{\Omega} D \epsilon^3\left((1-v)\gamma(\theta):\gamma(\psi) + \nu 
\mbox{div}(\theta)\mbox{div}(\psi)\right) dx \\
   & ~~~~~~~~~~~~~~ + \int_{\Omega}G\epsilon (\nabla u_3 - \theta)\cdot(\nabla 
v_3 - \psi)dx = \int_{\Omega} F_3v_3 + M.\psi dx,
 
-for all admissible test functions :math:`v_3 : \Omega \rightarrow \R,$ $\psi : 
\Omega \rightarrow \R^2` and where:
+for all admissible test functions :math:`v_3 : \Omega \rightarrow \rm 
I\hspace{-0.15em}R,$ $\psi : \Omega \rightarrow \rm I\hspace{-0.15em}R^2` and 
where:
 
 .. math::
-  & D = \Frac{E}{12(1-\nu^2)}, ~~ G = \Frac{E\kappa}{2(1+\nu)}, \\
+  & D = \dfrac{E}{12(1-\nu^2)}, ~~ G = \dfrac{E\kappa}{2(1+\nu)}, \\
   & \gamma(\theta) = (\nabla \theta + \nabla \theta^T)/2, \\
   & F_3 = \int_{-\epsilon/2}^{\epsilon/2} f_3dx_3 + g_3^+ + g_3^-, \\
   & M_{\alpha} = \epsilon(g^+_{\alpha} - g^-_{\alpha})/2 +  
\int_{-\epsilon/2}^{\epsilon/2} x_3 f_{\alpha}dx_3, \alpha \in \{1, 2\},
diff --git a/doc/sphinx/source/userdoc/model_Nitsche.rst 
b/doc/sphinx/source/userdoc/model_Nitsche.rst
index b0bee0e..950b49f 100644
--- a/doc/sphinx/source/userdoc/model_Nitsche.rst
+++ b/doc/sphinx/source/userdoc/model_Nitsche.rst
@@ -26,7 +26,7 @@ Note that the Neumann term :math:`G` will often depend on the 
variable :math:`u`
 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 (:math:`\Delta u`), applied on the 
variable :math:`u`, the Neumann term will be :math:`G = \Frac{\partial 
u}{\partial n}` where :math:`n` is the outward unit normal on the considered 
boundary.
+For instance, if there is a Laplace term (:math:`\Delta u`), applied on the 
variable :math:`u`, the Neumann term will be :math:`G = \dfrac{\partial 
u}{\partial n}` where :math:`n` is the outward unit normal on the considered 
boundary.
 If :math:`u` represents the displacements of a deformable body, the Neumann 
term will be :math:`G = \sigma(u)n`, where :math:`\sigma(u)` is the stress 
tensor depending on the consitutive law.
 Of course, in that case :math:`G` also depends on some material parameters.
 If additionally a mixed incompressibility brick is added with a variable 
:math:`p` denoting the pressure, the Neumann term on :math:`u` will depend on 
:math:`p` in the following way:
@@ -60,7 +60,7 @@ For instance if one wants to prescribe only the normal 
component, :math:`H` will
 Nitsche's method for prescribing this Dirichlet condition consists in adding 
the following term to the weak formulation of the problem
 
 .. math::
-  \int_{\Gamma_D} \Frac{1}{\gamma}(Hu-g-\gamma HG).(Hv) - 
\theta(Hu-g).(HD_uG[v])d\Gamma,
+  \int_{\Gamma_D} \dfrac{1}{\gamma}(Hu-g-\gamma HG).(Hv) - 
\theta(Hu-g).(HD_uG[v])d\Gamma,
 
 where :math:`\gamma` and :math:`\theta` are two parameters of Nitsche's method 
and :math:`v` is the test function corresponding to :math:`u`.
 The parameter :math:`\theta` can be chosen positive or negative. :math:`\theta 
= 1` corresponds to the more standard method which leads to a symmetric tangent 
term in standard situations, :math:`\theta = 0` corresponds to a non-symmetric 
method which has the advantage of a reduced number of terms and not requiring 
the second derivatives of :math:`G` in the nonlinear case, and :math:`\theta = 
-1` is a kind of skew-symmetric method which ensures an inconditonal coercivity 
(which means inde [...]
@@ -170,13 +170,13 @@ In order to simplify notations, let use denote by 
:math:`P_{n,\mathscr{F}}` the
 .. math::
        P_{n,\mathscr{F}}(x) = -(x.n)_- n + P_{B(0,\mathscr{F}(x.n)_-)}(x - 
(x.n)n)
 
-This application make the projection of the normal part of :math:`x` on 
:math:`\Reel_-` and the tangential part on the ball of center :math:`0` and 
radius :math:`\mathscr{F}(x.n)_-`, where :math:`\mathscr{F}` is the friction 
coefficient.
+This application make the projection of the normal part of :math:`x` on 
:math:`\rm I\hspace{-0.15em}R_-` and the tangential part on the ball of center 
:math:`0` and radius :math:`\mathscr{F}(x.n)_-`, where :math:`\mathscr{F}` is 
the friction coefficient.
 
 Using this, and considering that the sliding velocity is approximated by 
:math:`\alpha(u_{_T} - w_{_T})` where the expression of :math:`\alpha` and 
:math:`w_{_T}` depend on the time integration scheme used (see 
:ref:`weak_integral_contact_section`), Nitsche's term for contact with friction 
reads as:
 
 .. math::
        &-\int_{\Gamma_C} \theta \gamma G\cdot D_u G[v] d\Gamma \\
-       &+\int_{\Gamma_C} \gamma P_{n,\mathscr{F}}(G - \Frac{Au}{\gamma} + 
\Frac{gap}{\gamma}n + \Frac{\alpha w_{_T}}{\gamma})\cdot(\theta D_u G[v] - 
\Frac{v}{\gamma}) d\Gamma.
+       &+\int_{\Gamma_C} \gamma P_{n,\mathscr{F}}(G - \dfrac{Au}{\gamma} + 
\dfrac{gap}{\gamma}n + \dfrac{\alpha w_{_T}}{\gamma})\cdot(\theta D_u G[v] - 
\dfrac{v}{\gamma}) d\Gamma.
 
 where :math:`\Gamma_C` is the contact boundary, :math:`G` is the Neumann term 
which represents here :math:`\sigma n` the stress at the contact boundary and 
:math:`A` is the :math:`d\times d` matrix
 
diff --git a/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst 
b/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
index ff67172..d7bd8c4 100644
--- a/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
+++ b/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
@@ -195,7 +195,7 @@ The list of criteria:
 Nodal contact brick with projection
 +++++++++++++++++++++++++++++++++++
 
-Notations: :math:`\Omega \subset \Reel^d` denotes the reference configuration 
of a deformable body, possibly constituted by several unconnected parts (see  
:ref:`figure<ud-fig-masterslave>`). :math:`\Omega_t` is the deformed 
configuration and :math:`\varphi^h: \Omega \rightarrow \Omega_t` is the 
approximated deformation on a finite element space :math:`V^h`. The 
displacement  :math:`u^h: \Omega \rightarrow \Reel^d` is defined by 
:math:`\varphi^h(X) = X + u^h(X)`. A generic point of the r [...]
+Notations: :math:`\Omega \subset \rm I\hspace{-0.15em}R^d` denotes the 
reference configuration of a deformable body, possibly constituted by several 
unconnected parts (see  :ref:`figure<ud-fig-masterslave>`). :math:`\Omega_t` is 
the deformed configuration and :math:`\varphi^h: \Omega \rightarrow \Omega_t` 
is the approximated deformation on a finite element space :math:`V^h`. The 
displacement  :math:`u^h: \Omega \rightarrow \rm I\hspace{-0.15em}R^d` is 
defined by :math:`\varphi^h(X) = X + [...]
 
 
 
@@ -216,8 +216,8 @@ Considering only stationnary rigid obstacles and applying 
the principle of Alart
   \left\{\begin{array}{l}
   \mbox{Find } \varphi^h \in V^h \mbox{ such that } \\
   \displaystyle \delta J(\varphi^h)[\delta u^h] - \sum_{i \in I_{\text{def}}} 
\lambda_i \cdot (\delta u^h(X_i) - \delta u^h(Y_i)) - \sum_{i \in 
I_{\text{rig}}} \lambda_i \delta u^h(X_i) = 0 ~~~ \forall \delta u^h \in V^h, \\
-  \displaystyle \Frac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i 
+ r\left(g_i n_y - \alpha(\varphi^h(X_i) - \varphi^h(Y_i) - 
W_T(X_i)+W_T(Y_i)))\right)\right]= 0  ~~\forall i \in I_{\text{def}}, \\[1em]
-  \displaystyle \Frac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i 
+ r\left(g_i n_y - \alpha(\varphi^h(X_i) - W_T(X_i)))\right)\right]= 0  
~~\forall i \in I_{\text{rig}},
+  \displaystyle \dfrac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i 
+ r\left(g_i n_y - \alpha(\varphi^h(X_i) - \varphi^h(Y_i) - 
W_T(X_i)+W_T(Y_i)))\right)\right]= 0  ~~\forall i \in I_{\text{def}}, \\[1em]
+  \displaystyle \dfrac{1}{r} \left[\lambda_i + P_{n_y, {\mathscr F}}(\lambda_i 
+ r\left(g_i n_y - \alpha(\varphi^h(X_i) - W_T(X_i)))\right)\right]= 0  
~~\forall i \in I_{\text{rig}},
   \end{array}\right.
 
 where :math:`W_T, \alpha, P_{n_y, {\mathscr F}}` ... + tangent system
@@ -239,7 +239,7 @@ The following nonlinear operators are defined in the weak 
form language (see :re
 
     .. math::
 
-      n_{trans} = \Frac{(I+ \nabla u)^{-T} n}{\|(I+\nabla u)^{-T} n\|}
+      n_{trans} = \dfrac{(I+ \nabla u)^{-T} n}{\|(I+\nabla u)^{-T} n\|}
 
     with the following partial derivatives
 
@@ -247,7 +247,7 @@ The following nonlinear operators are defined in the weak 
form language (see :re
 
       \partial_{u} n_{trans}[\delta u] = -(I - n_{trans}\otimes n_{trans})(I+ 
\nabla u)^{-T}(\nabla \delta u)^T n_{trans}
 
-      \partial_{n} n_{trans}[\delta n] = \Frac{(I+ \nabla u)^{-T}\delta n - 
n_{trans}(n_{trans}\cdot \delta n)}{\|(I+\nabla u)^{-T} n\|}
+      \partial_{n} n_{trans}[\delta n] = \dfrac{(I+ \nabla u)^{-T}\delta n - 
n_{trans}(n_{trans}\cdot \delta n)}{\|(I+\nabla u)^{-T} n\|}
 
   - ``Coulomb_friction_coupled_projection(lambda, n, Vs, g, f, r)``
     where ``lambda`` is the contact force, ``n`` is a unit normal vector, 
``Vs``
@@ -272,15 +272,15 @@ The following nonlinear operators are defined in the weak 
form language (see :re
       \left\{\begin{array}{cl}
       0 & \mbox{for } \tau \le 0 \\
       \mathbf{T}_n & \mbox{for } \|q_{_T}\| \le \tau \\
-      \Frac{\tau}{\|q_{_T}\|}
-      \left(\mathbf{T}_n - \Frac{q_{_T}}{\|q_{_T}\|}\otimes 
\Frac{q_{_T}}{\|q_{_T}\|}
+      \dfrac{\tau}{\|q_{_T}\|}
+      \left(\mathbf{T}_n - \dfrac{q_{_T}}{\|q_{_T}\|}\otimes 
\dfrac{q_{_T}}{\|q_{_T}\|}
       \right) & \mbox{otherwise }
       \end{array} \right.
 
       \partial_{\tau} P_{B(n,\tau)}(q) =
       \left\{\begin{array}{cl}
       0 & \mbox{for } \tau \le 0 \mbox{ or } \|q_{_T}\| \le \tau \\
-      \Frac{q_{_T}}{\|q_{_T}\|} & \mbox{otherwise}
+      \dfrac{q_{_T}}{\|q_{_T}\|} & \mbox{otherwise}
       \end{array} \right.
 
       \partial_n P_{B(n,\tau)}(q) =
@@ -289,9 +289,9 @@ The following nonlinear operators are defined in the weak 
form language (see :re
       0 & \mbox{for } \tau \le 0 \\
       -q \cdot n~\mathbf{T}_n - n \otimes q_{_T}
       & \mbox{for } \|q_{_T}\| \le \tau \\
-      -\Frac{\tau}{\|q_{_T}\|}
+      -\dfrac{\tau}{\|q_{_T}\|}
       \left( q \cdot n
-      \left(\mathbf{T}_n - \Frac{q_{_T}}{\|q_{_T}\|}\otimes 
\Frac{q_{_T}}{\|q_{_T}\|}
+      \left(\mathbf{T}_n - \dfrac{q_{_T}}{\|q_{_T}\|}\otimes 
\dfrac{q_{_T}}{\|q_{_T}\|}
       \right)
       + n \otimes q_{_T}
       \right) & \mbox{otherwise.}
diff --git a/doc/sphinx/source/userdoc/model_linear_elasticity.rst 
b/doc/sphinx/source/userdoc/model_linear_elasticity.rst
index 119d7b9..196956d 100644
--- a/doc/sphinx/source/userdoc/model_linear_elasticity.rst
+++ b/doc/sphinx/source/userdoc/model_linear_elasticity.rst
@@ -36,13 +36,13 @@ Let us recall that the relation between the |Lame| 
coefficients an Young modulus
 
 .. math::
 
-   \lambda = \Frac{E\nu}{(1+\nu)(1-2\nu)}, ~~~ \mu = \Frac{E}{2(1+\nu)},
+   \lambda = \dfrac{E\nu}{(1+\nu)(1-2\nu)}, ~~~ \mu = \dfrac{E}{2(1+\nu)},
 
 except for the plane stress approximation (2D model) where
 
 .. math::
 
-   \lambda^* = \Frac{E\nu}{(1-\nu^2)}, ~~~ \mu = \Frac{E}{2(1+\nu)},
+   \lambda^* = \dfrac{E\nu}{(1-\nu^2)}, ~~~ \mu = \dfrac{E}{2(1+\nu)},
 
 
 The function which adds this brick to a model and parametrized with the |Lame| 
coefficients is::
diff --git a/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst 
b/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
index 1771bc0..de47f5b 100644
--- a/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
+++ b/doc/sphinx/source/userdoc/model_nonlinear_elasticity.rst
@@ -227,7 +227,7 @@ Incompressible material.
 
 The incompressibility constraint :math:`i_3( C) = 1` is handled with a 
Lagrange multiplier :math:`p` (the pressure)
 
-constraint: :math:`\sigma = -pI \Rightarrow {\hat{\hat{\sigma}}} = 
-p\nabla\Phi\nabla\Phi^{-T}\det\nabla\Phi`
+constraint: :math:`\sigma = -pI \rm I\hspace{-0.15em}Rightarrow 
{\hat{\hat{\sigma}}} = -p\nabla\Phi\nabla\Phi^{-T}\det\nabla\Phi`
 
 .. math::
 
diff --git a/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst 
b/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
index e5308df..fcc2405 100644
--- a/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
+++ b/doc/sphinx/source/userdoc/model_plasticity_small_strain.rst
@@ -24,7 +24,7 @@ We present a short introduction to small strain plasticity. 
We refer mainly to [
 Additive decomposition of the small strain tensor
 =================================================
 
-Let :math:`\Omega \subset \R^3` be the reference configuration of a deformable 
body and :math:`u : \Omega \rightarrow \R^3` be the displacement field. Small 
strain plasticity is based on the additive decomposition of the small strain 
tensor :math:`\varepsilon(u) = \Frac{\nabla u + \nabla u^T}{2}` in
+Let :math:`\Omega \subset \rm I\hspace{-0.15em}R^3` be the reference 
configuration of a deformable body and :math:`u : \Omega \rightarrow \rm 
I\hspace{-0.15em}R^3` be the displacement field. Small strain plasticity is 
based on the additive decomposition of the small strain tensor 
:math:`\varepsilon(u) = \dfrac{\nabla u + \nabla u^T}{2}` in
 
 .. math::
    \varepsilon(u) = \varepsilon^e + \varepsilon^p
@@ -37,7 +37,7 @@ Internal variables, free energy potential and elastic law
 We consider
 
 .. math::
-   \alpha : \Omega \rightarrow \R^{d_{\alpha}},
+   \alpha : \Omega \rightarrow \rm I\hspace{-0.15em}R^{d_{\alpha}},
 
 a vector field of :math:`d_{\alpha}` strain type internal variables 
(:math:`d_{\alpha} = 0` if no internal variables are considered). We consider 
also a free energy potential
 
@@ -47,7 +47,7 @@ a vector field of :math:`d_{\alpha}` strain type internal 
variables (:math:`d_{\
 such that corresponding stress type variables are determined by
 
 .. math::
-   \sigma = \Frac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e, 
\alpha), ~~~~ A =  \Frac{\partial \psi}{\partial \alpha}(\varepsilon^e, \alpha),
+   \sigma = \dfrac{\partial \psi}{\partial \varepsilon^e}(\varepsilon^e, 
\alpha), ~~~~ A =  \dfrac{\partial \psi}{\partial \alpha}(\varepsilon^e, 
\alpha),
 
 where :math:`\sigma` is the Cauchy stress tensor and :math:`A` the stress type 
internal variables. The plastic dissipation is given by
 
@@ -74,7 +74,7 @@ The surface :math:`f(\sigma, A) = 0` is the yield surface 
where the plastic defo
 
 Let us also consider the plastic potential :math:`\Psi(\sigma, A)`, (convex 
with respect to its both variables) which determines the plastic flow direction 
in the sense that the flow rule is defined as
 
-.. math:: \dot{\varepsilon}^p = \gamma \Frac{\partial \Psi}{\partial 
\sigma}(\sigma, A), ~~~~~~ \dot{\alpha} = -\gamma \Frac{\partial \Psi}{\partial 
A}(\sigma, A),
+.. math:: \dot{\varepsilon}^p = \gamma \dfrac{\partial \Psi}{\partial 
\sigma}(\sigma, A), ~~~~~~ \dot{\alpha} = -\gamma \dfrac{\partial 
\Psi}{\partial A}(\sigma, A),
 
 with the additional complementarity condition
 
@@ -90,7 +90,7 @@ The weak formulation of a dynamic elastoplastic problem can 
be written, for an a
 .. math::
 
    \left| \begin{array}{l}
-   \ds \int_{\Omega} \rho \ddot{u}\cdot v + \sigma : \nabla v dx =  
\int_{\Omega} f\dot v dx + \int_{\Gamma_N} g\dot v dx, \\
+    \int_{\Omega} \rho \ddot{u}\cdot v + \sigma : \nabla v dx =  \int_{\Omega} 
f\dot v dx + \int_{\Gamma_N} g\dot v dx, \\
    u(0,x) = u_0(x), ~~~\dot{u}(0) = \mathrm{v}_0(x), \\
    \varepsilon^p(0,x) = \varepsilon^p_0, ~~~ \alpha(0,x) = \alpha_0,
    \end{array} \right.
@@ -112,10 +112,10 @@ Let :math:`u_{n+1}` be the displacement at the considered 
time step and  :math:`
 
 The :math:`\theta`-scheme for the integration of the plastic flow rules reads 
as
 
-.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\Delta t 
\gamma_n \Frac{\partial \Psi}{\partial \sigma}(\sigma_{n}, A_{n}) + \theta 
\Delta t \gamma_{n+1} \Frac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, 
A_{n+1}),
+.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\Delta t 
\gamma_n \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n}, A_{n}) + \theta 
\Delta t \gamma_{n+1} \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, 
A_{n+1}),
   :label: thetascheme1
 
-.. math:: \alpha_{n+1} - \alpha_n = -(1-\theta)\Delta t \gamma_n 
\Frac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) - \theta\Delta t 
\gamma_{n+1} \Frac{\partial \Psi}{\partial A}(\sigma_{n+1}, A_{n+1}),
+.. math:: \alpha_{n+1} - \alpha_n = -(1-\theta)\Delta t \gamma_n 
\dfrac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) - \theta\Delta t 
\gamma_{n+1} \dfrac{\partial \Psi}{\partial A}(\sigma_{n+1}, A_{n+1}),
   :label: thetascheme2
 
 with the complementary condition
@@ -136,9 +136,9 @@ A solution would be to solve the whole problem with all the 
unknows, that is :ma
 with :math:`\eta_n, \zeta_{n}` the right hand side of equations 
:eq:`thetascheme1`, :eq:`thetascheme2`, i.e.
 
 .. math::
-   \zeta_n = \varepsilon^p_{n} + (1-\theta)\Delta t \gamma_n \Frac{\partial 
\Psi}{\partial \sigma}(\sigma_{n}, A_{n}) ,
+   \zeta_n = \varepsilon^p_{n} + (1-\theta)\Delta t \gamma_n \dfrac{\partial 
\Psi}{\partial \sigma}(\sigma_{n}, A_{n}) ,
 
-   \eta_n = \alpha_n - (1-\theta)\Delta t \gamma_n \Frac{\partial 
\Psi}{\partial A}(\sigma_{n}, A_{n})
+   \eta_n = \alpha_n - (1-\theta)\Delta t \gamma_n \dfrac{\partial 
\Psi}{\partial A}(\sigma_{n}, A_{n})
 
 This means in particular that :math:`(\varepsilon^p_{n+1}, \alpha_{n+1}) = 
({\mathscr E}^p(u_{n+1},  \zeta_n, \eta_n), {\mathscr A}(u_{n+1}, \zeta_{n}, 
\eta_n))` is the solution to equations :eq:`thetascheme1` and 
:eq:`thetascheme2`. Both these maps and their tangent moduli (usually called 
consistent tangent moduli) are then used in the global solve of the problem 
with a Newton method and for :math:`u_{n+1}` the unique remaining variable. The 
advantage of the return mapping strategy is t [...]
 
@@ -146,14 +146,14 @@ In |gf| we propose both the return mapping strategy and 
also an alternative stra
 
 First, we consider an additional (and optional) given function 
:math:`\alpha(\sigma_{n+1}, A_{n+1}) > 0` whose interest will appear later on 
(it will allow simple local inverses) and the new unknown scalar field
 
-.. math:: \xi_{n+1} = \Frac{\gamma_{n+1}}{\alpha(\sigma_{n+1}, A_{n+1})} ,
+.. math:: \xi_{n+1} = \dfrac{\gamma_{n+1}}{\alpha(\sigma_{n+1}, A_{n+1})} ,
 
 so that our two main unknows are now :math:`u_{n+1} \mbox{ and } \xi_{n+1}`. 
The discretized plastic flow rule integration now reads:
 
-.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = 
(1-\theta)\alpha(\sigma_n,A_n)\Delta t \xi_n \Frac{\partial \Psi}{\partial 
\sigma}(\sigma_{n}, A_{n}) + \theta \alpha(\sigma_{n+1},A_{n+1}) \Delta t 
\xi_{n+1} \Frac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, A_{n+1}),
+.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = 
(1-\theta)\alpha(\sigma_n,A_n)\Delta t \xi_n \dfrac{\partial \Psi}{\partial 
\sigma}(\sigma_{n}, A_{n}) + \theta \alpha(\sigma_{n+1},A_{n+1}) \Delta t 
\xi_{n+1} \dfrac{\partial \Psi}{\partial \sigma}(\sigma_{n+1}, A_{n+1}),
   :label: flowrule1
 
-.. math:: \alpha_{n+1} - \alpha_n = (1-\theta) \alpha(\sigma_n,A_n)\Delta t 
\xi_n \Frac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) + \theta 
\alpha(\sigma_{n+1},A_{n+1}) \Delta t \xi_{n+1} \Frac{\partial \Psi}{\partial 
A}(\sigma_{n+1}, A_{n+1}),
+.. math:: \alpha_{n+1} - \alpha_n = (1-\theta) \alpha(\sigma_n,A_n)\Delta t 
\xi_n \dfrac{\partial \Psi}{\partial A}(\sigma_{n}, A_{n}) + \theta 
\alpha(\sigma_{n+1},A_{n+1}) \Delta t \xi_{n+1} \dfrac{\partial \Psi}{\partial 
A}(\sigma_{n+1}, A_{n+1}),
   :label: flowrule2
 
 .. math:: f(\sigma_{n+1}, A_{n+1}) \le 0, ~~~ \xi_{n+1} \ge 0, ~~~ 
f(\sigma_{n+1}, A_{n+1}) \xi_{n+1} = 0.
@@ -172,17 +172,17 @@ where the pair :math:`(\varepsilon^p_{n+1}, \alpha_{n+1}) 
= (\tilde{\mathscr E}^
 
 Still :math:`u_{n+1} \mbox{ and } \xi_{n+1}` be given the stress 
:math:`\sigma_{n+1}` reads
 
-.. math:: \sigma_{n+1} = \Frac{\partial \psi^e}{\partial 
\varepsilon^e}(\varepsilon(u_{n+1}) -\varepsilon^p_{n+1}).
+.. math:: \sigma_{n+1} = \dfrac{\partial \psi^e}{\partial 
\varepsilon^e}(\varepsilon(u_{n+1}) -\varepsilon^p_{n+1}).
 
-.. math:: A_{n+1} = \Frac{\partial \psi^p}{\partial \alpha}(\alpha_{n+1}).
+.. math:: A_{n+1} = \dfrac{\partial \psi^p}{\partial \alpha}(\alpha_{n+1}).
 
 The complementarity equation :eq:`flowrule3` is then prescribed with the use 
of a well chosen complementarity function, as in [HA-WO2009]_ for :math:`r > 0` 
such as:
 
-.. math:: \ds \int_{\Omega} (\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, 
A_{n+1}))_+) \lambda dx = 0,   \forall \lambda
+.. math::  \int_{\Omega} (\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, 
A_{n+1}))_+) \lambda dx = 0,   \forall \lambda
 
 or
 
-.. math:: \ds \int_{\Omega} (f(\sigma_{n+1} + (-f(\sigma_{n+1}, A_{n+1}) - 
\xi_{n+1}/r)_+ , A_{n+1}) ) \lambda dx = 0,   \forall \lambda
+.. math::  \int_{\Omega} (f(\sigma_{n+1} + (-f(\sigma_{n+1}, A_{n+1}) - 
\xi_{n+1}/r)_+ , A_{n+1}) ) \lambda dx = 0,   \forall \lambda
 
 NOTE : The notation :math:`\Delta \xi_{n+1} = \Delta t \xi_{n+1}` is often 
used in the litterature. The choice here is to preserve the distinction between 
the two quantities, mainly because ot the possible use of adaptative time step 
: when the time step is changing, the value :math:`\xi_n` has to be multiplied 
by the new time step, so that it is preferable to store :math:`\xi_n` instead 
of :math:`\Delta \xi_{n}` when using the :math:`\theta`-scheme.
 
@@ -221,7 +221,7 @@ The nonzero :math:`\sigma_{3,3}` component of the stress 
tensor is given by
 
 Note that in the common case where isochoric plastic strain is assumed, one has
 
-.. math:: \mbox{ tr}(\varepsilon^p) = 0 ~~~~ \Rightarrow  ~~~ 
\varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).
+.. math:: \mbox{ tr}(\varepsilon^p) = 0 ~~~~ \rm I\hspace{-0.15em}Rightarrow  
~~~ \varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).
 
 
 
@@ -244,21 +244,21 @@ For an isotropic linearized elastic response, one has 
:math:`\sigma = \lambda \m
 
 with
 
-.. math:: \varepsilon^e_{3,3} = 
-\Frac{\lambda}{\lambda+2\mu}(\varepsilon^e_{1,1} + \varepsilon^e_{2,2})
+.. math:: \varepsilon^e_{3,3} = 
-\dfrac{\lambda}{\lambda+2\mu}(\varepsilon^e_{1,1} + \varepsilon^e_{2,2})
 
 so that
 
-.. math:: \bar{\sigma} = \lambda^* \mbox{tr}(\bar{\varepsilon}^e) + 
2\mu\bar{\varepsilon}^e ~~~\mbox{ with } \lambda^* = 
\Frac{2\mu\lambda}{\lambda+2\mu}
+.. math:: \bar{\sigma} = \lambda^* \mbox{tr}(\bar{\varepsilon}^e) + 
2\mu\bar{\varepsilon}^e ~~~\mbox{ with } \lambda^* = 
\dfrac{2\mu\lambda}{\lambda+2\mu}
    :label: plane_stress_iso
 
 Moreover
 
-.. math:: \|\mbox{Dev}(\sigma)\| = \left(\|\bar{\sigma}\|^2 - 
\Frac{1}{3}(\mbox{tr}(\bar{\sigma}))^2\right)^{1/2}.
+.. math:: \|\mbox{Dev}(\sigma)\| = \left(\|\bar{\sigma}\|^2 - 
\dfrac{1}{3}(\mbox{tr}(\bar{\sigma}))^2\right)^{1/2}.
    :label: plane_stress_dev
 
 Note that in the case where isochoric plastic strain is assumed, one still has
 
-.. math:: \mbox{ tr}(\varepsilon^p) = 0 ~~~~ \Rightarrow  ~~~ 
\varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).
+.. math:: \mbox{ tr}(\varepsilon^p) = 0 ~~~~ \rm I\hspace{-0.15em}Rightarrow  
~~~ \varepsilon^p_{3,3} = - (\varepsilon^p_{1,1} + \varepsilon^p_{2,2}).
 
 
 Some classical laws
@@ -276,22 +276,22 @@ Perfect isotropic associated elastoplasticity with 
Von-Mises criterion (Prandl-R
 
 There is no internal variables and we consider an isotropic elastic response. 
The flow rule reads
 
-.. math:: \dot{\varepsilon}^p = \gamma 
\Frac{\mbox{Dev}(\sigma)}{\|\mbox{Dev}(\sigma)\|}
+.. math:: \dot{\varepsilon}^p = \gamma 
\dfrac{\mbox{Dev}(\sigma)}{\|\mbox{Dev}(\sigma)\|}
 
 This corresponds to :math:`\Psi(\sigma) = f(\sigma) = \|\mbox{Dev}(\sigma)\| - 
\sqrt{\frac{2}{3}}\sigma_y`.
 
 
 The :math:`\theta`-scheme for the integration of the plastic flow rule reads:
 
-.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = 
(1-\theta)\alpha(\sigma_{n}) \Delta t \xi_n 
\Frac{\mbox{Dev}(\sigma_{n})}{\|\mbox{Dev}(\sigma_{n})\|} + 
\theta\alpha(\sigma_{n+1}) \Delta t \xi_{n+1} 
\Frac{\mbox{Dev}(\sigma_{n+1})}{\|\mbox{Dev}(\sigma_{n+1})\|}.
+.. math:: \varepsilon^p_{n+1} - \varepsilon^p_{n} = 
(1-\theta)\alpha(\sigma_{n}) \Delta t \xi_n 
\dfrac{\mbox{Dev}(\sigma_{n})}{\|\mbox{Dev}(\sigma_{n})\|} + 
\theta\alpha(\sigma_{n+1}) \Delta t \xi_{n+1} 
\dfrac{\mbox{Dev}(\sigma_{n+1})}{\|\mbox{Dev}(\sigma_{n+1})\|}.
 
-Choosing the factor :math:`\alpha(\sigma_{n}) = \|\mbox{Dev}(\sigma_{n})\|` 
and still with :math:`\xi_n = \Frac{\gamma_n}{\alpha(\sigma_{n})}` this gives 
the equation
+Choosing the factor :math:`\alpha(\sigma_{n}) = \|\mbox{Dev}(\sigma_{n})\|` 
and still with :math:`\xi_n = \dfrac{\gamma_n}{\alpha(\sigma_{n})}` this gives 
the equation
 
 .. math::  \varepsilon^p_{n+1} - \varepsilon^p_{n} = (1-\theta)\Delta t \xi_n 
\mbox{Dev}(\sigma_{n}) + \theta \Delta t \xi_{n+1} \mbox{Dev}(\sigma_{n+1}).
 
 Since :math:`\mbox{Dev}(\sigma_{n+1}) = 2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 
2\mu\varepsilon^p_{n+1}` this directly gives:
 
-.. math:: \tilde{\mathscr E}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}) 
= \zeta_n + \left(1-\Frac{1}{1+2\mu\theta\Delta 
t\xi_{n+1}}\right)(\mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n),
+.. math:: \tilde{\mathscr E}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}) 
= \zeta_n + \left(1-\dfrac{1}{1+2\mu\theta\Delta 
t\xi_{n+1}}\right)(\mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n),
 
 which is a linear expression with respect to :math:`u_{n+1}` (but not with 
respect to :math:`\xi_{n+1}`).
 
@@ -303,7 +303,7 @@ Moreover, :math:`\zeta_n` is defined by
 
 One has
 
-.. math:: \|\mbox{Dev}(\sigma_{n+1})\| = 
2\mu\|\mbox{Dev}(\varepsilon(u_{n+1})) -\varepsilon^p_{n+1}\| = 
\Frac{2\mu}{1+2\mu\theta\Delta t \xi_{n+1}}\|\mbox{Dev}(\varepsilon(u_{n+1})) - 
\zeta_n\|,
+.. math:: \|\mbox{Dev}(\sigma_{n+1})\| = 
2\mu\|\mbox{Dev}(\varepsilon(u_{n+1})) -\varepsilon^p_{n+1}\| = 
\dfrac{2\mu}{1+2\mu\theta\Delta t \xi_{n+1}}\|\mbox{Dev}(\varepsilon(u_{n+1})) 
- \zeta_n\|,
 
 Thus, denoting :math:`B = \mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n`, either
 
@@ -311,57 +311,57 @@ Thus, denoting :math:`B = 
\mbox{Dev}(\varepsilon(u_{n+1})) - \zeta_n`, either
 
 and :math:`\xi_{n+1} = 0`, i.e. we are in the elastic case, or  
:math:`\|\mbox{Dev}(\sigma_{n+1})\| =  \sqrt{\frac{2}{3}}` and one obtains
 
-.. math:: 1+2\mu\theta\Delta t \xi_{n+1} = 
\Frac{2\mu\|B\|}{\sqrt{\frac{2}{3}}\sigma_y},
+.. math:: 1+2\mu\theta\Delta t \xi_{n+1} = 
\dfrac{2\mu\|B\|}{\sqrt{\frac{2}{3}}\sigma_y},
 
 and thus
 
-.. math:: \varepsilon^p_{n+1} = \zeta_n + \left( 1 - 
\sqrt{\frac{2}{3}}\Frac{\sigma_y}{2\mu\|B\|}\right) B.
+.. math:: \varepsilon^p_{n+1} = \zeta_n + \left( 1 - 
\sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right) B.
 
 The two options can be summarized by
 
-.. math:: \varepsilon^p_{n+1} = {\mathscr E}^p(u_{n+1}, \zeta_{n}) = \zeta_n + 
\left( 1 - \sqrt{\frac{2}{3}}\Frac{\sigma_y}{2\mu\|B\|}\right)_+ B.
+.. math:: \varepsilon^p_{n+1} = {\mathscr E}^p(u_{n+1}, \zeta_{n}) = \zeta_n + 
\left( 1 - \sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right)_+ B.
 
 The multiplier :math:`\xi_{n+1}` (needed for the :math:`\theta`-scheme for 
:math:`\theta \ne 1`) is given by
 
-.. math:: \xi_{n+1} = \Frac{1}{\theta\Delta 
t}\left(\sqrt{\frac{3}{2}}\Frac{\|B\|}{\sigma_y} - \Frac{1}{2\mu}\right)_+.
+.. math:: \xi_{n+1} = \dfrac{1}{\theta\Delta 
t}\left(\sqrt{\frac{3}{2}}\dfrac{\|B\|}{\sigma_y} - \dfrac{1}{2\mu}\right)_+.
 
 
 **Plane strain approximation**
 
 The plane strain approximation has the same expression replacing the 3D strain 
tensors by the in-plane ones :math:`\bar{\varepsilon}^p` and  
:math:`\bar{\varepsilon}(u_{n+1})`.
 
-.. math:: \bar{\tilde{\mathscr E}}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, 
\bar{\zeta}_{n}) = \bar{\zeta}_n + \left(1-\Frac{1}{1+2\mu\theta\Delta 
t\xi_{n+1}}\right)(\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
\bar{\zeta}_n),
+.. math:: \bar{\tilde{\mathscr E}}^p(u_{n+1}, \theta \Delta t \xi_{n+1}, 
\bar{\zeta}_{n}) = \bar{\zeta}_n + \left(1-\dfrac{1}{1+2\mu\theta\Delta 
t\xi_{n+1}}\right)(\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
\bar{\zeta}_n),
 
-where :math:`\overline{\mbox{Dev}}(\bar{\varepsilon}) = \bar{\varepsilon} - 
\Frac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}` is the 2D restriction of the 3D 
deviator.
+where :math:`\overline{\mbox{Dev}}(\bar{\varepsilon}) = \bar{\varepsilon} - 
\dfrac{\mbox{tr}(\bar{\varepsilon})}{3} \bar{I}` is the 2D restriction of the 
3D deviator.
 
 Moreover, for the yield condition,
 
-.. math:: \|\mbox{Dev}(\sigma)\|^2 = 
4\mu^2\left(\|\overline{\mbox{Dev}}\bar{\varepsilon}(u) - 
\bar{\varepsilon}^p\|^2 + \left(\Frac{\mbox{tr}(\bar{\varepsilon}(u))}{3} 
-\mbox{tr}(\bar{\varepsilon}^p) \right)^2\right).
+.. math:: \|\mbox{Dev}(\sigma)\|^2 = 
4\mu^2\left(\|\overline{\mbox{Dev}}\bar{\varepsilon}(u) - 
\bar{\varepsilon}^p\|^2 + \left(\dfrac{\mbox{tr}(\bar{\varepsilon}(u))}{3} 
-\mbox{tr}(\bar{\varepsilon}^p) \right)^2\right).
 
 And for the elimination of the multiplier,
 
-.. math:: \bar{\mathscr E}^p(\bar{u}_{n+1}, \bar{\varepsilon}^p_{n}) = 
\bar{\zeta}^p_{n} + \left( 1 - 
\sqrt{\frac{2}{3}}\Frac{\sigma_y}{2\mu\|B\|}\right)_+ \bar{B}
+.. math:: \bar{\mathscr E}^p(\bar{u}_{n+1}, \bar{\varepsilon}^p_{n}) = 
\bar{\zeta}^p_{n} + \left( 1 - 
\sqrt{\frac{2}{3}}\dfrac{\sigma_y}{2\mu\|B\|}\right)_+ \bar{B}
 
-with :math:`\bar{B} = 
\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1}))-\bar{\zeta}_{n}` and 
:math:`\|B\|^2 = \|\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
\bar{\zeta}_n\|^2 + \left(\Frac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} 
-\mbox{tr}(\bar{\zeta}_n) \right)^2`.
+with :math:`\bar{B} = 
\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1}))-\bar{\zeta}_{n}` and 
:math:`\|B\|^2 = \|\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
\bar{\zeta}_n\|^2 + \left(\dfrac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} 
-\mbox{tr}(\bar{\zeta}_n) \right)^2`.
 
 **Plane stress approximation**
 
 For plane stress approximation, using :eq:`plane_stress_iso` we deduce from 
the expression of the 3D case
 
-.. math::  \bar{\varepsilon}^p_{n+1} = \Frac{1}{1+2\mu\theta\Delta 
\xi}\left(\bar{\zeta}_{n} +2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) 
- \Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1})) - 
\mbox{tr}(\bar{\varepsilon}_{n+1}^p))\bar{I}\right) \right),
+.. math::  \bar{\varepsilon}^p_{n+1} = \dfrac{1}{1+2\mu\theta\Delta 
\xi}\left(\bar{\zeta}_{n} +2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) 
- \dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1})) - 
\mbox{tr}(\bar{\varepsilon}_{n+1}^p))\bar{I}\right) \right),
 
-since :math:`\mbox{Dev}(\varepsilon(u)) = \varepsilon(u) - 
\Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u)) - 
\mbox{tr}(\bar{\varepsilon}^p))`. Of course, this relation still has to be 
inverted. Denoting :math:`\alpha = 1+2\mu\theta\Delta \xi`, :math:`\beta = 
\Frac{4\mu^2\theta\Delta \xi}{3\lambda+6\mu}` and :math:`C = \bar{\zeta}_{n} 
+2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) - 
\Frac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1}))))\bar{I}\right)`
 one [...]
+since :math:`\mbox{Dev}(\varepsilon(u)) = \varepsilon(u) - 
\dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u)) - 
\mbox{tr}(\bar{\varepsilon}^p))`. Of course, this relation still has to be 
inverted. Denoting :math:`\alpha = 1+2\mu\theta\Delta \xi`, :math:`\beta = 
\dfrac{4\mu^2\theta\Delta \xi}{3\lambda+6\mu}` and :math:`C = \bar{\zeta}_{n} 
+2\mu\theta\Delta \xi\left(\bar{\varepsilon}(u_{n+1}) - 
\dfrac{2\mu}{3(\lambda+2\mu)}(\mbox{tr}(\bar{\varepsilon}(u_{n+1}))))\bar{I}\right)`
  [...]
 
-.. math:: \bar{\varepsilon}^p_{n+1} = \Frac{\beta 
\mbox{tr}(C)}{\alpha(\alpha-2\beta)}\bar{I} + \Frac{1}{\alpha}C.
+.. math:: \bar{\varepsilon}^p_{n+1} = \dfrac{\beta 
\mbox{tr}(C)}{\alpha(\alpha-2\beta)}\bar{I} + \dfrac{1}{\alpha}C.
 
 Moreover, for the yield condition, expression :eq:`plane_stress_dev` can be 
used.
 
 Isotropic elastoplasticity with linear isotropic and kinematic hardening and 
Von-Mises criterion
 
================================================================================================
 
-We consider an isotropic elastic reponse and the internal variable 
:math:`\alpha : \Omega \rightarrow \R` being the accumulated plastic strain 
which satisfies
+We consider an isotropic elastic reponse and the internal variable 
:math:`\alpha : \Omega \rightarrow \rm I\hspace{-0.15em}R` being the 
accumulated plastic strain which satisfies
 
-.. math:: \dot{\alpha} = \sqrt{\Frac{2}{3}}\gamma
+.. math:: \dot{\alpha} = \sqrt{\dfrac{2}{3}}\gamma
 
 For :math:`H_i` the isotropic hardening modulus, the linear hardening consists 
in
 
@@ -378,67 +378,67 @@ for :math:`\sigma_{y0}` the initial uniaxial yield 
stress. The yield function (a
 
 where :math:`H_k` is the kinematic hardening modulus. The same computation as 
in the previous section leads to
 
-.. math:: \tilde{\mathscr E}^p(u_{n+1}, \theta\Delta t \xi_{n+1}, \zeta_n) = 
\zeta_n + \Frac{1}{2(\mu+H_k/3)}\left(1 - \Frac{1}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}\right)(2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n)
+.. math:: \tilde{\mathscr E}^p(u_{n+1}, \theta\Delta t \xi_{n+1}, \zeta_n) = 
\zeta_n + \dfrac{1}{2(\mu+H_k/3)}\left(1 - \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}\right)(2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n)
 
 
-.. math:: \begin{array}{rcl} \tilde{\mathscr A}(u_{n+1}, \theta \Delta t 
\xi_{n+1}, \zeta_{n}, \eta_n) &=& \eta_n + \sqrt{\Frac{2}{3}} \theta \Delta t 
\xi_{n+1}\|\mbox{Dev}(\sigma_{n+1} - \frac{2}{3}H_k\varepsilon^p_{n+1})\| \\ 
&=& \eta_n + \sqrt{\Frac{2}{3}} \theta \Delta t 
\xi_{n+1}\|2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 
2(\mu+H_k/3)\varepsilon^p_{n+1}\| \\ &=&  \eta_n + \sqrt{\Frac{2}{3}} 
\Frac{\theta \Delta t \xi_{n+1}}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}\|2\mu\mbox{Dev}(\varepsi [...]
+.. math:: \begin{array}{rcl} \tilde{\mathscr A}(u_{n+1}, \theta \Delta t 
\xi_{n+1}, \zeta_{n}, \eta_n) &=& \eta_n + \sqrt{\dfrac{2}{3}} \theta \Delta t 
\xi_{n+1}\|\mbox{Dev}(\sigma_{n+1} - \frac{2}{3}H_k\varepsilon^p_{n+1})\| \\ 
&=& \eta_n + \sqrt{\dfrac{2}{3}} \theta \Delta t 
\xi_{n+1}\|2\mu\mbox{Dev}(\varepsilon(u_{n+1})) - 
2(\mu+H_k/3)\varepsilon^p_{n+1}\| \\ &=&  \eta_n + \sqrt{\dfrac{2}{3}} 
\dfrac{\theta \Delta t \xi_{n+1}}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}\|2\mu\mbox{Dev}(\var [...]
 
 where :math:`\zeta_n` and :math:`\eta_n` are defined by
 
 .. math:: \zeta_n = \varepsilon^p_n+(1-\theta)\Delta t \xi_n 
(\mbox{Dev}(\sigma_n)-\frac{2}{3}H_k\varepsilon^n_p) = 
\varepsilon^p_n+(1-\theta)\Delta t \xi_n 
\left(2\mu\mbox{Dev}(\varepsilon(u_{n}))-2(\mu+H_k/3)\varepsilon^n_p\right),
 
-.. math:: \eta_n  = \alpha_n+(1-\theta)\sqrt{\Frac{2}{3}}\Delta t \xi_n 
\|\mbox{Dev}(\sigma_n)-\frac{2}{3}H_k\varepsilon^n_p\| =  
\alpha_n+(1-\theta)\sqrt{\Frac{2}{3}}\Delta t \xi_n 
\|2\mu\mbox{Dev}(\varepsilon(u_{n}))-2(\mu+H_k/3)\varepsilon^n_p\|.
+.. math:: \eta_n  = \alpha_n+(1-\theta)\sqrt{\dfrac{2}{3}}\Delta t \xi_n 
\|\mbox{Dev}(\sigma_n)-\frac{2}{3}H_k\varepsilon^n_p\| =  
\alpha_n+(1-\theta)\sqrt{\dfrac{2}{3}}\Delta t \xi_n 
\|2\mu\mbox{Dev}(\varepsilon(u_{n}))-2(\mu+H_k/3)\varepsilon^n_p\|.
 
 Note that the isotropic hardening modulus do not intervene in 
:math:`\tilde{\mathscr E}^p(u_{n+1}, \theta \Delta \xi, \varepsilon^p_{n})` but 
only in :math:`f(\sigma, A)`.
 
 **Elimination of the multiplier (for the return mapping approach)**
 
-Denoting :math:`\delta = \Frac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}`, 
:math:`\beta = \Frac{1-\delta}{2(\mu+H_k/3)}` and :math:`B = 
2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n` the expression for 
:math:`\varepsilon^p_{n+1}` and :math:`\alpha_{n+1}` becomes
+Denoting :math:`\delta = \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta t\xi_{n+1}}`, 
:math:`\beta = \dfrac{1-\delta}{2(\mu+H_k/3)}` and :math:`B = 
2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n` the expression for 
:math:`\varepsilon^p_{n+1}` and :math:`\alpha_{n+1}` becomes
 
-.. math:: \varepsilon^p_{n+1} = \zeta_n+\beta B, ~~~ \alpha_{n+1} = \eta_n + 
\sqrt{\Frac{2}{3}}\beta \|B\|,
+.. math:: \varepsilon^p_{n+1} = \zeta_n+\beta B, ~~~ \alpha_{n+1} = \eta_n + 
\sqrt{\dfrac{2}{3}}\beta \|B\|,
   :label: hardeningepsalp
 
 and the plastic constraint
 
-.. math:: \delta \|B\| \le \sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).
+.. math:: \delta \|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).
 
 Thus, either we are in the elastic case, i.e. :math:`\xi_{n+1} = 0, \delta = 
1` and
 
-.. math:: \|B\| \le \sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i \eta_n),
+.. math:: \|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \eta_n),
 
-or we are in the plastic case and :math:`\xi_{n+1} > 0, \delta < 1`, 
:math:`\delta \|B\| = \sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1})` and 
:math:`(1-\delta)` solves the equation
+or we are in the plastic case and :math:`\xi_{n+1} > 0, \delta < 1`, 
:math:`\delta \|B\| = \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1})` and 
:math:`(1-\delta)` solves the equation
 
-.. math:: \|B\| - (1-\delta)\|B\| = \sqrt{\Frac{2}{3}}\left(\sigma_{y0}+H_i 
\eta_n + \sqrt{\Frac{2}{3}} \Frac{H_i}{2(\mu+H_k/3)}(1-\delta)\|B\|\right),
+.. math:: \|B\| - (1-\delta)\|B\| = \sqrt{\dfrac{2}{3}}\left(\sigma_{y0}+H_i 
\eta_n + \sqrt{\dfrac{2}{3}} \dfrac{H_i}{2(\mu+H_k/3)}(1-\delta)\|B\|\right),
 
 which leads to
 
-.. math:: 1-\delta = 
\Frac{2(\mu+H_k/3)}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i
 \eta_n) \right)
+.. math:: 1-\delta = 
\dfrac{2(\mu+H_k/3)}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i
 \eta_n) \right)
 
 The two cases can be summarized by
 
-.. math:: \beta = 
\Frac{1}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i
 \eta_n) \right)_+
+.. math:: \beta = 
\dfrac{1}{\|B\|(2\mu+\frac{2}{3}(H_k+H_i))}\left(\|B\|-\sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i
 \eta_n) \right)_+
 
 which directly gives :math:`{\mathscr E}^p(u_{n+1}, \zeta_n, \eta_n)` and 
:math:`{\mathscr A}(u_{n+1}, \zeta_n, \eta_n)` thanks to :eq:`hardeningepsalp`. 
The multiplier :math:`\xi_{n+1}` being given by
 
-.. math:: \xi_{n+1} = \Frac{1}{(2(\mu+H_k/3))\theta\Delta 
t}(\Frac{1}{\delta}-1) = \Frac{1}{\theta\Delta 
t}~\Frac{\beta}{1-2(\mu+H_k/3)\beta}.
+.. math:: \xi_{n+1} = \dfrac{1}{(2(\mu+H_k/3))\theta\Delta 
t}(\dfrac{1}{\delta}-1) = \dfrac{1}{\theta\Delta 
t}~\dfrac{\beta}{1-2(\mu+H_k/3)\beta}.
 
 
 **Plane strain approximation**
 
-Still denoting  :math:`\delta = \Frac{1}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}`, :math:`\beta = \Frac{1-\delta}{2(\mu+H_k/3)}`, :math:`B = 
2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n` and 
:math:`\overline{B} = 
2\mu\overline{Dev}(\bar{\varepsilon}(u_{n+1}))-2(\mu+H_k/3)\bar{\zeta}_n` its 
in-plane part, one has
+Still denoting  :math:`\delta = \dfrac{1}{1+2(\mu+H_k/3)\theta\Delta 
t\xi_{n+1}}`, :math:`\beta = \dfrac{1-\delta}{2(\mu+H_k/3)}`, :math:`B = 
2\mu\mbox{Dev}(\varepsilon(u_{n+1}))-2(\mu+H_k/3)\zeta_n` and 
:math:`\overline{B} = 
2\mu\overline{Dev}(\bar{\varepsilon}(u_{n+1}))-2(\mu+H_k/3)\bar{\zeta}_n` its 
in-plane part, one has
 
 .. math:: \bar{\tilde{\mathscr E}}^p(u_{n+1}, \theta\Delta t \xi_{n+1}, 
\bar{\zeta}_n) = \bar{\zeta}_n + \beta \overline{B},
 
 
-.. math:: \tilde{\mathscr A}(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, 
\eta_n) = \eta_n + \sqrt{\Frac{2}{3}}\beta\|B\|,
+.. math:: \tilde{\mathscr A}(u_{n+1}, \theta \Delta t \xi_{n+1}, \zeta_{n}, 
\eta_n) = \eta_n + \sqrt{\dfrac{2}{3}}\beta\|B\|,
 
 with
 
-.. math:: \|B\|^2 = \|2\mu\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
2(\mu+H_k/3)\bar{\zeta}_n\|^2 + 
\left(2\mu\Frac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} 
-2(\mu+H_k/3)\mbox{tr}(\bar{\zeta}_n) \right)^2.
+.. math:: \|B\|^2 = \|2\mu\overline{\mbox{Dev}}(\bar{\varepsilon}(u_{n+1})) - 
2(\mu+H_k/3)\bar{\zeta}_n\|^2 + 
\left(2\mu\dfrac{\mbox{tr}(\bar{\varepsilon}(u_{n+1}))}{3} 
-2(\mu+H_k/3)\mbox{tr}(\bar{\zeta}_n) \right)^2.
 
 The yield condition still reads
 
-.. math:: \delta \|B\| \le \sqrt{\Frac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).
+.. math:: \delta \|B\| \le \sqrt{\dfrac{2}{3}}(\sigma_{y0}+H_i \alpha_{n+1}).
 
 and for the elimination of the multiplier, :math:`\beta` has the same 
expression as in the previous section adapting the value of :math:`\|B\|`. The 
expressions of :math:`\bar{\zeta}_n` and :math:`\eta_n` have to be adapted 
accoringly.
 
@@ -451,14 +451,14 @@ Souza-Auricchio elastoplasticity law (for shape memory 
alloys)
 
 See for instance [GR-ST2015]_ for the justification of the construction of 
this flow rule. A Von-Mises stress criterion together with an isotropic elastic 
response, no internal variables and a special type of kinematic hardening is 
considered with a constraint :math:`\|\varepsilon^p\| \le c_3`. The plastic 
potential and yield function have the form
 
-.. math:: \Psi(\sigma) = f(\sigma)  = \left\|\mbox{Dev}\left(\sigma - 
c_1\Frac{\varepsilon^p}{\|\varepsilon^p\|} - c_2\varepsilon^p - \delta 
\Frac{\varepsilon^p}{\|\varepsilon^p\|}\right)\right\| - 
\sqrt{\frac{2}{3}}\sigma_{y},
+.. math:: \Psi(\sigma) = f(\sigma)  = \left\|\mbox{Dev}\left(\sigma - 
c_1\dfrac{\varepsilon^p}{\|\varepsilon^p\|} - c_2\varepsilon^p - \delta 
\dfrac{\varepsilon^p}{\|\varepsilon^p\|}\right)\right\| - 
\sqrt{\frac{2}{3}}\sigma_{y},
 
 with the complementarity condition
 
 .. math::
    \delta \ge 0, \|\varepsilon^p\| \le c_3,  \delta (\|\varepsilon^p\|-c_3) = 
0,
 
-where :math:`c_1, c_2 \mbox{ and } c_3` are some physical parameters. Note 
that :math:`\Frac{\varepsilon^p}{\|\varepsilon^p\|}` has to be understood to be 
the whole unit ball for :math:`\varepsilon^p = 0`.
+where :math:`c_1, c_2 \mbox{ and } c_3` are some physical parameters. Note 
that :math:`\dfrac{\varepsilon^p}{\|\varepsilon^p\|}` has to be understood to 
be the whole unit ball for :math:`\varepsilon^p = 0`.
 
 
 to be done ...
@@ -478,7 +478,7 @@ Generic brick
 
 There are two versions of the generic brick. A first one when the plastic 
multiplier is kept as a variable of the problem where the added term is of the 
form:
 
-.. math:: \int_{\Omega} \sigma_{n+1} : \nabla \delta u dx +  \ds \int_{\Omega} 
(\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, A_{n+1}))_+) \delta\xi dx = 0,
+.. math:: \int_{\Omega} \sigma_{n+1} : \nabla \delta u dx +   \int_{\Omega} 
(\xi_{n+1} - (\xi_{n+1} + r f(\sigma_{n+1}, A_{n+1}))_+) \delta\xi dx = 0,
 
 with :math:`r > 0` having a specific value chosen by the brick (in terms of 
the elasticity coefficients), and when the return mapping strategy is selected 
(plastic multiplier is just a data), just the added term:
 
@@ -624,7 +624,7 @@ computes the new stress constraint values supported by the 
material after a load
 
 .. math::
 
-   u^{n+1} \Rightarrow u^n \ \ \ \ \ \textrm{ and } \ \ \ \ \ \sigma^{n+1} 
\Rightarrow \sigma^n
+   u^{n+1} \rm I\hspace{-0.15em}Rightarrow u^n \ \ \ \ \ \textrm{ and } \ \ \ 
\ \ \sigma^{n+1} \rm I\hspace{-0.15em}Rightarrow \sigma^n
 
 Then, :math:`u^n` and :math:`\sigma^n` contains the new values computed and 
one can restart the process.
 
diff --git a/doc/sphinx/source/userdoc/model_solvers.rst 
b/doc/sphinx/source/userdoc/model_solvers.rst
index 0afba8b..e42461a 100644
--- a/doc/sphinx/source/userdoc/model_solvers.rst
+++ b/doc/sphinx/source/userdoc/model_solvers.rst
@@ -39,7 +39,7 @@ where ``1E-7`` is the relative tolerance for the stopping 
criterion, `1` is the
 
   \min\left( \|F(u)\|_1 / \max(L, 10^{-25}) ~, ~~ \|h\|_1 / \max(\|u\|_1, 
10^{-25})\right) < \varepsilon
 
-where :math:`F(u)` is the residual vector, :math:`\|\cdot\|_1` is the 
classical 1-norm in :math:`\R^n`, :math:`h` is the search direction given by 
Newton's algorithm, :math:`L` is the norm of an estimated external loads 
(coming from source term and Dirichlet bricks) and :math:`u` is the current 
state of the searched variable. The maximum taken with :math:`10^{-25}` is to 
avoid pathological cases when :math:`L` and/or :math:`u` are vanishing.
+where :math:`F(u)` is the residual vector, :math:`\|\cdot\|_1` is the 
classical 1-norm in :math:`\rm I\hspace{-0.15em}R^n`, :math:`h` is the search 
direction given by Newton's algorithm, :math:`L` is the norm of an estimated 
external loads (coming from source term and Dirichlet bricks) and :math:`u` is 
the current state of the searched variable. The maximum taken with 
:math:`10^{-25}` is to avoid pathological cases when :math:`L` and/or :math:`u` 
are vanishing.
 
 
 



reply via email to

[Prev in Thread] Current Thread [Next in Thread]