Skip to content

Solver

fenics-constitutive: Interfaces for solver of own constitutive models following models interface for dolfinx

__all__ module-attribute 

__all__ = ['CDMSolver', 'CorotationalIncrSmallStrainProblem', 'IncrSmallStrainProblem', 'critical_timestep', 'diagonal_inverted_mass']

CDMSolver

Class to solve the incremental small strain problem using the central difference method. This class will also determine the critical timestep . The user may however prescribe a timestep in the IncrSmallStrainProblem to be used instead. The solver will then use

with being the safety factor.

Parameters:

Name Type Description Default

problem

IncrSmallStrainProblem

The incremental small strain problem to be solved.

 required

density

list[float] | float

The density of the material. This is either a single value for a homogenous domain or a list of values for each submesh.

 required

u0

Function

The initial displacement. This is used to set the initial condition for the solver.

 required

v0

Function

The initial velocity. This is used to set the initial condition for the solver.

 required

safety_factor

float

The safety factor for the time step. The time step is set to the critical time step multiplied by the safety factor. This should be a value between 0 and 1.

 required
Source code in src/fenics_constitutive/solver/central_difference_method.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
class CDMSolver:
    r"""
    Class to solve the incremental small strain problem using the central difference method.
    This class will also determine the critical timestep $\Delta t_\mathrm{crit}$. The user may however prescribe a
    timestep in the `IncrSmallStrainProblem` $\Delta t_\mathrm{problem}$ to be used instead. The solver will then use

    $$
    \Delta t =\min\{s\cdot\Delta t_\mathrm{crit}, \Delta t_\mathrm{problem}\}
    $$
    with $s$ being the safety factor.

    Args:
        problem: The incremental small strain problem to be solved.
        density: The density of the material. This is either a single 
            value for a homogenous domain or a list of values for each submesh.
        u0: The initial displacement. This is used to set the initial condition for the solver.
        v0: The initial velocity. This is used to set the initial condition for the solver.
        safety_factor: The safety factor for the time step. The time step is set to the critical 
            time step multiplied by the safety factor. This should be a value between 0 and 1.
    """

    def __init__(
        self,
        problem: IncrSmallStrainProblem,
        density: list[float] | float,
        u0: df.fem.Function,
        v0: df.fem.Function,
        safety_factor: float,
    ) -> None:
        self.problem = problem

        density = [density] if isinstance(density, float) else density

        laws = [(law.law, law.cells) for law in problem._law_on_submeshs]

        assert len(density) == len(laws)

        self.del_t_crit = critical_timestep(laws, density, u0)

        self.problem.sim_time.set_timestep(self.del_t_crit.min() * safety_factor)

        self.problem.incr_disp.current.x.array[:] = u0.x.array[:]
        self.problem.incr_disp.current.x.scatter_forward()
        self.problem.incr_disp.previous.x.array[:] = u0.x.array[:]
        self.problem.incr_disp.previous.x.scatter_forward()

        self.f = cast(df.fem.Function, df.fem.Function(v0.function_space))
        self.a = cast(df.fem.Function, df.fem.Function(v0.function_space))
        self.v = v0
        cells = [law[1] for law in laws]
        self.M_inv = diagonal_inverted_mass(
            v0.function_space, density, cells
        )

    def step(self) -> None:
        r"""Advance the solution by $\Delta t$"""

        df.fem.assemble_vector(self.f.x.array, self.problem.L)
        self.f.x.scatter_reverse(df.la.InsertMode.add)

        self.a.x.array[:] = self.M_inv.x.array * self.f.x.array
        self.a.x.scatter_forward()

        self.v.x.array[:] += self.problem.sim_time.dt * self.a.x.array
        self.v.x.scatter_forward()

        self.problem.incr_disp.current.x.array[:] += self.problem.sim_time.dt * self.v.x.array
        for bc in self.problem.bcs:
            bc.set(self.problem.incr_disp.current.x.array)
        self.problem.incr_disp.current.x.scatter_forward()

        for bc in self.problem.bcs:
            bc.set(self.v.x.array, self.problem.incr_disp.previous.x.array, 1.0/self.problem.sim_time.dt)
        self.v.x.scatter_forward()

        self.problem.form_without_petsc(evaluate_tangent=False)

