Models

fenics-constitutive: Interfaces for constitutive models for dolfinx

__all__ module-attribute

__all__ = ['LinearElasticityModel', 'MisesPlasticityLinearHardening3D', 'SpringKelvinModel', 'SpringMaxwellModel', 'VonMises3D']

IncrSmallStrainModel

Bases: ABC

Interface for incremental small strain models.

Source code in src/fenics_constitutive/models/interfaces.py

class IncrSmallStrainModel(ABC):
    """
    Interface for incremental small strain models.
    """

    @abstractmethod
    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray | None,
        history: dict[str, np.ndarray] | None,
    ) -> None:
        r"""
        Evaluate the constitutive model and overwrite the stress, tangent and history.

        Args:
            t: The current global time $t_n$.
            del_t: The time increment $\Delta t$. The time at the end of the increment is $t_{n+1}=t_n+\Delta t$.
            grad_del_u: The gradient of the increment of the displacement field $\nabla\delta$ with $\delta=u_{n+1}-u_n$.
            stress: The current stress in Mandel notation.
            tangent: The tangent compatible with Mandel notation. `None` if you don't need the tangent like in explicit dynamics.
            history: The history variable(s).
        """

    @property
    @abstractmethod
    def constraint(self) -> StressStrainConstraint:
        """
        The constraint for the stresses or the strains.

        Returns:
            The constraint.
        """

    @property
    def stress_strain_dim(self) -> int:
        """
        The stress-strain dimension that the model is implemented for.

        Returns:
            The stress-strain dimension.
        """
        return self.constraint.stress_strain_dim

    @property
    def geometric_dim(self) -> int:
        """
        The geometric dimension that the model is implemented for.

        Returns:
            The geometric dimension.
        """
        return self.constraint.geometric_dim

    @property
    @abstractmethod
    def history_dim(self) -> dict[str, int | tuple[int, int]] | None:
        """
        The dimensions of history variable(s). This is needed to tell the solver which quadrature
        spaces or arrays to build. If it is not none, a dictionary is returned with the name of the
        history variable as key and the dimension of the history variable as value.

        Returns:
            The dimension of the history variable(s).
        """

constraint abstractmethod property

constraint: StressStrainConstraint

The constraint for the stresses or the strains.

Returns:

Type	Description
StressStrainConstraint	The constraint.

geometric_dim property

geometric_dim: int

The geometric dimension that the model is implemented for.

Returns:

Type	Description
int	The geometric dimension.

history_dim abstractmethod property

history_dim: dict[str, int | tuple[int, int]] | None

The dimensions of history variable(s). This is needed to tell the solver which quadrature spaces or arrays to build. If it is not none, a dictionary is returned with the name of the history variable as key and the dimension of the history variable as value.

Returns:

Type	Description
dict[str, int | tuple[int, int]] | None	The dimension of the history variable(s).

stress_strain_dim property

stress_strain_dim: int

The stress-strain dimension that the model is implemented for.

Returns:

Type	Description
int	The stress-strain dimension.

evaluate abstractmethod

evaluate(t: float, del_t: float, grad_del_u: ndarray, stress: ndarray, tangent: ndarray | None, history: dict[str, ndarray] | None) -> None

Evaluate the constitutive model and overwrite the stress, tangent and history.

Parameters:

Name	Type	Description	Default
t	float	The current global time $t_n$.	required
del_t	float	The time increment $\Delta t$. The time at the end of the increment is $t_{n+1}=t_n+\Delta t$.	required
grad_del_u	ndarray	The gradient of the increment of the displacement field $\nabla\delta$ with $\delta=u_{n+1}-u_n$.	required
stress	ndarray	The current stress in Mandel notation.	required
tangent	ndarray | None	The tangent compatible with Mandel notation. None if you don't need the tangent like in explicit dynamics.	required
history	dict[str, ndarray] | None	The history variable(s).	required

Source code in src/fenics_constitutive/models/interfaces.py

@abstractmethod
def evaluate(
    self,
    t: float,
    del_t: float,
    grad_del_u: np.ndarray,
    stress: np.ndarray,
    tangent: np.ndarray | None,
    history: dict[str, np.ndarray] | None,
) -> None:
    r"""
    Evaluate the constitutive model and overwrite the stress, tangent and history.

    Args:
        t: The current global time $t_n$.
        del_t: The time increment $\Delta t$. The time at the end of the increment is $t_{n+1}=t_n+\Delta t$.
        grad_del_u: The gradient of the increment of the displacement field $\nabla\delta$ with $\delta=u_{n+1}-u_n$.
        stress: The current stress in Mandel notation.
        tangent: The tangent compatible with Mandel notation. `None` if you don't need the tangent like in explicit dynamics.
        history: The history variable(s).
    """

LinearElasticityModel

Bases: IncrSmallStrainModel

A linear elastic material model which has been implemented for all constraints.

Parameters:

Name	Type	Description	Default
parameters	dict[str, float]	Material parameters. Must contain "E" for the Youngs modulus and "nu" for the Poisson ratio.	required
constraint	StressStrainConstraint	Constraint type.	required

Source code in src/fenics_constitutive/models/linear_elasticity_model.py

