Merge pull request CEED#1602 from CEED/jeremy/solids-simplify

solids - reduce to minimal example set, see Ratel for full

examples/solids/README.md

@@ -72,7 +72,7 @@ As an alternative example exploiting {code}`-dm_plex_box_faces`, we consider a {
 Sides 1 through 6 are rotated around $x$-axis:
 
 ```
-./elasticity -problem FSInitial-NH1 -E 1 -nu 0.3 -num_steps 40 -snes_linesearch_type cp -dm_plex_box_faces 4,4,4 -bc_clamp 1,2,3,4,5,6 -bc_clamp_1_rotate 0,0,1,0,.3 -bc_clamp_2_rotate 0,0,1,0,.3 -bc_clamp_3_rotate 0,0,1,0,.3 -bc_clamp_4_rotate 0,0,1,0,.3 -bc_clamp_5_rotate 0,0,1,0,.3 -bc_clamp_6_rotate 0,0,1,0,.3
+./elasticity -problem FS-NH -E 1 -nu 0.3 -num_steps 40 -snes_linesearch_type cp -dm_plex_box_faces 4,4,4 -bc_clamp 1,2,3,4,5,6 -bc_clamp_1_rotate 0,0,1,0,.3 -bc_clamp_2_rotate 0,0,1,0,.3 -bc_clamp_3_rotate 0,0,1,0,.3 -bc_clamp_4_rotate 0,0,1,0,.3 -bc_clamp_5_rotate 0,0,1,0,.3 -bc_clamp_6_rotate 0,0,1,0,.3
 ```
 
 :::{note}
@@ -103,7 +103,7 @@ The command line options just shown are the minimum requirements to run the mini
   -
 
 * - `-problem`
-  - Problem to solve (`Linear`, `SS-NH`, `FSInitial-NH1`, etc.)
+  - Problem to solve (`Linear`, `FS-NH`, `FS-MR`, etc.)
   - `Linear`
 
 * - `-forcing`

examples/solids/elasticity.c

@@ -25,9 +25,9 @@
 //
 // Sample meshes can be found at https://github.com/jeremylt/ceedSampleMeshes
 //
-//TESTARGS(name="linear elasticity, MMS")                                 -ceed {ceed_resource} -test -degree 3 -nu 0.3 -E 1 -dm_plex_box_faces 3,3,3
-//TESTARGS(name="Neo-Hookean hyperelasticity, initial configuration 1")   -ceed {ceed_resource} -test -problem FSInitial-NH1 -E 2.8 -nu 0.4 -degree 2 -dm_plex_box_faces 2,2,2 -num_steps 1 -bc_clamp 6 -bc_traction 5 -bc_traction_5 0,0,-.5 -expect_final_strain_energy 2.124627916174e-01
-//TESTARGS(name="Mooney-Rivlin hyperelasticity, initial configuration 1") -ceed {ceed_resource} -test -problem FSInitial-MR1 -mu_1 .5 -mu_2 .5 -nu 0.4 -degree 2 -dm_plex_box_faces 2,2,2 -num_steps 1 -bc_clamp 6 -bc_traction 5 -bc_traction_5 0,0,-.5 -expect_final_strain_energy 2.339138880207e-01
+//TESTARGS(name="linear elasticity, MMS")        -ceed {ceed_resource} -test -degree 3 -nu 0.3 -E 1 -dm_plex_box_faces 3,3,3
+//TESTARGS(name="Neo-Hookean hyperelasticity")   -ceed {ceed_resource} -test -problem FS-NH -E 2.8 -nu 0.4 -degree 2 -dm_plex_box_faces 2,2,2 -num_steps 1 -bc_clamp 6 -bc_traction 5 -bc_traction_5 0,0,-.5 -expect_final_strain_energy 2.124627916174e-01
+//TESTARGS(name="Mooney-Rivlin hyperelasticity") -ceed {ceed_resource} -test -problem FS-MR -mu_1 .5 -mu_2 .5 -nu 0.4 -degree 2 -dm_plex_box_faces 2,2,2 -num_steps 1 -bc_clamp 6 -bc_traction 5 -bc_traction_5 0,0,-.5 -expect_final_strain_energy 2.339138880207e-01
 
 /// @file
 /// CEED elasticity example using PETSc with DMPlex

examples/solids/index.md

@@ -526,254 +526,3 @@ In the case where complete linearization is preferred, note the symmetry $\maths
 Along with 6 entries for $\bm S$, this totals 27 entries of overhead compared to computing everything from $\bm F$.
 This compares with 13 entries of overhead for direct storage of $\{ \bm S, \bm C^{-1}, \log J \}$, which is sufficient for the Neo-Hookean model to avoid all but matrix products.
 :::
-
-(problem-hyperelasticity-finite-strain-current-configuration)=
-
-## Hyperelasticity in current configuration
-
-In the preceeding discussion, all equations have been formulated in the initial configuration.
-This may feel convenient in that the computational domain is clearly independent of the solution, but there are some advantages to defining the equations in the current configuration.
-
-1. Body forces (like gravity), traction, and contact are more easily defined in the current configuration.
-2. Mesh quality in the initial configuration can be very bad for large deformation.
-3. The required storage and numerical representation can be smaller in the current configuration.
-
-Most of the benefit in case 3 can be attained solely by moving the Jacobian representation to the current configuration {cite}`davydov2020matrix`, though residual evaluation may also be slightly faster in current configuration.
-There are multiple commuting paths from the nonlinear weak form in initial configuration {eq}`hyperelastic-weak-form-initial` to the Jacobian weak form in current configuration {eq}`jacobian-weak-form-current`.
-One may push forward to the current configuration and then linearize or linearize in initial configuration and then push forward, as summarized below.
-
-$$
-\begin{CD}
-  {\overbrace{\nabla_X \bm{v} \tcolon \bm{FS}}^{\text{Initial Residual}}}
-  @>{\text{push forward}}>{}>
-  {\overbrace{\nabla_x \bm{v} \tcolon \bm{\tau}}^{\text{Current Residual}}} \\
-  @V{\text{linearize}}V{\begin{smallmatrix} \diff\bm F = \nabla_X\diff\bm u \\ \diff\bm S(\diff\bm E) \end{smallmatrix}}V
-  @V{\begin{smallmatrix} \diff\nabla_x\bm v = -\nabla_x\bm v \nabla_x \diff\bm u \\ \diff\bm\tau(\diff\bm\epsilon) \end{smallmatrix}}V{\text{linearize}}V \\
-  {\underbrace{\nabla_X\bm{v}\tcolon \Big(\diff\bm{F}\bm{S} + \bm{F}\diff\bm{S}\Big)}_\text{Initial Jacobian}}
-  @>{\text{push forward}}>{}>
-  {\underbrace{\nabla_x\bm{v}\tcolon \Big(\diff\bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T \Big)}_\text{Current Jacobian}}
-\end{CD}
-$$ (initial-current-linearize)
-
-We will follow both paths for consistency and because both intermediate representations may be useful for implementation.
-
-### Push forward, then linearize
-
-The first term of {eq}`hyperelastic-weak-form-initial` can be rewritten in terms of the symmetric Kirchhoff stress tensor
-$\bm{\tau}=J\bm{\sigma}=\bm{P}\bm{F}^T = \bm F \bm S \bm F^T$ as
-
-$$
-\nabla_X \bm{v} \tcolon \bm{P} = \nabla_X \bm{v} \tcolon \bm{\tau}\bm{F}^{-T} = \nabla_X \bm{v}\bm{F}^{-1} \tcolon \bm{\tau} = \nabla_x \bm{v} \tcolon \bm{\tau}
-$$
-
-therefore, the weak form in terms of $\bm{\tau}$ and $\nabla_x$ with integral over $\Omega_0$ is
-
-$$
-\int_{\Omega_0}{\nabla_x \bm{v} \tcolon \bm{\tau}} \, dV
- - \int_{\Omega_0}{\bm{v} \cdot \rho_0 \bm{g}} \, dV
- - \int_{\partial \Omega_0}{\bm{v}\cdot(\bm{P}\cdot\hat{\bm{N}})} \, dS
- = 0, \quad \forall \bm v \in \mathcal V.
-$$ (hyperelastic-weak-form-current)
-
-#### Linearize in current configuration
-
-To derive a Newton linearization of {eq}`hyperelastic-weak-form-current`, first we define
-
-$$
-\nabla_x \diff \bm{u} = \nabla_X \diff \bm{u} \  \bm{F}^{-1} = \diff \bm{F} \bm{F}^{-1}
-$$ (nabla_xdu)
-
-and $\bm{\tau}$ for Neo-Hookean materials as the push forward of {eq}`neo-hookean-stress`
-
-$$
-\bm{\tau} = \bm{F}\bm{S}\bm{F}^T = \mu (\bm{b} - \bm I_3) + \lambda \log J \bm{I}_3,
-$$ (tau-neo-hookean)
-
-where $\bm{b} = \bm{F} \bm{F}^T$, is the left Cauchy-Green tensor.
-Then by expanding the directional derivative of $\nabla_x \bm{v} \tcolon \bm{\tau}$, we arrive at
-
-$$
-\diff \ (\nabla_x \bm{v} \tcolon \bm{\tau}) = \diff \ (\nabla_x \bm{v})\tcolon \bm{\tau} + \nabla_x \bm{v} \tcolon \diff \bm{\tau} .
-$$ (hyperelastic-linearization-current1)
-
-The first term of {eq}`hyperelastic-linearization-current1` can be written as
-
-$$
-\begin{aligned} \diff \ (\nabla_x \bm{v})\tcolon \bm{\tau} &= \diff \ (\nabla_X \bm{v} \bm{F}^{-1})\tcolon \bm{\tau} = \Big(\underbrace{\nabla_X (\diff \bm{v})}_{0}\bm{F}^{-1} +  \nabla_X \bm{v}\diff \bm{F}^{-1}\Big)\tcolon \bm{\tau}\\   &= \Big(-\nabla_X \bm{v} \bm{F}^{-1}\diff\bm{F}\bm{F}^{-1}\Big)\tcolon \bm{\tau}=\Big(-\nabla_x \bm{v} \diff\bm{F}\bm{F}^{-1}\Big)\tcolon \bm{\tau}\\   &= \Big(-\nabla_x \bm{v} \nabla_x \diff\bm{u} \Big)\tcolon \bm{\tau}= -\nabla_x \bm{v}\tcolon\bm{\tau}(\nabla_x \diff\bm{u})^T \,, \end{aligned}
-$$
-
-where we have used $\diff \bm{F}^{-1}=-\bm{F}^{-1} \diff \bm{F} \bm{F}^{-1}$ and {eq}`nabla_xdu`.
-Using this and {eq}`hyperelastic-linearization-current1` in {eq}`hyperelastic-weak-form-current` yields the weak form in the current configuration
-
-$$
-\int_{\Omega_0} \nabla_x \bm v \tcolon \Big(\diff\bm\tau - \bm\tau (\nabla_x \diff\bm u)^T \Big) = \text{rhs}.
-$$ (jacobian-weak-form-current)
-
-In the following, we will sometimes make use of the incremental strain tensor in the current configuration,
-
-$$
-\diff\bm\epsilon \equiv \frac{1}{2}\Big(\nabla_x \diff\bm{u} + (\nabla_x \diff\bm{u})^T   \Big) .
-$$
-
-:::{dropdown} Deriving $\diff\bm\tau$ for Neo-Hookean material
-To derive a useful expression of $\diff\bm\tau$ for Neo-Hookean materials, we will use the representations
-
-$$
-\begin{aligned}
-\diff \bm{b} &= \diff \bm{F} \bm{F}^T + \bm{F} \diff \bm{F}^T \\
-&= \nabla_x \diff \bm{u} \ \bm{b} + \bm{b} \ (\nabla_x \diff \bm{u})^T \\
-&= (\nabla_x \diff\bm u)(\bm b - \bm I_3) + (\bm b - \bm I_3) (\nabla_x \diff\bm u)^T + 2 \diff\bm\epsilon
-\end{aligned}
-$$
-
-and
-
-$$
-\begin{aligned} \diff\ (\log J) &= \frac{\partial \log J}{\partial \bm{b}}\tcolon \diff \bm{b} = \frac{\partial J}{J\partial \bm{b}}\tcolon \diff \bm{b}=\frac{1}{2}\bm{b}^{-1}\tcolon \diff \bm{b} \\ &= \frac 1 2 \bm b^{-1} \tcolon \Big(\nabla_x \diff\bm u \ \bm b + \bm b (\nabla_x \diff\bm u)^T \Big) \\ &= \trace (\nabla_x \diff\bm u) \\ &= \trace \diff\bm\epsilon . \end{aligned}
-$$
-
-Substituting into {eq}`tau-neo-hookean` gives
-
-$$
-\begin{aligned}
-\diff \bm{\tau} &= \mu \diff \bm{b} + \lambda \trace (\diff\bm\epsilon) \bm I_3 \\
-&= \underbrace{2 \mu \diff\bm\epsilon + \lambda \trace (\diff\bm\epsilon) \bm I_3 - 2\lambda \log J \diff\bm\epsilon}_{\bm F \diff\bm S \bm F^T} \\
-&\quad + (\nabla_x \diff\bm u)\underbrace{\Big( \mu (\bm b - \bm I_3) + \lambda \log J \bm I_3 \Big)}_{\bm\tau} \\
-&\quad + \underbrace{\Big( \mu (\bm b - \bm I_3) + \lambda \log J \bm I_3 \Big)}_{\bm\tau}  (\nabla_x \diff\bm u)^T ,
-\end{aligned}
-$$ (dtau-neo-hookean)
-
-where the final expression has been identified according to
-
-$$
-\diff\bm\tau = \diff\ (\bm F \bm S \bm F^T) = (\nabla_x \diff\bm u) \bm\tau + \bm F \diff\bm S \bm F^T + \bm\tau(\nabla_x \diff\bm u)^T.
-$$
-:::
-
-Collecting terms, we may thus opt to use either of the two forms
-
-$$
-\begin{aligned}
-\diff \bm{\tau} -\bm{\tau}(\nabla_x \diff\bm{u})^T &= (\nabla_x \diff\bm u)\bm\tau + \bm F \diff\bm S \bm F^T \\
-&= (\nabla_x \diff\bm u)\bm\tau + \lambda \trace(\diff\bm\epsilon) \bm I_3 + 2(\mu - \lambda \log J) \diff\bm\epsilon,
-\end{aligned}
-$$ (cur_simp_Jac)
-
-with the last line showing the especially compact representation available for Neo-Hookean materials.
-
-### Linearize, then push forward
-
-We can move the derivatives to the current configuration via
-
-$$
-\nabla_X \bm v \!:\! \diff\bm P = (\nabla_X \bm v) \bm F^{-1} \!:\! \diff \bm P \bm F^T = \nabla_x \bm v \!:\! \diff\bm P \bm F^T
-$$
-
-and expand
-
-$$
-\begin{aligned}
-\diff\bm P \bm F^T &= \diff\bm F \bm S \bm F^T + \bm F \diff\bm S \bm F^T \\
-&= \underbrace{\diff\bm F \bm F^{-1}}_{\nabla_x \diff\bm u} \underbrace{\bm F \bm S \bm F^T}_{\bm\tau} + \bm F \diff\bm S \bm F^T .
-\end{aligned}
-$$
-
-:::{dropdown} Representation of $\bm F \diff\bm S \bm F^T$ for Neo-Hookean materials
-Now we push {eq}`eq-neo-hookean-incremental-stress` forward via
-
-$$
-\begin{aligned}
-\bm F \diff\bm S \bm F^T &= \lambda (\bm C^{-1} \!:\! \diff\bm E) \bm F \bm C^{-1} \bm F^T
-  + 2 (\mu - \lambda \log J) \bm F \bm C^{-1} \diff\bm E \, \bm C^{-1} \bm F^T \\
-    &= \lambda (\bm C^{-1} \!:\! \diff\bm E) \bm I_3 + 2 (\mu - \lambda \log J) \bm F^{-T} \diff\bm E \, \bm F^{-1} \\
-    &= \lambda \operatorname{trace}(\nabla_x \diff\bm u) \bm I_3 + 2 (\mu - \lambda \log J) \diff\bm \epsilon
-\end{aligned}
-$$
-
-where we have used
-
-$$
-\begin{aligned}
-\bm C^{-1} \!:\! \diff\bm E &= \bm F^{-1} \bm F^{-T} \!:\! \bm F^T \diff\bm F \\
-&= \operatorname{trace}(\bm F^{-1} \bm F^{-T} \bm F^T \diff \bm F) \\
-&= \operatorname{trace}(\bm F^{-1} \diff\bm F) \\
-&= \operatorname{trace}(\diff \bm F \bm F^{-1}) \\
-&= \operatorname{trace}(\nabla_x \diff\bm u)
-\end{aligned}
-$$
-
-and
-
-$$
-\begin{aligned}
-\bm F^{-T} \diff\bm E \, \bm F^{-1} &= \frac 1 2 \bm F^{-T} (\bm F^T \diff\bm F + \diff\bm F^T \bm F) \bm F^{-1} \\
-&= \frac 1 2 (\diff \bm F \bm F^{-1} + \bm F^{-T} \diff\bm F^T) \\
-&= \frac 1 2 \Big(\nabla_x \diff\bm u + (\nabla_x\diff\bm u)^T \Big) \equiv \diff\bm\epsilon.
-\end{aligned}
-$$
-:::
-
-Collecting terms, the weak form of the Newton linearization for Neo-Hookean materials in the current configuration is
-
-$$
-\int_{\Omega_0} \nabla_x \bm v \!:\! \Big( (\nabla_x \diff\bm u) \bm\tau + \lambda \operatorname{trace}(\diff\bm\epsilon)\bm I_3 + 2(\mu - \lambda\log J)\diff \bm\epsilon \Big) dV = \text{rhs},
-$$ (jacobian-weak-form-current2)
-
-which equivalent to Algorithm 2 of {cite}`davydov2020matrix` and requires only derivatives with respect to the current configuration. Note that {eq}`cur_simp_Jac` and {eq}`jacobian-weak-form-current2` have recovered the same representation
-using different algebraic manipulations.
-
-:::{tip}
-We define a second order *Green-Euler* strain tensor (cf. Green-Lagrange strain {eq}`eq-green-lagrange-strain`) as
-
-$$
-\bm e = \frac 1 2 \Big(\bm{b} - \bm{I}_3 \Big) = \frac 1 2 \Big( \nabla_X \bm{u} + (\nabla_X \bm{u})^T + \nabla_X \bm{u} \, (\nabla_X \bm{u})^T \Big).
-$$ (green-euler-strain)
-
-Then, the Kirchhoff stress tensor {eq}`tau-neo-hookean` can be written as
-
-$$
-\bm \tau = \lambda \log J \bm I_{3} + 2\mu \bm e,
-$$ (tau-neo-hookean-stable)
-
-which is more numerically stable for small strain, and thus preferred for computation. Note that the $\log J$ is computed via `log1p` {eq}`log1p`, as we discussed in the previous tip.
-:::
-
-### Jacobian representation
-
-We have implemented four storage variants for the Jacobian in our finite strain hyperelasticity. In each case, some variables are computed during residual evaluation and used during Jacobian application.
-
-:::{list-table} Four algorithms for Jacobian action in finite strain hyperelasticity problem
-:header-rows: 1
-:widths: auto
-
-* - Option `-problem`
-  - Static storage
-  - Computed storage
-  - \# scalars
-  - Equations
-
-
-* - `FSInitial-NH1`
-  - $\nabla_{X} \hat X, \operatorname{det}\nabla_{\hat X} X$
-  - $\nabla_X \bm u$
-  - 19
-  - {eq}`eq-diff-P` {eq}`eq-neo-hookean-incremental-stress`
-
-* - `FSInitial-NH2`
-  - $\nabla_{X} \hat X, \operatorname{det}\nabla_{\hat X} X$
-  - $\nabla_X \bm u, \bm C^{-1}, \lambda \log J$
-  - 26
-  - {eq}`eq-diff-P` {eq}`eq-neo-hookean-incremental-stress`
-
-* - `FSCurrent-NH1`
-  - $\nabla_{X} \hat X, \operatorname{det}\nabla_{\hat X} X$
-  - $\nabla_X \bm u$
-  - 19
-  - {eq}`jacobian-weak-form-current` {eq}`nabla_xdu`
-
-* - `FSCurrent-NH2`
-  - $\operatorname{det}\nabla_{\hat X} X$
-  - $\nabla_x \hat X, \bm \tau, \lambda \log J$
-  - 17
-  - {eq}`jacobian-weak-form-current` {eq}`jacobian-weak-form-current2`
-:::

examples/solids/problems/cl-problems.h

@@ -7,22 +7,7 @@
 #pragma once
 
 // Problem options
-typedef enum {
-  ELAS_LINEAR        = 0,
-  ELAS_SS_NH         = 1,
-  ELAS_FSInitial_NH1 = 2,
-  ELAS_FSInitial_NH2 = 3,
-  ELAS_FSCurrent_NH1 = 4,
-  ELAS_FSCurrent_NH2 = 5,
-  ELAS_FSInitial_MR1 = 6
-} problemType;
-static const char *const problemTypes[]        = {"Linear",        "SS-NH",         "FSInitial-NH1", "FSInitial-NH2", "FSCurrent-NH1",
-                                                  "FSCurrent-NH2", "FSInitial-MR1", "problemType",   "ELAS_",         0};
-static const char *const problemTypesForDisp[] = {
-    "Linear elasticity",
-    "Hyperelasticity small strain, Neo-Hookean",
-    "Hyperelasticity finite strain Initial configuration Neo-Hookean w/ dXref_dxinit, Grad(u) storage",
-    "Hyperelasticity finite strain Initial configuration Neo-Hookean w/ dXref_dxinit, Grad(u), C_inv, constant storage",
-    "Hyperelasticity finite strain Current configuration Neo-Hookean w/ dXref_dxinit, Grad(u) storage",
-    "Hyperelasticity finite strain Current configuration Neo-Hookean w/ dXref_dxcurr, tau, constant storage",
-    "Hyperelasticity finite strain Initial configuration Moony-Rivlin w/ dXref_dxinit, Grad(u) storage"};
+typedef enum { ELAS_LINEAR = 0, ELAS_FS_NH = 2, ELAS_FS_MR = 2 } problemType;
+static const char *const problemTypes[]        = {"Linear", "FS-NH", "FS-MR", "problemType", "ELAS_", 0};
+static const char *const problemTypesForDisp[] = {"Linear elasticity", "Hyperelasticity finite strain Initial configuration Neo-Hookean",
+                                                  "Hyperelasticity finite strain Initial configuration Moony-Rivlin"};