step 

step() -> None

Advance the solution by

Source code in src/fenics_constitutive/solver/central_difference_method.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def step(self) -> None:
    r"""Advance the solution by $\Delta t$"""

    df.fem.assemble_vector(self.f.x.array, self.problem.L)
    self.f.x.scatter_reverse(df.la.InsertMode.add)

    self.a.x.array[:] = self.M_inv.x.array * self.f.x.array
    self.a.x.scatter_forward()

    self.v.x.array[:] += self.problem.sim_time.dt * self.a.x.array
    self.v.x.scatter_forward()

    self.problem.incr_disp.current.x.array[:] += self.problem.sim_time.dt * self.v.x.array
    for bc in self.problem.bcs:
        bc.set(self.problem.incr_disp.current.x.array)
    self.problem.incr_disp.current.x.scatter_forward()

    for bc in self.problem.bcs:
        bc.set(self.v.x.array, self.problem.incr_disp.previous.x.array, 1.0/self.problem.sim_time.dt)
    self.v.x.scatter_forward()

    self.problem.form_without_petsc(evaluate_tangent=False)

CorotationalIncrSmallStrainProblem

Bases: IncrSmallStrainProblem

Corotational extension of IncrSmallStrainProblem that adds objective stress-rate rotation and mesh-update steps.

This subclass reuses all core functionality of the base solver and overrides only the initialization and form method.

Parameters:

Name Type Description Default

laws

A list of tuples where the first element is the constitutive law and the second element is the cells for the submesh. If only one law is provided, it is assumed that the domain is homogenous. The cell indices should be local indices of the MPI-process.

 required

u

The displacement field. This is the unknown in the nonlinear problem.

 required

bcs

The Dirichlet boundary conditions.

 required

q_degree

The quadrature degree (Polynomial degree which the quadrature rule needs to integrate exactly).

 required

del_t

The maximal allowed time increment between steps. This is relevant if an explicit solver is used or if the material model depends on the time, e.g. viscoelasticity. The actually used timestep may be overwritten by the solver if convergence issues occur or if the critical timestep is smaller.

 required

external_forces

The external forces applied to the system. This should be a list of ufl Forms which are subtracted from the residual. The user can use this to apply body forces or Neumann boundary conditions.

 required

form_compiler_options

The options for the form compiler.

 required

jit_options

The options for the JIT compiler.

 required