class LinearElasticityModel(IncrSmallStrainModel):
    """
    A linear elastic material model which has been implemented for all constraints.

    Args:
        parameters: Material parameters. Must contain "E" for the Youngs modulus and "nu" for the Poisson ratio.
        constraint: Constraint type.
    """

    def __init__(
        self, parameters: dict[str, float], constraint: StressStrainConstraint
    ):
        self._constraint = constraint
        E = parameters["E"]
        nu = parameters["nu"]
        self.D = get_elastic_tangent(E, nu, constraint)

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray | None,
        history: np.ndarray | dict[str, np.ndarray] | None,
    ) -> None:
        # Unused: t, del_t, history
        assert (
            grad_del_u.size // (self.geometric_dim**2)
            == stress.size // self.stress_strain_dim
        )
        n_gauss = grad_del_u.size // (self.geometric_dim**2)
        if tangent is not None:
            assert n_gauss == tangent.size // (self.stress_strain_dim**2)
            tangent[:] = np.tile(self.D.flatten(), n_gauss)

        mandel_view = stress.reshape(-1, self.stress_strain_dim)
        strain_increment = strain_from_grad_u(grad_del_u, self.constraint)
        mandel_view += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D

    @property
    def constraint(self) -> StressStrainConstraint:
        return self._constraint

    @property
    def history_dim(self) -> None:
        return None

MisesPlasticityLinearHardening3D

Bases: IncrSmallStrainModel

A von Mises plasticity model with linear hardening for 3D stress states.

This class implements the von Mises yield criterion with linear isotropic hardening. The yield function is defined as: $f = \sqrt{3/2 \cdot s:s} - \sigma_y$, where:

• $s$ is the deviatoric stress tensor
• $\sigma_y = y_0 + h \cdot \alpha$ is the current yield stress
• $\alpha$ is the equivalent plastic strain

Parameters:

Name	Type	Description	Default
parameters	dict[str, ndarray]	A dictionary containing: "mu": Shear modulus "kappa": Bulk modulus "y_0": Initial yield stress "h": Linear hardening modulus	required

Source code in src/fenics_constitutive/models/rust_models.py

@fenics_constitutive_wrapper(PyMisesPlasticity3D)
class MisesPlasticityLinearHardening3D(IncrSmallStrainModel):
    r"""
    A von Mises plasticity model with linear hardening for 3D stress states.

    This class implements the von Mises yield criterion with linear isotropic hardening.
    The yield function is defined as: $f = \sqrt{3/2 \cdot s:s} - \sigma_y$, where:

    - $s$ is the deviatoric stress tensor
    - $\sigma_y = y_0 + h \cdot \alpha$ is the current yield stress
    - $\alpha$ is the equivalent plastic strain

    Args:
        parameters (dict[str, np.ndarray]): A dictionary containing:

            - "mu": Shear modulus
            - "kappa": Bulk modulus
            - "y_0": Initial yield stress
            - "h": Linear hardening modulus
    """

PlaneStrainFrom3D

Bases: IncrSmallStrainModel

Convert a 3D model to a plane strain model. This is achieved by copying the relevant 2D components to a 3D array, calling the 3D model, and then copying the 3D components back to the 2D array. Since this converter creates new arrays for the 3D components, it is (probably) not suitable for large-scale simulations due to the memory consumption.

Parameters:

Name	Type	Description	Default
model	IncrSmallStrainModel	3D model to convert to plane strain.	required

Attributes:

Name	Type	Description
model	IncrSmallStrainModel	3D model to convert to plane strain.
stress_3d	ndarray	3D stress array.
tangent_3d	ndarray	3D tangent array.
grad_del_u_3d	ndarray	3D array of gradient of displacement increment.

Source code in src/fenics_constitutive/models/utils.py

class PlaneStrainFrom3D(IncrSmallStrainModel):
    """
    Convert a 3D model to a plane strain model. This is achieved by copying the
    relevant 2D components to a 3D array, calling the 3D model, and then copying the
    3D components back to the 2D array. Since this converter creates new arrays for
    the 3D components, it is (probably) not suitable for large-scale simulations due
    to the memory consumption.

    Args:
        model: 3D model to convert to plane strain.

    Attributes:
        model (IncrSmallStrainModel): 3D model to convert to plane strain.
        stress_3d (np.ndarray): 3D stress array.
        tangent_3d (np.ndarray): 3D tangent array.
        grad_del_u_3d (np.ndarray): 3D array of gradient of displacement increment.
    """

    def __init__(self, model: IncrSmallStrainModel) -> None:
        assert model.constraint == StressStrainConstraint.FULL
        self.model = model
        self.stress_3d = None
        self.tangent_3d = None
        self.grad_del_u_3d = None

    @property
    def constraint(self) -> StressStrainConstraint:
        return StressStrainConstraint.PLANE_STRAIN

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray | None,
        history: dict[str, np.ndarray] | None,
    ) -> None:
        n_gauss = int(grad_del_u.size / 4)
        current_tangent = None
        if tangent is not None:
            self.tangent_3d = (
                np.zeros(6 * 6 * n_gauss)
                if self.tangent_3d is None
                else self.tangent_3d
            )
            current_tangent = self.tangent_3d

        self.stress_3d = (
            np.zeros(6 * n_gauss) if self.stress_3d is None else self.stress_3d
        )
        self.grad_del_u_3d = (
            np.zeros(9 * n_gauss) if self.grad_del_u_3d is None else self.grad_del_u_3d
        )

        self._grad_del_u_to_3d(grad_del_u)
        self._stress_to_3d(stress)

        self.model.evaluate(
            t, del_t, self.grad_del_u_3d, self.stress_3d, current_tangent, history
        )
        if tangent is not None:
            self._tangent_to_2d(tangent)
        self._stress_to_2d(stress)

    @property
    def history_dim(self) -> dict[str, int | tuple[int, int]] | None:
        return self.model.history_dim

    @df.common.timed("model-conversion-wrapper")
    def _grad_del_u_to_3d(
        self,
        grad_del_u_2d: np.ndarray,
    ) -> None:
        # grad_del_u_2d: 0 1
        #                2 3
        # grad_del_u_3d: 0 1 2
        #                3 4 5
        #                6 7 8
        # We map the 0, 1, 2,3 components of the 2D grad_del_u to the 3D grad_del_u
        # components 0, 1, 3, 4
        self.grad_del_u_3d.reshape(-1, 9)[:, 0:2] = grad_del_u_2d.reshape(-1, 4)[:, 0:2]
        self.grad_del_u_3d.reshape(-1, 9)[:, 3:5] = grad_del_u_2d.reshape(-1, 4)[:, 2:4]

    @df.common.timed("model-conversion-wrapper")
    def _stress_to_3d(
        self,
        stress_2d: np.ndarray,
    ) -> None:
        # Copy the 0, 1, 2, 3 components of the 2D stress to the 3D stress
        self.stress_3d.reshape(-1, 6)[:, 0:4] = stress_2d.reshape(-1, 4)

    @df.common.timed("model-conversion-wrapper")
    def _stress_to_2d(self, stress_2d: np.ndarray) -> None:
        # Copy the 0, 1, 2, 3 components of the 3D stress to the 2D stress
        stress_2d.reshape(-1, 4)[:] = self.stress_3d.reshape(-1, 6)[:, 0:4]

    @df.common.timed("model-conversion-wrapper")
    def _tangent_to_2d(self, tangent_2d: np.ndarray) -> None:
        # tangent_2d: 0 1 2 3
        #             4 5 6 7
        #             8 9 10 11
        #             12 13 14 15
        # tangent_3d: 0 1 2 3 4 5
        #             6 7 8 9 10 11
        #             12 13 14 15 16 17
        #             18 19 20 21 22 23
        #             24 25 26 27 28 29
        #             30 31 32 33 34 35
        # We map the first 4x4 block of the 3D tangent to the 2D tangent
        view_2d = tangent_2d.reshape(-1, 16)
        view_3d = self.tangent_3d.reshape(-1, 36)

        view_2d[:, 0:4] = view_3d[:, 0:4]
        view_2d[:, 4:8] = view_3d[:, 6:10]
        view_2d[:, 8:12] = view_3d[:, 12:16]
        view_2d[:, 12:16] = view_3d[:, 18:22]

