Skip to content

Commit

Permalink
Few dfixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@bleyerj
bleyerj committed Oct 2, 2024
1 parent 77c40f4 commit f55d4ed
Show file tree
Hide file tree
Showing 3 changed files with 116 additions and 76 deletions.
2 changes: 1 addition & 1 deletion demos/jax/elastoplasticity/Hosford_yield_surface.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ elastic_model = LinearElasticIsotropic(E, nu)
a = 10.0
sig0 = 500.0
material = GeneralIsotropicHardening(
elastic_model, lambda sig: sig0, lambda sig: Hosford_stress(sig, a=a)
elastic_model, lambda p: sig0, lambda sig: Hosford_stress(sig, a=a)
)
```

Expand Down
45 changes: 29 additions & 16 deletions demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,11 +62,14 @@ H0.setEntryName("HardeningSlope");
};
```

## `FEniCS` implementation
## `FEniCSx` implementation

We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight.

<img src="finite_strain_plasticity_solution.png" width="500">
```{image} finite_strain_plasticity_solution.png
:align: center
:width: 500px
```

```{code-cell} ipython3
import numpy as np
Expand All @@ -84,12 +87,16 @@ from dolfinx_materials.utils import (
nonsymmetric_tensor_to_vector,
)
comm = MPI.COMM_WORLD
rank = comm.rank
current_path = os.getcwd()
length, width, height = 1.0, 0.04, 0.1
nx, ny, nz = 30, 5, 8
nx, ny, nz = 20, 4, 6
domain = mesh.create_box(
MPI.COMM_WORLD,
comm,
[(0, -width / 2, -height / 2.0), (length, width / 2, height / 2.0)],
[nx, ny, nz],
cell_type=mesh.CellType.tetrahedron,
Expand Down Expand Up @@ -134,10 +141,11 @@ material = MFrontMaterial(
"LogarithmicStrainPlasticity",
material_properties={"YieldStrength": 250e6, "HardeningSlope": 1e6},
)
print(material.behaviour.getBehaviourType())
print(material.behaviour.getKinematic())
print(material.gradient_names, material.gradient_sizes)
print(material.flux_names, material.flux_sizes)
if rank == 0:
print(material.behaviour.getBehaviourType())
print(material.behaviour.getKinematic())
print(material.gradient_names, material.gradient_sizes)
print(material.flux_names, material.flux_sizes)
```