Source code in src/fenics_constitutive/solver/corotational_solver.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
class CorotationalIncrSmallStrainProblem(IncrSmallStrainProblem):
    """
    Corotational extension of IncrSmallStrainProblem that adds objective
    stress-rate rotation and mesh-update steps.

    This subclass reuses all core functionality of the base solver and
    overrides only the initialization and `form` method.

    Args:
        laws: A list of tuples where the first element is the constitutive law and the second
            element is the cells for the submesh. If only one law is provided, it is assumed
            that the domain is homogenous. The cell indices should be local indices of the MPI-process.
        u: The displacement field. This is the unknown in the nonlinear problem.
        bcs: The Dirichlet boundary conditions.
        q_degree: The quadrature degree (Polynomial degree which the quadrature rule needs to integrate exactly).
        del_t: The maximal allowed time increment between steps. This is relevant if an explicit solver is used or if the 
            material model depends on the time, e.g. viscoelasticity. The actually used timestep may be overwritten by
            the solver if convergence issues occur or if the critical timestep is smaller.
        external_forces: The external forces applied to the system. This should be a list of ufl Forms which are subtracted from the residual. 
            The user can use this to apply body forces or Neumann boundary conditions.
        form_compiler_options: The options for the form compiler.
        jit_options: The options for the JIT compiler.
    """

    def __init__(self, *args, **kwargs) -> None:

        super().__init__(*args, **kwargs)

        self.mesh_updater = MeshUpdater(self.incr_disp)
        self._check_spatial_dimension_3d()

        # Replace each base LawOnSubMesh with its corotational variant
        # while preserving all submesh data and function references
        self._law_on_submeshs = [
            CorotationalLawOnSubMesh(
                law=lom.law,
                cells=lom.cells,
                displacement_gradient_fn=lom.displacement_gradient_fn,
                stress=lom.stress,
                local_tangent=lom.local_tangent,
                submesh_map=lom.submesh_map,
                history=lom.history,
            )
            for lom in self._law_on_submeshs
        ]

    @df.common.timed("constitutive-form-evaluation-corotational")
    def form(self, x: PETSc.Vec) -> None:
        """
        Assemble the nonlinear form using a corotational midpoint configuration.

        This override performs a mesh update to the midpoint configuration
        and applies a corotational stress rotation before evaluating the constitutive
        laws. The mesh is updated to final configuration prior to scattering stress and tangent fields.
        """
        NonlinearProblem.form(self, x)

        self.mesh_updater.move_to_midpoint()

        self.incr_disp.update_current(x)

        for law in self._law_on_submeshs:
            law.evaluate(self.sim_time, self.incr_disp, self.stress, self.tangent)

        self.mesh_updater.move_to_final()

        self.stress.scatter_current()
        self.tangent.x.scatter_forward()


    def _check_spatial_dimension_3d(self) -> None:
        """
        Ensure that the solver is only used for 3D problems.

        The current corotational implementation assumes 3x3 tensors and a
        6-component Mandel stress (3 normal + 3 shear). For 2D / 1D
        constitutive models a separate implementation is required.
        """
        gdim = self._u.function_space.mesh.geometry.dim
        if gdim != 3:
            msg = (
                f"CorotationalIncrSmallStrainProblem currently supports only 3D "
                f"geometries (gdim=3). Got gdim={gdim}."
            )
            raise NotImplementedError(
                msg
            )

form 

form(x: Vec) -> None

Assemble the nonlinear form using a corotational midpoint configuration.

This override performs a mesh update to the midpoint configuration and applies a corotational stress rotation before evaluating the constitutive laws. The mesh is updated to final configuration prior to scattering stress and tangent fields.

Source code in src/fenics_constitutive/solver/corotational_solver.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
@df.common.timed("constitutive-form-evaluation-corotational")
def form(self, x: PETSc.Vec) -> None:
    """
    Assemble the nonlinear form using a corotational midpoint configuration.

    This override performs a mesh update to the midpoint configuration
    and applies a corotational stress rotation before evaluating the constitutive
    laws. The mesh is updated to final configuration prior to scattering stress and tangent fields.
    """
    NonlinearProblem.form(self, x)

    self.mesh_updater.move_to_midpoint()

    self.incr_disp.update_current(x)

    for law in self._law_on_submeshs:
        law.evaluate(self.sim_time, self.incr_disp, self.stress, self.tangent)

    self.mesh_updater.move_to_final()

    self.stress.scatter_current()
    self.tangent.x.scatter_forward()

IncrSmallStrainProblem

Bases: NonlinearProblem

A nonlinear problem for incremental small strain models. To be used with the dolfinx NewtonSolver.

Parameters:

Name Type Description Default

laws

list[tuple[IncrSmallStrainModel, ndarray]] | IncrSmallStrainModel

A list of tuples where the first element is the constitutive law and the second element is the cells for the submesh. If only one law is provided, it is assumed that the domain is homogenous. The cell indices should be local indices of the MPI-process.

 required

u

Function

The displacement field. This is the unknown in the nonlinear problem.

 required