SpringKelvinModel

Bases: IncrSmallStrainModel

viscoelastic model based on 1D Three Parameter Model with spring and Kelvin body in row

                       |--- E_1: spring ---|
   --- E_0: spring  ---|                   |--
                       |--- eta: damper ---|

with deviatoric assumptions for 3D generalization (volumetric part of visco strain == 0 damper just working on deviatoric part) time integration: backward Euler

Parameters:

Name	Type	Description	Default
parameters	dict[str, float]	Material parameters. Must contain "E0" for the elastic Youngs modulus, "E1" for the viscous modulus and "tau" for the relaxation time.	required
constraint	StressStrainConstraint	Constraint type.	required

Source code in src/fenics_constitutive/models/spring_kelvin_model.py

class SpringKelvinModel(IncrSmallStrainModel):
    """viscoelastic model based on 1D Three Parameter Model with spring and Kelvin body in row

                               |--- E_1: spring ---|
           --- E_0: spring  ---|                   |--
                               |--- eta: damper ---|

    with deviatoric assumptions for 3D generalization (volumetric part of visco strain == 0 damper just working on deviatoric part)
    time integration: backward Euler

    Args:
        parameters: Material parameters. Must contain "E0" for the elastic Youngs modulus, "E1" for the viscous modulus and "tau" for the relaxation time.
        constraint: Constraint type.
    """

    def __init__(
        self, parameters: dict[str, float], constraint: StressStrainConstraint
    ):
        self._constraint = constraint
        self.E0 = parameters["E0"]  # elastic modulus
        self.E1 = parameters["E1"]  # visco modulus
        self.tau = parameters[
            "tau"
        ]  # relaxation time == eta/(2 mu1) for 1D case eta/E1
        if constraint == StressStrainConstraint.UNIAXIAL_STRESS:
            self.nu = 0.0
        else:
            self.nu = parameters["nu"]  # Poisson's ratio

        self.D_0 = get_elastic_tangent(self.E0, self.nu, constraint)
        self.I2 = get_identity(self.stress_strain_dim, constraint)
        self.mu0, self.lam0 = lame_parameters(self.E0, self.nu)
        self.mu1, _ = lame_parameters(self.E1, self.nu)

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray,
        history: np.ndarray | dict[str, np.ndarray] | None,
    ) -> None:
        _ = t
        assert (
            grad_del_u.size // (self.geometric_dim**2)
            == stress.size // self.stress_strain_dim
            == tangent.size // (self.stress_strain_dim**2)
        )
        n_gauss = grad_del_u.size // (self.geometric_dim**2)
        mandel_view = stress.reshape(-1, self.stress_strain_dim)
        strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(
            -1, self.stress_strain_dim
        )
        if history is None:
            msg = "history must not be None"
            raise ValueError(msg)
        strain_visco_n = history["strain_visco"].reshape(-1, self.stress_strain_dim)
        strain_n = history["strain"].reshape(-1, self.stress_strain_dim)
        I2 = np.tile(self.I2, n_gauss).reshape(-1, self.stress_strain_dim)
        tr_eps = np.sum(strain_increment[:, : self.geometric_dim], axis=1)[
            :, np.newaxis
        ]
        assert del_t > 0, "Time step must be defined and positive."
        factor = 1 / del_t + 1 / self.tau + self.mu0 / (self.tau * self.mu1)
        _deps_visko = (
            1
            / factor
            * (
                1 / (self.tau * 2 * self.mu1) * mandel_view
                - 1 / self.tau * strain_visco_n
                + self.mu0 / (self.tau * self.mu1) * strain_increment
                + self.lam0 / (self.tau * 2 * self.mu1) * tr_eps * I2
            )
        )
        mandel_view += strain_increment @ self.D_0 - 2 * self.mu0 * _deps_visko
        D = (1 - self.mu0 / (self.tau * self.mu1 * factor)) * self.D_0
        tangent[:] = np.tile(D.flatten(), n_gauss)
        strain_visco_n += _deps_visko
        strain_n += strain_increment

    @property
    def constraint(self) -> StressStrainConstraint:
        return self._constraint

    @property
    def history_dim(self) -> dict[str, int | tuple[int, int]] | None:
        return {
            "strain_visco": self.stress_strain_dim,
            "strain": self.stress_strain_dim,
        }