In this large-strain setting, the `QuadratureMapping` acts from the deformation gradient $\boldsymbol{F}=\boldsymbol{I}+\nabla\boldsymbol{u}$ to the first Piola-Kirchhoff stress $\boldsymbol{P}$. We must therefore register the deformation gradient as `Identity(3)+grad(u)`.
Expand Down Expand Up @@ -178,7 +186,7 @@ Finally, we setup the nonlinear problem, the corresponding Newton solver and sol
problem = NonlinearMaterialProblem(qmap, Res, Jac, u, bcs)
newton = NewtonSolver(MPI.COMM_WORLD)
newton = NewtonSolver(comm)
newton.rtol = 1e-4
newton.atol = 1e-4
newton.convergence_criterion = "residual"
Expand All @@ -203,15 +211,19 @@ for i, t in enumerate(load_steps[1:]):
converged, it = problem.solve(newton, print_solution=False)
print(f"Increment {i+1} converged in {it} iterations.")
if rank == 0:
print(f"Increment {i+1} converged in {it} iterations.")
p0 = qmap.project_on("EquivalentPlasticStrain", ("DG", 0))
vtk.write_function(u, t)
vtk.write_function(p0, t)
w = u.sub(2).collapse()
results[i + 1, 0] = max(np.abs(w.vector.array))
local_max = max(np.abs(w.vector.array))
# Perform the reduction to get the global maximum on rank 0
global_max = comm.reduce(local_max, op=MPI.MAX, root=0)
results[i + 1, 0] = global_max
results[i + 1, 1] = t
vtk.close()
```
Expand All @@ -220,11 +232,12 @@ During the load incrementation, we monitor the evolution of the maximum vertical
The load-displacement curve exhibits a classical elastoplastic behavior rapidly followed by a stiffening behavior due to membrane catenary effects.

```{code-cell} ipython3
plt.figure()
plt.plot(results[:, 0], results[:, 1], "-oC3")
plt.xlabel("Displacement")
plt.ylabel("Load")
plt.show()
if rank==0:
plt.figure()
plt.plot(results[:, 0], results[:, 1], "-oC3")
plt.xlabel("Displacement")
plt.ylabel("Load")
plt.show()
```

## References
Expand Down
145 changes: 86 additions & 59 deletions demos/mfront/finite_strain_elastoplasticity/finite_strain_elastoplasticity.py
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,38 @@
# %% [markdown]
# ---
# jupyter:
# jupytext:
# text_representation:
# extension: .py
# format_name: light
# format_version: '1.5'
# jupytext_version: 1.16.1
# kernelspec:
# display_name: Python 3
# language: python
# name: python3
# ---

# # Finite-strain elastoplasticity within the logarithmic strain framework
#
# This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in [@miehe_anisotropic_2002].
# This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in {cite:p}`miehe2002anisotropic`.
#
# ```{tip}
# This simulation is a bit heavy to run so we suggest running it in parallel.
# ```
#
# ## Logarithmic strains
#
# This framework expresses constitutive relations between the Hencky strain measure $\boldsymbol{H} = \dfrac{1}{2}\log (\boldsymbol{F}^T\cdot\boldsymbol{F})$ and its dual stress measure $\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split additively into many contributions (elastic, plastic, thermal, swelling, etc.) e.g. $\boldsymbol{H}=\boldsymbol{H}^e+\boldsymbol{H}^p$. Its trace is also linked with the volume change $J=\exp(\operatorname{tr}(\boldsymbol{H}))$. As a result, the deformation gradient $\boldsymbol{F}$ is used for expressing the Hencky strain $\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\boldsymbol{H},\boldsymbol{T})$-pair and the dual stress $\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\boldsymbol{\sigma}$ or Piola-Kirchhoff stresses.
# This framework expresses constitutive relations between the Hencky strain measure $\boldsymbol{H} = \dfrac{1}{2}\log (\boldsymbol{F}^\text{T}\cdot\boldsymbol{F})$ and its dual stress measure $\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split additively into many contributions (elastic, plastic, thermal, swelling, etc.) e.g. $\boldsymbol{H}=\boldsymbol{H}^e+\boldsymbol{H}^p$. Its trace is also linked with the volume change $J=\exp(\operatorname{tr}(\boldsymbol{H}))$. As a result, the deformation gradient $\boldsymbol{F}$ is used for expressing the Hencky strain $\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\boldsymbol{H},\boldsymbol{T})$-pair and the dual stress $\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\boldsymbol{\sigma}$ or Piola-Kirchhoff stresses.
#
# ## `MFront` implementation
#
# The logarithmic strain framework discussed in the previous paragraph consists merely as a pre-processing and a post-processing stages of the behaviour integration. The pre-processing stage compute the logarithmic strain and its increment and the post-processing stage inteprets the stress resulting from the behaviour integration as the dual stress $\boldsymbol{T}$ and convert it to the Cauchy stress.
# The logarithmic strain framework discussed in the previous paragraph consists merely as a pre-processing and a post-processing stages of the behavior integration. The pre-processing stage compute the logarithmic strain and its increment and the post-processing stage inteprets the stress resulting from the behavior integration as the dual stress $\boldsymbol{T}$ and convert it to the Cauchy stress.
#
# `MFront` provides the `@StrainMeasure` keyword that allows to specify which strain measure is used by the behaviour. When choosing the `Hencky` strain measure, `MFront` automatically generates those pre- and post-processing stages, allowing the user to focus on the behaviour integration.
# `MFront` provides the `@StrainMeasure` keyword that allows to specify which strain measure is used by the behavior. When choosing the `Hencky` strain measure, `MFront` automatically generates those pre- and post-processing stages, allowing the user to focus on the behavior integration.
#
# This leads to the following implementation (see the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html) for details about the various implementations available):
#
# ```cpp
# ```
# @DSL Implicit;
#
# @Behaviour LogarithmicStrainPlasticity;
Expand Down Expand Up @@ -45,15 +62,16 @@
# };
# ```
#
# ## `FEniCS` implementation
# ## `FEniCSx` implementation
#
# We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight.
#
#
# <img src="finite_strain_plasticity_solution.png" width="500">
#
# ```{image} finite_strain_plasticity_solution.png
# :align: center
# :width: 500px
# ```