bcs

list[DirichletBC]

The Dirichlet boundary conditions.

 required

q_degree

int

The quadrature degree (Polynomial degree which the quadrature rule needs to integrate exactly).

 required

del_t

float

The maximal allowed time increment between steps. This is relevant if an explicit solver is used or if the material model depends on the time, e.g. viscoelasticity. The actually used timestep may be overwritten by the solver if convergence issues occur or if the critical timestep is smaller.

 1.0

external_forces

list[Form] | None

The external forces applied to the system. This should be a list of ufl Forms which are subtracted from the residual. The user can use this to apply body forces or Neumann boundary conditions.

 None

form_compiler_options

dict | None

The options for the form compiler.

 None

jit_options

dict | None

The options for the JIT compiler.

 None
Source code in src/fenics_constitutive/solver/_solver.py 
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
class IncrSmallStrainProblem(NonlinearProblem):
    """
    A nonlinear problem for incremental small strain models. To be used with
    the dolfinx NewtonSolver.

    Args:
        laws: A list of tuples where the first element is the constitutive law and the second
            element is the cells for the submesh. If only one law is provided, it is assumed
            that the domain is homogenous. The cell indices should be local indices of the MPI-process.
        u: The displacement field. This is the unknown in the nonlinear problem.
        bcs: The Dirichlet boundary conditions.
        q_degree: The quadrature degree (Polynomial degree which the quadrature rule needs to integrate exactly).
        del_t: The maximal allowed time increment between steps. This is relevant if an explicit solver is used or if the 
            material model depends on the time, e.g. viscoelasticity. The actually used timestep may be overwritten by
            the solver if convergence issues occur or if the critical timestep is smaller.
        external_forces: The external forces applied to the system. This should be a list of ufl Forms which are subtracted from the residual. 
            The user can use this to apply body forces or Neumann boundary conditions.
        form_compiler_options: The options for the form compiler.
        jit_options: The options for the JIT compiler.
    """

    def __init__(
        self,
        laws: list[tuple[IncrSmallStrainModel, np.ndarray]] | IncrSmallStrainModel,
        u: df.fem.Function,
        bcs: list[df.fem.DirichletBC],
        q_degree: int,
        del_t: float = 1.0,
        external_forces: list[ufl.Form] | None = None,
        form_compiler_options: dict | None = None,
        jit_options: dict | None = None,
    ) -> None:
        mesh = u.function_space.mesh
        map_c = mesh.topology.index_map(mesh.topology.dim)
        num_cells = map_c.size_local + map_c.num_ghosts
        if isinstance(laws, IncrSmallStrainModel):
            laws = [(laws, np.arange(0, num_cells, dtype=np.int32))]

        constraint = laws[0][0].constraint
        assert all(
            law[0].constraint == constraint for law in laws
        ), "All laws must have the same constraint"

        element_spaces = ElementSpaces.create(mesh, constraint, q_degree)
        self.stress = IncrementalStress(element_spaces.stress_vector_space)
        self.tangent = fn_for(element_spaces.stress_tensor_space(mesh))

        self._law_on_submeshs: list[LawOnSubMesh] = []
        self.sim_time = SimulationTime(dt_max=del_t)

        self._law_on_submeshs = [
            create_law_on_submesh(law, local_cells, element_spaces)
            for law, local_cells in laws
        ]

        u_, du = ufl.TestFunction(u.function_space), ufl.TrialFunction(u.function_space)

        self.metadata = {"quadrature_degree": q_degree, "quadrature_scheme": "default"}
        self.dxm = ufl.dx(metadata=self.metadata)

        R_form = (
            ufl.inner(ufl_mandel_strain(u_, constraint), self.stress.current) * self.dxm
        )

        if external_forces is not None:
            R_form -= sum(external_forces)

        dR_form = (
            ufl.inner(
                ufl_mandel_strain(du, constraint),
                ufl.dot(self.tangent, ufl_mandel_strain(u_, constraint)),
            )
            * self.dxm
        )

        self.incr_disp = IncrementalDisplacement(u, q_degree)
        super().__init__(
                R_form,
                self.incr_disp.current,
                bcs=bcs,
                J=dR_form,
                form_compiler_options=form_compiler_options if form_compiler_options is not None else {},
                jit_options=jit_options if jit_options is not None else {}
            )

    @df.common.timed("constitutive-form-evaluation")
    def form(self, x: PETSc.Vec) -> None:
        """This function is called before the residual or Jacobian is
        computed. This is usually used to update ghost values, but here
        we use it to update the stress, tangent and history.

        Args:
            x: The vector containing the latest solution

        """
        super().form(x)

        self.incr_disp.update_current(x)
        self.form_without_petsc(evaluate_tangent=True)

    def form_without_petsc(self, evaluate_tangent: bool)-> None:
        """
        This function should be used when the solver does not require PETSc. We assume that the current solution
        is stored in `self.incr_displ`.

        Args:
            evaluate_tangent: Whether to evaluate the tangent. If `False`, the old values in `self.tangent` will be used. 
        """
        tangent = self.tangent if evaluate_tangent else None
        for law in self._law_on_submeshs:
            law.evaluate(self.sim_time, self.incr_disp, self.stress, tangent)

        self.stress.scatter_current()
        self.tangent.x.scatter_forward()

    def update(self) -> None:
        """
        Update the current displacement, stress and history.
        Any sensor evaluations that require you to know the current and 
        the previous state should be done before calling this function, 
        as it will update the previous state to the current state.
        """
        self.incr_disp.update_previous()
        self.stress.update_previous()

        for law in self._law_on_submeshs:
            law.update_history()

        self.sim_time.advance()

    # -------------------------------------------------------------------
    # NOTE: The following properties are used for backward compatibility
    # -------------------------------------------------------------------

    @property
    def _time(self) -> float:
        return self.sim_time.current

    @_time.setter
    def _time(self, value: float) -> None:
        self.sim_time.current = value

    @property
    def _del_t(self) -> float:
        return self.sim_time.dt

    @_del_t.setter
    def _del_t(self, value: float) -> None:
        self.sim_time.dt = value

    @property
    def _u(self) -> df.fem.Function:
        return self.incr_disp.current

    @property
    def _u0(self) -> df.fem.Function:
        return self.incr_disp.previous

    @property
    def stress_0(self) -> df.fem.Function:
        return self.stress.previous

    @property
    def stress_1(self) -> df.fem.Function:
        return self.stress.current

    @property
    def _history_0(self) -> list[dict[str, Function] | None]:
        """Return a list of history_0 dicts for all laws (for backward compatibility)."""

        def _history_or_none(law) -> dict[str, Function] | None:
            return law.history.history_0 if law.history else None

        return [_history_or_none(law) for law in self._law_on_submeshs]

    @property
    def _history_1(self) -> list[dict[str, Function] | None]:
        """Return a list of history_1 dicts for all laws (for backward compatibility)."""

        def _history_or_none(law) -> dict[str, Function] | None:
            return law.history.history_1 if law.history else None

        return [_history_or_none(law) for law in self._law_on_submeshs]

    @property
    def _del_grad_u(self) -> list[Function]:
        """Return a list of inc_disp_grad Functions for all laws (for backward compatibility)."""
        return [law.displacement_gradient_fn for law in self._law_on_submeshs]