SpringMaxwellModel

Bases: IncrSmallStrainModel

viscoelastic model based on 1D Three Parameter Model with spring and Maxwell body in parallel

     |----------- E_0: spring  ----------|
   --|                                   |--
     |--- E_1: spring --- eta: damper ---|

with deviatoric assumptions for 3D generalization (volumetric part of visco strain == 0 damper just working on deviatoric part) time integration: backward Euler

Parameters:

Name	Type	Description	Default
parameters	dict[str, float]	Material parameters. Must contain "E0" for the elastic Youngs modulus, "E1" for the viscous modulus and "tau" for the relaxation time.	required
constraint	StressStrainConstraint	Constraint type.	required

Source code in src/fenics_constitutive/models/spring_maxwell_model.py

class SpringMaxwellModel(IncrSmallStrainModel):
    """viscoelastic model based on 1D Three Parameter Model with spring and Maxwell body in parallel

             |----------- E_0: spring  ----------|
           --|                                   |--
             |--- E_1: spring --- eta: damper ---|

    with deviatoric assumptions for 3D generalization (volumetric part of visco strain == 0 damper just working on deviatoric part)
    time integration: backward Euler

    Args:
        parameters: Material parameters. Must contain "E0" for the elastic Youngs modulus, "E1" for the viscous modulus and "tau" for the relaxation time.
        constraint: Constraint type.

    """

    def __init__(
        self, parameters: dict[str, float], constraint: StressStrainConstraint
    ):
        self._constraint = constraint
        self.E0 = parameters["E0"]  # elastic modulus
        self.E1 = parameters["E1"]  # visco modulus
        self.tau = parameters["tau"]  # relaxation time
        if constraint == StressStrainConstraint.UNIAXIAL_STRESS:
            self.nu = 0.0
        else:
            self.nu = parameters["nu"]  # Poisson's ratio

        self.D_0 = get_elastic_tangent(self.E0, self.nu, constraint)
        self.D_1 = get_elastic_tangent(self.E1, self.nu, constraint)
        self.mu1, _ = lame_parameters(self.E1, self.nu)

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray,
        history: np.ndarray | dict[str, np.ndarray] | None,
    ) -> None:
        _ = t
        assert (
            grad_del_u.size // (self.geometric_dim**2)
            == stress.size // self.stress_strain_dim
            == tangent.size // (self.stress_strain_dim**2)
        )

        # reshape gauss point arrays
        n_gauss = grad_del_u.size // (self.geometric_dim**2)
        mandel_view = stress.reshape(-1, self.stress_strain_dim)

        strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(
            -1, self.stress_strain_dim
        )
        if history is None:
            msg = "history must not be None"
            raise ValueError(msg)
        strain_visco_n = history["strain_visco"].reshape(-1, self.stress_strain_dim)
        strain_n = history["strain"].reshape(-1, self.stress_strain_dim)

        assert del_t > 0, "Time step must be defined and positive."

        strain_total = strain_n + strain_increment
        factor = 1 / del_t + 1 / self.tau
        _deps_visko = (
            1
            / factor
            * (
                1 / (self.tau * 2 * self.mu1) * strain_total @ self.D_1
                - 1 / self.tau * strain_visco_n
            )
        )

        dstress = strain_increment @ (self.D_0 + self.D_1) - 2 * self.mu1 * _deps_visko
        mandel_view += dstress
        if tangent is not None:
            assert grad_del_u.size // (self.geometric_dim**2) == tangent.size // (
                self.stress_strain_dim**2
            )
            D = self.D_0 + (1 - 1 / (self.tau * factor)) * self.D_1
            tangent[:] = np.tile(D.flatten(), n_gauss)

        strain_visco_n += _deps_visko
        strain_n += strain_increment

    @property
    def constraint(self) -> StressStrainConstraint:
        return self._constraint

    @property
    def history_dim(self) -> dict[str, int | tuple[int, int]] | None:
        return {
            "strain_visco": self.stress_strain_dim,
            "strain": self.stress_strain_dim,
        }

StressStrainConstraint

Bases: Enum

Enum for the model constraint.

The constraint can either be: UNIAXIAL_STRAIN, UNIAXIAL_STRESS, PLANE_STRAIN, PLANE_STRESS, FULL

Source code in src/fenics_constitutive/models/interfaces.py

