Cleanup and docs for constitutive models

docs/src/api_reference.md

@@ -22,6 +22,12 @@ CKIMaterial
 NoCorrection
 EnergySurfaceCorrection
 ZEMSilling
+LinearElastic
+NeoHooke
+MooneyRivlin
+SaintVenantKirchhoff
+linear_kernel
+cubic_b_spline_kernel
 ```
 
 ## Discretization

src/Peridynamics.jl

@@ -12,7 +12,10 @@ import LibGit2, Dates
 export BBMaterial, OSBMaterial, CMaterial, BACMaterial, CKIMaterial
 
 # CMaterial related types
-export NeoHookeNonlinear, SaintVenantKirchhoff, ZEMSilling
+export LinearElastic, NeoHooke, MooneyRivlin, SaintVenantKirchhoff, ZEMSilling
+
+# Kernels
+export linear_kernel, cubic_b_spline_kernel
 
 # Systems related types
 export NoCorrection, EnergySurfaceCorrection

src/discretization/kernels.jl

@@ -1,8 +1,35 @@
-# peridynamic kernels / influence functions for arbitrary bond systems
+@doc raw"""
+    linear_kernel(δ, L)
 
+A linear kernel function ``\omega`` (also called influence function) used for weighting the
+bonds in a family. The kernel function is defined as
+```math
+\omega = \frac{\delta}{L} \; ,
+```
+with the horizon ``\delta`` and the initial bond length ``L``.
+"""
 @inline linear_kernel(δ, L) = δ / L
 
-@inline function cubic_b_spline(δ, L)
+@doc raw"""
+    cubic_b_spline_kernel(δ, L)
+
+A cubic B-spline kernel function ``\omega`` used for weighting the bonds in a family.
+The kernel function is defined as
+```math
+\begin{aligned}
+\xi &= \frac{L}{\delta} \; , \\
+\omega &= \left\{
+    \begin{array}{ll}
+        \frac{2}{3} - 4 \xi^2 + 4 \xi^3 & \quad \text{if} \; 0 < \xi \leq 0.5 \; , \\[3pt]
+        \frac{4}{3} - 4 \xi + 4 \xi^2 - \frac{4}{3} \xi^3 & \quad \text{if} \; 0.5 < \xi \leq 1 \; , \\[3pt]
+        0 & \quad \text{else} \; ,
+    \end{array}
+    \right.
+\end{aligned}
+```
+with the horizon ``\delta`` and the initial bond length ``L``.
+"""
+@inline function cubic_b_spline_kernel(δ, L)
     ξ = L / δ
     if 0 < ξ ≤ 0.5
         return 2/3 - 4 * ξ^2 + 4 * ξ^3

src/physics/constitutive_models.jl

@@ -1,3 +1,19 @@
+@doc raw"""
+    LinearElastic
+
+Linear elastic constitutive model that can be specified when using a [`CMaterial`](@ref) and
+[`BACMaterial`](@ref).
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\boldsymbol{P} = \mathbb{C} : \boldsymbol{E} \; ,
+```
+with the elastic stiffness tensor  ``\mathbb{C}`` and the Green-Lagrange strain tensor
+``\boldsymbol{E}`` with
+```math
+\boldsymbol{E} = \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I}
+                             \right) \; .
+```
+"""
 struct LinearElastic <: AbstractConstitutiveModel end
 
 function first_piola_kirchhoff(::LinearElastic, storage::AbstractStorage,
@@ -19,6 +35,25 @@ function get_hooke_matrix(nu, λ, μ)
     return CVoigt
 end
 
+@doc raw"""
+    NeoHooke
+
+Neo-Hookean constitutive model that can be specified when using a [`CMaterial`](@ref) and
+[`BACMaterial`](@ref).
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\begin{aligned}
+\boldsymbol{C} &= \boldsymbol{F}^{\top} \boldsymbol{F} \; , \\
+\boldsymbol{S} &= \mu \left( \boldsymbol{I} - \boldsymbol{C}^{-1} \right)
+    + \lambda \log(J) \boldsymbol{C}^{-1} \; , \\
+\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
+\end{aligned}
+```
+with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor
+``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second
+Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and second Lamé parameters
+``\lambda`` and ``\mu``.
+"""
 struct NeoHooke <: AbstractConstitutiveModel end
 
 function first_piola_kirchhoff(::NeoHooke, storage::AbstractStorage,
@@ -30,9 +65,33 @@ function first_piola_kirchhoff(::NeoHooke, storage::AbstractStorage,
     return P
 end
 
-struct NeoHookeNonlinear <: AbstractConstitutiveModel end
+@doc raw"""
+    MooneyRivlin
+
+Mooney-Rivlin constitutive model that can be specified when using a [`CMaterial`](@ref) and
+[`BACMaterial`](@ref).
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\begin{aligned}
+\boldsymbol{C} &= \boldsymbol{F}^{\top} \boldsymbol{F} \; , \\
+\boldsymbol{S} &= G \left( \boldsymbol{I} - \frac{1}{3} \mathrm{tr}(\boldsymbol{C})
+                           \boldsymbol{C}^{-1} \right) \cdot J^{-\frac{2}{3}}
+                + \frac{K}{4} \left( J^2 - J^{-2} \right) \boldsymbol{C}^{-1} \; , \\
+\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
+\end{aligned}
+```
+with the deformation gradient ``\boldsymbol{F}``, the right Cauchy-Green deformation tensor
+``\boldsymbol{C}``, the Jacobian ``J = \mathrm{det}(\boldsymbol{F})``, the second
+Piola-Kirchhoff stress ``\boldsymbol{S}``, the shear modulus ``G``, and the
+bulk modulus ``K``.
 
-function first_piola_kirchhoff(::NeoHookeNonlinear, storage::AbstractStorage,
+# Error handling
+If the Jacobian ``J`` is smaller than the machine precision `eps()` or a `NaN`, the first
+Piola-Kirchhoff stress tensor is defined as ``\boldsymbol{P} = \boldsymbol{0}``.
+"""
+struct MooneyRivlin <: AbstractConstitutiveModel end
+
+function first_piola_kirchhoff(::MooneyRivlin, storage::AbstractStorage,
                                params::AbstractPointParameters, F::SMatrix{3,3,T,9}) where T
     J = det(F)
     J < eps() && return zero(SMatrix{3,3,T,9})
@@ -45,6 +104,25 @@ function first_piola_kirchhoff(::NeoHookeNonlinear, storage::AbstractStorage,
     return P
 end
 
+@doc raw"""
+    SaintVenantKirchhoff
+
+Saint-Venant-Kirchhoff constitutive model that can be specified when using a
+[`CMaterial`](@ref) and [`BACMaterial`](@ref).
+The first Piola-Kirchhoff stress ``\boldsymbol{P}`` is given by
+```math
+\begin{aligned}
+\boldsymbol{E} &= \frac{1}{2} \left( \boldsymbol{F}^{\top} \boldsymbol{F} - \boldsymbol{I}
+                              \right) \; , \\
+\boldsymbol{S} &= \lambda \, \mathrm{tr}(\boldsymbol{E}) \, \boldsymbol{I}
+                + 2 \mu \boldsymbol{E} \; , \\
+\boldsymbol{P} &= \boldsymbol{F} \, \boldsymbol{S} \; ,
+\end{aligned}
+```
+with the deformation gradient ``\boldsymbol{F}``, the Green-Lagrange strain tensor
+``\boldsymbol{E}``, the second Piola-Kirchhoff stress ``\boldsymbol{S}``, and the first and
+second Lamé parameters ``\lambda`` and ``\mu``.
+"""
 struct SaintVenantKirchhoff <: AbstractConstitutiveModel end
 
 function first_piola_kirchhoff(::SaintVenantKirchhoff, storage::AbstractStorage,

src/physics/correspondence.jl

@@ -5,16 +5,24 @@ A material type used to assign the material of a [`Body`](@ref) with the local c
 consistent (correspondence) formulation of non-ordinary state-based peridynamics.
 
 # Keywords
-- `kernel::Function`: Kernel function used for weighting the interactions between points.
-    (default: `linear_kernel`)
-- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior.
-    (default: `LinearElastic()`)
+- `kernel::Function`: Kernel function used for weighting the interactions between points. \\
+    (default: `linear_kernel`) \\
+    The following kernels can be used:
+    - [`linear_kernel`](@ref)
+    - [`cubic_b_spline_kernel`](@ref)
+- `model::AbstractConstitutiveModel`: Constitutive model defining the material behavior. \\
+    (default: `LinearElastic()`) \\
+    The following models can be used:
+    - [`LinearElastic`](@ref)
+    - [`NeoHooke`](@ref)
+    - [`MooneyRivlin`](@ref)
+    - [`SaintVenantKirchhoff`](@ref)
 - `zem::AbstractZEMStabilization`: Zero-energy mode stabilization. The
-    stabilization algorithm of Silling (2017) is used as default.
+    stabilization algorithm of Silling (2017) is used as default. \\
     (default: `ZEMSilling()`)
 - `maxdmg::Float64`: Maximum value of damage a point is allowed to obtain. If this value is
     exceeded, all bonds of that point are broken because the deformation gradient would then
-    possibly contain `NaN` values.
+    possibly contain `NaN` values. \\
     (default: `0.85`)
 
 !!! note "Stability of fracture simulations"

test/integration/b_int_correspondence.jl

@@ -4,7 +4,7 @@
                     0.0 0.0 0.0 1.0 2.0]
     volume = fill(1.0, 5)
     δ = 1.5
-    body = Body(CMaterial(model=NeoHookeNonlinear()), ref_position, volume)
+    body = Body(CMaterial(model=MooneyRivlin()), ref_position, volume)
     material!(body, horizon=δ, rho=1, E=1, nu=0.25, Gc=1.0)
     failure_permit!(body, false)
 
@@ -36,7 +36,7 @@ end
                     0.0 0.0 0.0 1.0 2.0]
     volume = fill(1.0, 5)
     δ = 1.5
-    body = Body(CMaterial(model=NeoHookeNonlinear()), ref_position, volume)
+    body = Body(CMaterial(model=MooneyRivlin()), ref_position, volume)
     point_set!(body, :a, [1])
     point_set!(body, :b, [2,3,4,5])
     material!(body, :a, horizon=δ, rho=1, E=1, nu=0.25, Gc=1.0)

test/integration/symmetry_correspondence.jl

@@ -8,7 +8,7 @@
     pos = hcat(([x;y;z] for x in grid for y in grid for z in grid)...)
     n_points = size(pos, 2)
     vol = fill(Δx^3, n_points)
-    body = Body(CMaterial(model=NeoHookeNonlinear()), pos, vol)
+    body = Body(CMaterial(model=MooneyRivlin()), pos, vol)
     failure_permit!(body, false)
     material!(body, horizon=3.015Δx, rho=7850, E=210e9, nu=0.25, Gc=1)
     point_set!(z -> z > width/2 - 0.6Δx, body, :set_a)