form 

form(x: Vec) -> None

This function is called before the residual or Jacobian is computed. This is usually used to update ghost values, but here we use it to update the stress, tangent and history.

Parameters:

Name Type Description Default

x

Vec

The vector containing the latest solution

 required
Source code in src/fenics_constitutive/solver/_solver.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
@df.common.timed("constitutive-form-evaluation")
def form(self, x: PETSc.Vec) -> None:
    """This function is called before the residual or Jacobian is
    computed. This is usually used to update ghost values, but here
    we use it to update the stress, tangent and history.

    Args:
        x: The vector containing the latest solution

    """
    super().form(x)

    self.incr_disp.update_current(x)
    self.form_without_petsc(evaluate_tangent=True)

form_without_petsc 

form_without_petsc(evaluate_tangent: bool) -> None

This function should be used when the solver does not require PETSc. We assume that the current solution is stored in self.incr_displ.

Parameters:

Name Type Description Default

evaluate_tangent

bool

Whether to evaluate the tangent. If False, the old values in self.tangent will be used.

 required
Source code in src/fenics_constitutive/solver/_solver.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
def form_without_petsc(self, evaluate_tangent: bool)-> None:
    """
    This function should be used when the solver does not require PETSc. We assume that the current solution
    is stored in `self.incr_displ`.

    Args:
        evaluate_tangent: Whether to evaluate the tangent. If `False`, the old values in `self.tangent` will be used. 
    """
    tangent = self.tangent if evaluate_tangent else None
    for law in self._law_on_submeshs:
        law.evaluate(self.sim_time, self.incr_disp, self.stress, tangent)

    self.stress.scatter_current()
    self.tangent.x.scatter_forward()