class StressStrainConstraint(Enum):
    """
    Enum for the model constraint.

    The constraint can either be: UNIAXIAL_STRAIN,
    UNIAXIAL_STRESS, PLANE_STRAIN, PLANE_STRESS, FULL

    """

    UNIAXIAL_STRAIN = 1
    UNIAXIAL_STRESS = 2
    PLANE_STRAIN = 3
    PLANE_STRESS = 4
    FULL = 5

    @property
    def stress_strain_dim(self) -> int:
        """
        The stress-strain dimension of the constraint.

        Returns:
            The stress-strain dimension.
        """
        match self:
            case StressStrainConstraint.UNIAXIAL_STRAIN:
                return 1
            case StressStrainConstraint.UNIAXIAL_STRESS:
                return 1
            case StressStrainConstraint.PLANE_STRAIN:
                return 4
            case StressStrainConstraint.PLANE_STRESS:
                return 4
            case StressStrainConstraint.FULL:
                return 6
            case _:
                msg = f"Constraint {self} not supported"
                raise Exception(msg)

    @property
    def geometric_dim(self) -> int:
        """
        The geometric dimension for the constraint.

        Returns:
            The geometric dimension.
        """
        match self:
            case StressStrainConstraint.UNIAXIAL_STRAIN:
                return 1
            case StressStrainConstraint.UNIAXIAL_STRESS:
                return 1
            case StressStrainConstraint.PLANE_STRAIN:
                return 2
            case StressStrainConstraint.PLANE_STRESS:
                return 2
            case StressStrainConstraint.FULL:
                return 3
            case _:
                msg = f"Constraint {self} not supported"
                raise Exception(msg)

geometric_dim property

geometric_dim: int

The geometric dimension for the constraint.

Returns:

Type	Description
int	The geometric dimension.

stress_strain_dim property

stress_strain_dim: int

The stress-strain dimension of the constraint.

Returns:

Type	Description
int	The stress-strain dimension.

UniaxialStrainFrom3D

Bases: IncrSmallStrainModel

Convert a 3D model to a uniaxial strain model. This is achieved by copying the relevant 1D components to a 3D array, calling the 3D model, and then copying the 3D components back to the 1D array. Since this converter creates new arrays for the 3D components, it is (probably) not suitable for large-scale simulations due to the memory consumption.

Parameters:

Name	Type	Description	Default
model	IncrSmallStrainModel	3D model to convert to uniaxial strain.	required

Attributes:

Name	Type	Description
model	IncrSmallStrainModel	3D model to convert to uniaxial strain.
stress_3d	ndarray	3D stress array.
tangent_3d	ndarray	3D tangent array.
grad_del_u_3d	ndarray	3D array of gradient of displacement increment.

Source code in src/fenics_constitutive/models/utils.py

class UniaxialStrainFrom3D(IncrSmallStrainModel):
    """
    Convert a 3D model to a uniaxial strain model. This is achieved by copying the
    relevant 1D components to a 3D array, calling the 3D model, and then copying the
    3D components back to the 1D array. Since this converter creates new arrays for
    the 3D components, it is (probably) not suitable for large-scale simulations due
    to the memory consumption.

    Args:
        model: 3D model to convert to uniaxial strain.

    Attributes:
        model (IncrSmallStrainModel): 3D model to convert to uniaxial strain.
        stress_3d (np.ndarray): 3D stress array.
        tangent_3d (np.ndarray): 3D tangent array.
        grad_del_u_3d (np.ndarray): 3D array of gradient of displacement increment.
    """

    def __init__(self, model: IncrSmallStrainModel) -> None:
        assert model.constraint == StressStrainConstraint.FULL
        self.model = model
        self.stress_3d = None
        self.tangent_3d = None
        self.grad_del_u_3d = None

    @property
    def constraint(self) -> StressStrainConstraint:
        return StressStrainConstraint.UNIAXIAL_STRAIN

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray | None,
        history: dict[str, np.ndarray] | None,
    ) -> None:
        current_tangent = None
        if tangent is not None:
            self.tangent_3d = (
                np.zeros(6 * 6 * len(grad_del_u))
                if self.tangent_3d is None
                else self.tangent_3d
            )
            current_tangent = self.tangent_3d

        self.stress_3d = (
            np.zeros(6 * len(grad_del_u)) if self.stress_3d is None else self.stress_3d
        )
        self.grad_del_u_3d = (
            np.zeros(9 * len(grad_del_u))
            if self.grad_del_u_3d is None
            else self.grad_del_u_3d
        )

        self._grad_del_u_to_3d(grad_del_u)
        self._stress_to_3d(stress)

        self.model.evaluate(
            t, del_t, self.grad_del_u_3d, self.stress_3d, current_tangent, history
        )
        if tangent is not None:
            self._tangent_to_1d(tangent)
        self._stress_to_1d(stress)

    @property
    def history_dim(self) -> dict[str, int | tuple[int, int]] | None:
        return self.model.history_dim

    @df.common.timed("model-conversion-wrapper")
    def _grad_del_u_to_3d(self, grad_del_u_1d: np.ndarray) -> None:
        # Copy the 11 component of the 1D grad_del_u to the 3D grad_del_u
        self.grad_del_u_3d.reshape(-1, 9)[:, 0] = grad_del_u_1d

    @df.common.timed("constitutive-model-conversion-wrapper")
    def _stress_to_3d(self, stress_1d: np.ndarray) -> None:
        # Copy the 11 component of the 1D stress to the 3D stress
        self.stress_3d.reshape(-1, 6)[:, 0] = stress_1d

    @df.common.timed("model-conversion-wrapper")
    def _stress_to_1d(self, stress_1d: np.ndarray) -> None:
        # Copy the 11 component of the 3D stress to the 1D stress
        stress_1d[:] = self.stress_3d.reshape(-1, 6)[:, 0]

    @df.common.timed("model-conversion-wrapper")
    def _tangent_to_1d(self, tangent_1d: np.ndarray) -> None:
        # Copy the 11 component of the 3D tangent to the 1D tangent
        tangent_1d[:] = self.tangent_3d.reshape(-1, 36)[:, 0]

VonMises3D

Bases: IncrSmallStrainModel

Von Mises Plasticity model with non-linear isotropic hardening. Computation of trial stress state is entirely deviatoric. Volumetric part is added later when the stress increment for the current time step is calculated.

Following are the elastic potential, plastic potential and yield surface accordingly