# %%
# +
import numpy as np
import matplotlib.pyplot as plt
import os
Expand All @@ -69,12 +87,16 @@
nonsymmetric_tensor_to_vector,
)


comm = MPI.COMM_WORLD
rank = comm.rank

current_path = os.getcwd()

length, width, height = 1.0, 0.04, 0.1
nx, ny, nz = 30, 5, 8
nx, ny, nz = 20, 4, 6
domain = mesh.create_box(
MPI.COMM_WORLD,
comm,
[(0, -width / 2, -height / 2.0), (length, width / 2, height / 2.0)],
[nx, ny, nz],
cell_type=mesh.CellType.tetrahedron,
Expand Down Expand Up @@ -109,26 +131,25 @@ def right(x):
du = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
u = fem.Function(V, name="Displacement")
# -

# %% [markdown]
# The `MFrontMaterial` instance is loaded from the `MFront` `LogarithmicStrainPlasticity` behaviour. This behaviour is a finite-strain behaviour (`material.is_finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\boldsymbol{F}$. By default, a `MFront` behaviour always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piola-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the "flux" associated to the "gradient" $\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [`MFront` conventions on tensors](http://tfel.sourceforge.net/tensors.html).
# The `MFrontMaterial` instance is loaded from the `MFront` `LogarithmicStrainPlasticity` behavior. This behavior is a finite-strain behavior (`material.is_finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\boldsymbol{F}$. By default, a `MFront` behavior always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piola-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the "flux" associated to the "gradient" $\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [`MFront` conventions on tensors](https://thelfer.github.io/tfel/web/tensors.html).

# %%
material = MFrontMaterial(
os.path.join(current_path, "src/libBehaviour.so"),
"LogarithmicStrainPlasticity",
material_properties={"YieldStrength": 250e6, "HardeningSlope": 1e6},
)
print(material.behaviour.getBehaviourType())
print(material.behaviour.getKinematic())
print(material.gradient_names, material.gradient_sizes)
print(material.flux_names, material.flux_sizes)
if rank == 0:
print(material.behaviour.getBehaviourType())
print(material.behaviour.getKinematic())
print(material.gradient_names, material.gradient_sizes)
print(material.flux_names, material.flux_sizes)

# %% [markdown]
# The `QuadratureMap` object must therefore register the deformation gradient as `Identity(3)+grad(u)`.

# In this large-strain setting, the `QuadratureMapping` acts from the deformation gradient $\boldsymbol{F}=\boldsymbol{I}+\nabla\boldsymbol{u}$ to the first Piola-Kirchhoff stress $\boldsymbol{P}$. We must therefore register the deformation gradient as `Identity(3)+grad(u)`.

# %%
# +
def F(u):
return nonsymmetric_tensor_to_vector(ufl.Identity(gdim) + ufl.grad(u))

Expand All @@ -139,75 +160,81 @@ def dF(u):

qmap = QuadratureMap(domain, 2, material)
qmap.register_gradient("DeformationGradient", F(u))
# -

# %%
sig = qmap.fluxes["FirstPiolaKirchhoffStress"]
Res = (ufl.dot(sig, dF(v)) - ufl.dot(selfweight, v)) * qmap.dx
# We will work in a Total Lagrangian formulation, writing the weak form of equilibrium on the reference configuration $\Omega_0$, thereby defining the nonlinear residual weak form as:
# Find $\boldsymbol{u}\in V$ such that:
#
# $$
# \int_{\Omega_0} \boldsymbol{P}(\boldsymbol{F}(\boldsymbol{u})):\nabla \boldsymbol{v} \,\text{d}\Omega - \int_{\Omega_0} \boldsymbol{f}\cdot\boldsymbol{v}\,\text{d}\Omega = 0 \quad \forall \boldsymbol{v}\in V
# $$
# where $\boldsymbol{f}$ is the self-weight.
#
# The corresponding Jacobian form is computed via automatic differentiation. As for the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behavior is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the `MFront` behavior file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.

PK1 = qmap.fluxes["FirstPiolaKirchhoffStress"]
Res = (ufl.dot(PK1, dF(v)) - ufl.dot(selfweight, v)) * qmap.dx
Jac = qmap.derivative(Res, u, du)

# %% [markdown]
# The loading is then defined and, as for the [small-strain elastoplasticity example](https://thelfer.github.io/mgis/web/mgis_fenics_small_strain_elastoplasticity.html), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behaviour is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the `MFront` behaviour file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.
#
# Finally, we setup a few parameters of the Newton non-linear solver.
# Finally, we setup the nonlinear problem, the corresponding Newton solver and solve the load-stepping problem.

# %%
# + tags=["hide-output"]
problem = NonlinearMaterialProblem(qmap, Res, Jac, u, bcs)

newton = NewtonSolver(MPI.COMM_WORLD)
newton = NewtonSolver(comm)
newton.rtol = 1e-4
newton.atol = 1e-4
newton.convergence_criterion = "incremental"
newton.convergence_criterion = "residual"
newton.report = True

# Set solver options
ksp = newton.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

Nincr = 30
load_steps = np.linspace(0.0, 1.0, Nincr + 1)

vtk_u = io.VTKFile(domain.comm, "results_u.pvd", "w")
vtk_p = io.VTKFile(domain.comm, "results_p.pvd", "w")
vtk = io.VTKFile(domain.comm, f"results/{material.name}.pvd", "w")
results = np.zeros((Nincr + 1, 2))
for i, t in enumerate(load_steps[1:]):
selfweight.value[-1] = -50e6 * t

converged, it = problem.solve(newton, print_solution=False)

print(f"Increment {i+1} converged in {it} iterations.")
if rank == 0:
print(f"Increment {i+1} converged in {it} iterations.")

p0 = qmap.project_on("EquivalentPlasticStrain", ("DG", 0))

vtk_u.write_function(u, t)
vtk_p.write_function(p0, t)
vtk.write_function(u, t)
vtk.write_function(p0, t)

w = u.sub(2).collapse()
results[i + 1, 0] = max(np.abs(w.vector.array))
local_max = max(np.abs(w.vector.array))
# Perform the reduction to get the global maximum on rank 0
global_max = MPI.COMM_WORLD.reduce(local_max, op=MPI.MAX, root=0)
results[i + 1, 0] = global_max
results[i + 1, 1] = t
vtk_u.close()
vtk_p.close()
vtk.close()
# -

# %% [markdown]
# During the load incrementation, we monitor the evolution of the maximum vertical downwards displacement.
#
# This simulation is a bit heavy to run so we suggest running it in parallel:
# ```bash
# mpirun -np 4 python3 finite_strain_elastoplasticity.py
# ```

# %% [markdown]
# The load-displacement curve exhibits a classical elastoplastic behaviour rapidly followed by a stiffening behaviour due to membrane catenary effects.
# The load-displacement curve exhibits a classical elastoplastic behavior rapidly followed by a stiffening behavior due to membrane catenary effects.

# %%
plt.figure()
plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel("Displacement")
plt.ylabel("Load")
plt.show()
if rank==0:
plt.figure()
plt.plot(results[:, 0], results[:, 1], "-oC3")
plt.xlabel("Displacement")
plt.ylabel("Load")
plt.show()

# %%
# ## References
#
# ```{bibliography}
# :filter: docname in docnames
# ```

0 comments on commit f55d4ed

Please sign in to comment.