update 

update() -> None

Update the current displacement, stress and history. Any sensor evaluations that require you to know the current and the previous state should be done before calling this function, as it will update the previous state to the current state.

Source code in src/fenics_constitutive/solver/_solver.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def update(self) -> None:
    """
    Update the current displacement, stress and history.
    Any sensor evaluations that require you to know the current and 
    the previous state should be done before calling this function, 
    as it will update the previous state to the current state.
    """
    self.incr_disp.update_previous()
    self.stress.update_previous()

    for law in self._law_on_submeshs:
        law.update_history()

    self.sim_time.advance()

critical_timestep 

critical_timestep(laws: list[tuple[IncrSmallStrainModel, ndarray]], density: list[float], u: Function, method: str = 'cdm', h: float | None = None) -> np.ndarray

Determines the critical timesteps for all submeshes. This assumes that the constitutive law returns a linear elastic tangent for . The input is not verified for consistency as this function is supposed to be called in the CDMSolver or any other solver.

Parameters:

Name Type Description Default

laws

list[tuple[IncrSmallStrainModel, ndarray]]

A list of tuples containing the constitutive law and the corresponding cells for each submesh.

 required

density

list[float]

A list of densities for each submesh.

 required

u

Function

The current displacement. This is used to determine the geometric dimension and the function space of the problem.

 required

method

str

The time integration method. This is used to determine the factor for the critical time step. Currently only "cdm" is implemented, which corresponds to the central difference method. The factor for the central difference method is 2, which means that the time step should be halved.

 'cdm'

h

float | None

The mesh size. If this is not provided, the mesh size will be determined based on the mesh and the cells for each submesh. This can be used to override the mesh size if the user already knows h for example from the mesh creation.

 None

Returns:

Type Description
ndarray