$$\begin{aligned} & \hat{\psi}_e\left(\varepsilon_e\right) = \frac{1}{2} \kappa {e_e^2}+\mu \varepsilon{_e^{\prime}}: \varepsilon{_e^{\prime}} \\ & \hat{\psi}_p(\alpha)=\left(y_{\infty}-y_0\right)\left(-\frac{1}{\omega}+\alpha+\frac{1}{\omega} \exp (- \omega \alpha)\right) \\ & \hat{\phi}(\boldsymbol{\sigma}, \beta)=\left\|\boldsymbol{\sigma}^{\prime}\right\|-\sqrt{\frac{2} {3}}\left(y_0+\beta\right) \quad \text { with } \quad \beta:=\partial_\alpha \hat{\psi}_p(\alpha) \end{aligned}$$

Parameters:

Name	Type	Description	Default
param	dict[str, float]	Must contain following material parameters: p_ka : bulk modulus, p_mu : shear modulus, p_y0 : initial yield stress, p_y00 : final yield stress, p_w : saturation parameter	required

Source code in src/fenics_constitutive/models/mises_plasticity_isotropic_hardening.py

class VonMises3D(IncrSmallStrainModel):
    r"""
    Von Mises Plasticity model with non-linear isotropic hardening.
    Computation of trial stress state is entirely deviatoric. Volumetric part is added later
    when the stress increment for the current time step is calculated.

    Following are the elastic potential, plastic potential and yield surface accordingly

    $$
    \begin{aligned}
    & \hat{\psi}_e\left(\varepsilon_e\right) = \frac{1}{2} \kappa {e_e^2}+\mu \varepsilon{_e^{\prime}}: 
    \varepsilon{_e^{\prime}} \\
    & \hat{\psi}_p(\alpha)=\left(y_{\infty}-y_0\right)\left(-\frac{1}{\omega}+\alpha+\frac{1}{\omega} \exp (- 
    \omega \alpha)\right) \\
    & \hat{\phi}(\boldsymbol{\sigma}, \beta)=\left\|\boldsymbol{\sigma}^{\prime}\right\|-\sqrt{\frac{2} 
    {3}}\left(y_0+\beta\right) \quad \text { with } \quad \beta:=\partial_\alpha \hat{\psi}_p(\alpha)
    \end{aligned}
    $$

    Args:
           param: Must contain following material parameters: p_ka :  bulk modulus, p_mu : shear modulus, p_y0 : initial yield stress, p_y00 : final yield stress, p_w : saturation parameter
    """

    def __init__(self, param: dict[str, float]):
        self.xioi = np.array(
            [
                [1, 1, 1, 0, 0, 0],  # 1dyadic1 tensor
                [1, 1, 1, 0, 0, 0],
                [1, 1, 1, 0, 0, 0],
                [0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0, 0],
            ]
        )

        self.I2 = get_identity(self.stress_strain_dim, self.constraint)

        self.I4 = np.eye(
            self.stress_strain_dim, dtype=np.float64
        )  # Identity of rank 4 tensor
        self.xpp = self.I4 - (1 / 3) * self.xioi  # Projection tensor of rank 4

        self.p_ka = param["p_ka"]  # kappa - bulk modulus
        self.p_mu = param["p_mu"]  # mu - shear modulus
        self.p_y0 = param["p_y0"]  # y0 - initial yield stress
        self.p_y00 = param["p_y00"]  # y00 - final yield stress
        self.p_w = param["p_w"]  # w - saturation parameter

    def evaluate(
        self,
        t: float,
        del_t: float,
        grad_del_u: np.ndarray,
        stress: np.ndarray,
        tangent: np.ndarray | None,
        history: np.ndarray | dict[str, np.ndarray] | None = None,
    ) -> None:
        stress_view = stress.reshape(-1, self.stress_strain_dim)
        strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(
            -1, self.stress_strain_dim
        )
        eps_n = history["eps_n"].reshape(-1, self.stress_strain_dim)
        alpha = history["alpha"]

        for n, eps in enumerate(strain_increment):
            tr_eps = np.sum(eps[:3])  # trace of strain
            eps_dev = eps - tr_eps * self.I2 / 3  # deviatoric strain

            # deviatoric trial internal forces (trial stress)
            del_sigtr = 2 * self.p_mu * eps_dev  # incremental trial stress
            stress_n_dev = (
                stress_view[n] - np.sum(stress_view[n][:3]) * self.I2 / 3
            )  # total deviatoric stress in last time step
            sigtr = (
                stress_n_dev + del_sigtr
            )  # total (deviatoric) trial stress in current time step

            # norm of deviatoric trial internal forces
            sigtrn = np.sqrt(np.dot(sigtr, sigtr))

            # trial yield criterion
            phitr = sigtrn - np.sqrt(2 / 3) * (
                self.p_y0
                + (self.p_y00 - self.p_y0) * (1 - np.exp(-self.p_w * alpha[n]))
            )

            # CHECK YIELD CRITERION
            # elastic - plastic step
            if phitr > 0:
                # initialization
                gamma_0 = 1
                gamma_1 = 0
                xr = 1
                it = 0
                tol = 1e-12
                tol_rel = 1e-8
                nmax = 100
                # flow direction
                xn = sigtr / sigtrn

                # UPDATE PLASTIC MULTIPLIER VIA NEWTON-RAPHSON SCHEME
                def f(x):
                    return (
                        sigtrn
                        - 2 * self.p_mu * x
                        - np.sqrt(2 / 3)
                        * (
                            self.p_y0
                            + (self.p_y00 - self.p_y0)
                            * (1 - np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x)))
                        )
                    )

                def df(x):
                    return -2 * self.p_mu - (2 / 3) * (
                        self.p_y00 - self.p_y0
                    ) * self.p_w * np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x))

                # start Newton iteration
                while np.abs(xr) > tol and abs(gamma_1 - gamma_0) > tol_rel * abs(
                    gamma_1
                ):
                    gamma_0 = gamma_1
                    it = it + 1
                    # compute residium
                    xr = f(gamma_0)
                    # compute tangent
                    xg = df(gamma_0)
                    # update plastic flow
                    gamma_1 = gamma_0 - xr / xg
                    # exit Newton algorithm for iteration > nmax
                    if it > nmax:
                        msg = "Newton-Raphson method did not converge for plastic multiplier."
                        raise RuntimeError(msg)
                    # end of Newton iterration

                # compute tangent with converged gamma
                xg = df(gamma_1)

                # algorithmic parameters
                xc1 = -1 / xg
                xc2 = gamma_1 / sigtrn

                # ELASTIC STEP
            else:
                xn = np.zeros(6)
                gamma_1 = 0
                xc1 = 0
                xc2 = 0

            # update eps^p_n+1 and alpha
            eps_n[n] += gamma_1 * xn
            alpha[n] += np.sqrt(2 / 3) * gamma_1

            # determine incremental elastic-plastic stresses (with volumetric part)
            sh = self.p_ka * tr_eps * self.I2 + del_sigtr - 2 * self.p_mu * gamma_1 * xn
            # update total stresses
            stress_view[n] += sh

            if tangent is not None:
                tangent_view = tangent.reshape(-1, self.stress_strain_dim**2)
                # determine elastic-plastic moduli
                aah = (
                    self.p_ka * self.xioi
                    + 2 * self.p_mu * (1 - 2 * self.p_mu * xc2) * self.xpp
                    + 4 * self.p_mu * self.p_mu * (xc2 - xc1) * np.outer(xn, xn)
                )
                tangent_view[n] = aah.flatten()

    # def update(self) -> None:
    #    pass

    @property
    def constraint(self) -> StressStrainConstraint:
        return StressStrainConstraint.FULL

    @property
    def history_dim(self) -> int:
        return {"eps_n": self.constraint.stress_strain_dim, "alpha": 1}

get_elastic_tangent

get_elastic_tangent(E: float, nu: float, constraint: StressStrainConstraint) -> np.ndarray

Get the linear elastic tangent based on the stress-strain constraint. Args: E (float): Young's modulus. nu (float): Poisson's ratio. constraint (StressStrainConstraint): Stress-strain constraint (.FULL, PLANE_STRAIN, PLANE_STRESS, UNIAXIAL_STRAIN, UNIAXIAL_STRESS).

Returns:

Type	Description
ndarray	np.ndarray: Elastic tangent matrix.

Source code in src/fenics_constitutive/models/utils.py

def get_elastic_tangent(
    E: float, nu: float, constraint: StressStrainConstraint
) -> np.ndarray:
    """Get the linear elastic tangent based on the stress-strain constraint.
    Args:
        E (float): Young's modulus.
        nu (float): Poisson's ratio.
        constraint (StressStrainConstraint): Stress-strain constraint (.FULL, PLANE_STRAIN, PLANE_STRESS, UNIAXIAL_STRAIN, UNIAXIAL_STRESS).

    Returns:
        np.ndarray: Elastic tangent matrix.
    """

    mu, lam = lame_parameters(E, nu)
    match constraint:
        case StressStrainConstraint.FULL:
            # see https://en.wikipedia.org/wiki/Hooke%27s_law
            D = np.array(
                [
                    [2.0 * mu + lam, lam, lam, 0.0, 0.0, 0.0],
                    [lam, 2.0 * mu + lam, lam, 0.0, 0.0, 0.0],
                    [lam, lam, 2.0 * mu + lam, 0.0, 0.0, 0.0],
                    [0.0, 0.0, 0.0, 2.0 * mu, 0.0, 0.0],
                    [0.0, 0.0, 0.0, 0.0, 2.0 * mu, 0.0],
                    [0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * mu],
                ]
            )
        case StressStrainConstraint.PLANE_STRAIN:
            # We assert that the strain is being provided with 0 in the z-direction
            # see https://en.wikipedia.org/wiki/Hooke%27s_law
            D = np.array(
                [
                    [2.0 * mu + lam, lam, lam, 0.0],
                    [lam, 2.0 * mu + lam, lam, 0.0],
                    [lam, lam, 2.0 * mu + lam, 0.0],
                    [0.0, 0.0, 0.0, 2.0 * mu],
                ]
            )
        case StressStrainConstraint.PLANE_STRESS:
            # We do not make any assumptions about strain in the z-direction
            # This matrix just multiplies the z component by 0.0 which results
            # in a plane stress state
            # see https://en.wikipedia.org/wiki/Hooke%27s_law
            D = (
                E
                / (1 - nu**2.0)
                * np.array(
                    [
                        [1.0, nu, 0.0, 0.0],
                        [nu, 1.0, 0.0, 0.0],
                        [0.0, 0.0, 0.0, 0.0],
                        [0.0, 0.0, 0.0, (1.0 - nu)],
                    ]
                )
            )
        case StressStrainConstraint.UNIAXIAL_STRAIN:
            # see https://csmbrannon.net/2012/08/02/distinction-between-uniaxial-stress-and-uniaxial-strain/
            C = E * (1.0 - nu) / ((1.0 + nu) * (1.0 - 2.0 * nu))
            D = np.array([[C]])

        case StressStrainConstraint.UNIAXIAL_STRESS:
            # see https://csmbrannon.net/2012/08/02/distinction-between-uniaxial-stress-and-uniaxial-strain/
            D = np.array([[E]])

        case _:
            msg = "Constraint not implemented"
            raise NotImplementedError(msg)

    return D