An array containing the critical timesteps of all submeshes.
Source code in src/fenics_constitutive/solver/central_difference_method.py 
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def critical_timestep(
    laws: list[tuple[IncrSmallStrainModel, np.ndarray]],
    density: list[float],
    u: df.fem.Function,
    method: str = "cdm",
    h: float | None = None,
) -> np.ndarray:
    r"""
    Determines the critical timesteps for all submeshes. This assumes that the constitutive law
    returns a linear elastic tangent for $\sigma=0,\varepsilon=0$. The input is not verified for
    consistency as this function is supposed to be called in the CDMSolver or any other solver.

    Args:
        laws: A list of tuples containing the constitutive law and the corresponding cells for each submesh.
        density: A list of densities for each submesh.
        u: The current displacement. This is used to determine the geometric dimension and the function space of the problem.
        method: The time integration method. This is used to determine the factor for the critical time step. Currently only 
            "cdm" is implemented, which corresponds to the central difference method. The factor for the central difference method is 2,
            which means that the time step should be halved.
        h: The mesh size. If this is not provided, the mesh size will be determined based on the mesh and the cells for each submesh.
            This can be used to override the mesh size if the user already knows `h` for example from the mesh creation.

    Returns:
        An array containing the critical timesteps of all submeshes.
    """
    method_to_factor = {"cdm": 2}
    factor = method_to_factor[method]
    mesh = u.function_space.mesh

    del_t: list[float] = []
    for (law, cells), density_ in zip(laws, density):
        tangent = np.zeros((law.stress_strain_dim, law.stress_strain_dim))
        stress = np.zeros(law.stress_strain_dim)
        grad_del_u = np.zeros((law.geometric_dim, law.geometric_dim))
        history = (
            {key: np.zeros(dim) for (key, dim) in law.history_dim.items()}
            if law.history_dim is not None
            else None
        )
        law.evaluate(
            0.0,
            1e-12,
            grad_del_u,
            stress,
            tangent.reshape(
                -1,
            ),
            history,
        )
        assert np.linalg.norm(tangent) > 0.0, (
            "The constitutive law must return a non-zero tangent"
        )
        h_ = mesh.h(mesh.topology.dim, cells) if h is None else np.array([h])
        h_min = h_.min()
        omega = _max_frequency_one_element(h_min, u, tangent, density_, law.constraint)
        del_t.append(factor / omega)
    return np.array(del_t)

diagonal_inverted_mass 

diagonal_inverted_mass(function_space: FunctionSpace, density: list[float], cells: list[ndarray]) -> df.fem.Function

Determine the inverse of the diagonal mass matrix. For intervals, quadrilaterals and hexahedra, the Gauss-Lobatto-Legendre quadrature is used to compute the mass matrix, which results in a diagonal mass matrix. For other cell types, this function is not implemented and an exception is raised.

Parameters:

Name Type Description Default

function_space

FunctionSpace

The function space for which the mass matrix is computed.

 required

density

list[float]

The density of the material. This is either a single value for a homogenous domain or a list of values for each submesh.

 required

cells

list[ndarray]

The cells corresponding to each submesh. This is used to assign the density to the correct cells in the mass matrix assembly.

 required

Returns:

Type Description
Function

The inverted diagonal mass matrix as a function object.
Source code in src/fenics_constitutive/solver/central_difference_method.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
def diagonal_inverted_mass(
    function_space: df.fem.FunctionSpace, density: list[float], cells: list[np.ndarray]
) -> df.fem.Function:
    """
    Determine the inverse of the diagonal mass matrix. For intervals, quadrilaterals and hexahedra, the
    Gauss-Lobatto-Legendre quadrature is used to compute the mass matrix, which results in a diagonal mass matrix.
    For other cell types, this function is not implemented and an exception is raised.

    Args:
        function_space: The function space for which the mass matrix is computed.
        density: The density of the material. This is either a single value for a homogenous domain or a list of values for each submesh.
        cells: The cells corresponding to each submesh. This is used to assign the density to the correct cells in the mass matrix assembly.

    Returns:
        The inverted diagonal mass matrix as a function object.
    """
    mesh_cell = function_space.mesh.ufl_cell().cellname()
    basix_cell = basix.CellType[mesh_cell]
    if basix_cell in [
        basix.CellType.interval,
        basix.CellType.quadrilateral,
        basix.CellType.hexahedron,
    ]:
        # do gll integration
        # todo:adapt for higher order elements
        p_degree_to_q_degree = {1: 1, 2: 2}
        geo_dim = function_space.mesh.geometry.dim
        V_degree = function_space.ufl_element().degree
        action_fn = cast(df.fem.Function, df.fem.Function(function_space))
        action_fn.x.array[:] = 1.0
        action_fn.x.scatter_forward()
        q_degree = p_degree_to_q_degree[V_degree]

        metadata = {"quadrature_degree": q_degree, "quadrature_scheme": "gll"}
        dxm = ufl.dx(metadata=metadata)
        if len(density) > 1:
            density_space = df.fem.functionspace(function_space.mesh, ("DG", 0))
            density_fn = cast(df.fem.Function, df.fem.Function(density_space))
            for density_, cells_ in zip(density, cells):
                density_fn.x.array[cells_] = density_
            density_fn.x.scatter_forward()
        else:
            density_fn = df.fem.Constant(function_space.mesh, density[0])

        u_ = ufl.TestFunction(function_space)
        v_ = ufl.TrialFunction(function_space)
        mass_form = cast(
            ufl.Form, ufl.action(density_fn * ufl.inner(u_, v_) * dxm, action_fn)
        )
        M_action = cast(df.fem.Function, df.fem.Function(function_space))
        df.fem.assemble_vector(M_action.x.array, df.fem.form(mass_form))
        M_action.x.scatter_reverse(df.la.InsertMode.add)
    else:
        raise Exception(
            "Only implemented for intervals, quadrilaterals and hexahedral elements"
        )
    M_action.x.array[:] = 1.0 / M_action.x.array[:]
    M_action.x.scatter_forward()
    return M_action