get_identity

get_identity(stress_strain_dim: int, constraint: StressStrainConstraint) -> np.ndarray

Get the identity tensor based on the stress-strain constraint. Args: stress_strain_dim (int): Dimension of the stress/strain vector. constraint (StressStrainConstraint): Stress-strain constraint (.FULL, PLANE_STRAIN, PLANE_STRESS, UNIAXIAL_STRAIN, UNIAXIAL_STRESS) Returns: np.ndarray: Identity tensor.

Source code in src/fenics_constitutive/models/utils.py

def get_identity(
    stress_strain_dim: int, constraint: StressStrainConstraint
) -> np.ndarray:
    """Get the identity tensor based on the stress-strain constraint.
    Args:
        stress_strain_dim (int): Dimension of the stress/strain vector.
        constraint (StressStrainConstraint): Stress-strain constraint (.FULL, PLANE_STRAIN, PLANE_STRESS, UNIAXIAL_STRAIN, UNIAXIAL_STRESS)
    Returns:
        np.ndarray: Identity tensor.
    """

    match constraint:
        case StressStrainConstraint.FULL:
            I2 = np.zeros(stress_strain_dim, dtype=np.float64)
            I2[0:3] = 1.0
        case StressStrainConstraint.PLANE_STRAIN:
            # We assert that the strain is being provided with 0 in the z-direction
            I2 = np.zeros(stress_strain_dim, dtype=np.float64)
            I2[0:3] = 1.0
        case StressStrainConstraint.PLANE_STRESS:
            I2 = np.zeros(stress_strain_dim, dtype=np.float64)
            I2[0:2] = 1.0
        case StressStrainConstraint.UNIAXIAL_STRAIN:
            I2 = np.zeros(stress_strain_dim, dtype=np.float64)
            I2[0] = 1.0
        case StressStrainConstraint.UNIAXIAL_STRESS:
            I2 = np.zeros(stress_strain_dim, dtype=np.float64)
            I2[0] = 1.0

        case _:
            msg = "Constraint not implemented"
            raise NotImplementedError(msg)

    return I2

lame_parameters

lame_parameters(E: float, nu: float) -> tuple[float, float]

Compute Lame parameters (mu, lam) from Young's modulus E and Poisson's ratio nu.

Source code in src/fenics_constitutive/models/utils.py

def lame_parameters(E: float, nu: float) -> tuple[float, float]:
    """Compute Lame parameters (mu, lam) from Young's modulus E and Poisson's ratio nu."""
    mu = E / (2.0 * (1.0 + nu))
    lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
    return mu, lam

strain_from_grad_u

strain_from_grad_u(grad_u: ndarray, constraint: StressStrainConstraint) -> np.ndarray

Compute the strain in Mandel notation from the gradient of displacement (or increments of both quantities).

Parameters:

Name	Type	Description	Default
grad_u	ndarray	Gradient of displacement field.	required
constraint	StressStrainConstraint	Constraint that the model is implemented for.	required

Returns:

Type	Description
ndarray	Numpy array containing the strain for all IPs.

Source code in src/fenics_constitutive/models/utils.py

def strain_from_grad_u(
    grad_u: np.ndarray, constraint: StressStrainConstraint
) -> np.ndarray:
    """
    Compute the strain in Mandel notation from the gradient of displacement (or increments of both
    quantities).

    Args:
        grad_u: Gradient of displacement field.
        constraint: Constraint that the model is implemented for.

    Returns:
        Numpy array containing the strain for all IPs.
    """
    strain_dim = constraint.stress_strain_dim
    n_gauss = int(grad_u.size / (constraint.geometric_dim**2))
    strain = np.zeros(strain_dim * n_gauss)
    grad_u_view = grad_u.reshape(-1, constraint.geometric_dim**2)
    strain_view = strain.reshape(-1, strain_dim)

    match constraint:
        case StressStrainConstraint.UNIAXIAL_STRAIN:
            strain_view[:, 0] = grad_u_view[:, 0]
        case StressStrainConstraint.UNIAXIAL_STRESS:
            strain_view[:, 0] = grad_u_view[:, 0]
        case StressStrainConstraint.PLANE_STRAIN:
            """
            Full tensor:

            0 1
            2 3

            Mandel vector:
            f = 1 / sqrt(2)
            0 3 "0" f*(1+2)
            """
            strain_view[:, 0] = grad_u_view[:, 0]
            strain_view[:, 1] = grad_u_view[:, 3]
            strain_view[:, 2] = 0.0
            strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 2])
        case StressStrainConstraint.PLANE_STRESS:
            """
            Full tensor:

            0 1
            2 3

            Mandel vector:
            f = 1 / sqrt(2)
            0 3 "0" f*(1+2)
            """
            strain_view[:, 0] = grad_u_view[:, 0]
            strain_view[:, 1] = grad_u_view[:, 3]
            strain_view[:, 2] = 0.0
            strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 2])
        case StressStrainConstraint.FULL:
            """
            Full tensor:

            0 1 2
            3 4 5
            6 7 8

            Mandel vector:
            f = 1 / sqrt(2)
            0 4 8 f*(1+3) f*(5+7) f*(2+6)
            """
            strain_view[:, 0] = grad_u_view[:, 0]
            strain_view[:, 1] = grad_u_view[:, 4]
            strain_view[:, 2] = grad_u_view[:, 8]
            strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 3])
            strain_view[:, 4] = 1 / 2**0.5 * (grad_u_view[:, 2] + grad_u_view[:, 6])
            strain_view[:, 5] = 1 / 2**0.5 * (grad_u_view[:, 5] + grad_u_view[:, 7])
        case _:
            msg = "Constraint not supported."
            raise NotImplementedError(msg)
    return strain