ufl_mandel_strain 

ufl_mandel_strain(u: Expr, constraint: StressStrainConstraint) -> ufl.core.expr.Expr

Compute the strain in Mandel notation from the displacement field.

Parameters:

Name Type Description Default

u

Expr

Displacement field.

 required

constraint

StressStrainConstraint

Constraint that the model is implemented for.

 required

Returns:

Type Description
Expr

Vector-valued UFL expression of the strain in Mandel notation.
Source code in src/fenics_constitutive/solver/utils.py 
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
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
def ufl_mandel_strain(
    u: ufl.core.expr.Expr, constraint: StressStrainConstraint
) -> ufl.core.expr.Expr:
    """
    Compute the strain in Mandel notation from the displacement field.

    Args:
        u: Displacement field.
        constraint: Constraint that the model is implemented for.

    Returns:
        Vector-valued UFL expression of the strain in Mandel notation.
    """
    shape = len(u.ufl_shape)
    geometric_dim = u.ufl_shape[0] if shape > 0 else 1
    assert geometric_dim == constraint.geometric_dim
    match constraint:
        case StressStrainConstraint.UNIAXIAL_STRAIN:
            return ufl.nabla_grad(u)
        case StressStrainConstraint.UNIAXIAL_STRESS:
            return ufl.nabla_grad(u)
        case StressStrainConstraint.PLANE_STRAIN:
            grad_u = ufl.nabla_grad(u)
            return ufl.as_vector(
                [
                    grad_u[0, 0],
                    grad_u[1, 1],
                    0.0,
                    1 / 2**0.5 * (grad_u[0, 1] + grad_u[1, 0]),
                ]
            )
        case StressStrainConstraint.PLANE_STRESS:
            grad_u = ufl.nabla_grad(u)
            return ufl.as_vector(
                [
                    grad_u[0, 0],
                    grad_u[1, 1],
                    0.0,
                    1 / 2**0.5 * (grad_u[0, 1] + grad_u[1, 0]),
                ]
            )
        case StressStrainConstraint.FULL:
            grad_u = ufl.nabla_grad(u)
            return ufl.as_vector(
                [
                    grad_u[0, 0],
                    grad_u[1, 1],
                    grad_u[2, 2],
                    1 / 2**0.5 * (grad_u[0, 1] + grad_u[1, 0]),
                    1 / 2**0.5 * (grad_u[0, 2] + grad_u[2, 0]),
                    1 / 2**0.5 * (grad_u[1, 2] + grad_u[2, 1]),
                ]
            )