First test for solver and linear elasticity

examples/linear_elasticity/linear_elasticity_model.py

@@ -69,15 +69,15 @@ def evaluate(
         history: np.ndarray | dict[str, np.ndarray] | None,
     ) -> None:
         assert (
-            grad_del_u.size / (self.geometric_dim**2)
-            == mandel_stress.size / self.stress_strain_dim
-            == tangent.size / (self.stress_strain_dim**2)
+            grad_del_u.size // (self.geometric_dim**2)
+            == mandel_stress.size // self.stress_strain_dim
+            == tangent.size // (self.stress_strain_dim**2)
         )
-        n_gauss = grad_del_u.shape / (self.geometric_dim**2)
-
+        n_gauss = grad_del_u.size // (self.geometric_dim**2)
+        mandel_view = mandel_stress.reshape(-1, self.stress_strain_dim)
         strain_increment = strain_from_grad_u(grad_del_u, self.constraint)
-        mandel_stress += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D
-        tangent[:] = np.tile(self.D, n_gauss)
+        mandel_view += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D
+        tangent[:] = np.tile(self.D.flatten(), n_gauss)
 
     @property
     def constraint(self) -> Constraint:

examples/linear_elasticity/test_elasticity.py

@@ -0,0 +1,49 @@
+from mpi4py import MPI
+import dolfinx as df
+from linear_elasticity_model import LinearElasticityModel
+from fenics_constitutive import Constraint, IncrSmallStrainProblem
+import numpy as np
+from dolfinx.nls.petsc import NewtonSolver
+
+youngs_modulus = 42.
+poissons_ratio = 0.3
+
+def test_uniaxial_stress():
+    mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
+    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    u = df.fem.Function(V)
+    law = LinearElasticityModel(
+        parameters={"E": youngs_modulus, "nu": poissons_ratio},
+        constraint=Constraint.UNIAXIAL_STRESS,
+    )
+
+    def left_boundary(x):
+        return np.isclose(x[0], 0.0)
+
+    def right_boundary(x):
+        return np.isclose(x[0], 1.0)
+
+    displacement = df.fem.Constant(mesh, 0.01)
+    dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary)
+    dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary)
+    bc_left = df.fem.dirichletbc(df.fem.Constant(mesh,0.0), dofs_left, V)
+    bc_right = df.fem.dirichletbc(displacement, dofs_right, V)
+
+
+    problem = IncrSmallStrainProblem(
+        law,
+        u,
+        [bc_left, bc_right],
+    )
+
+    solver = NewtonSolver(MPI.COMM_WORLD, problem)
+    n, converged = solver.solve(u)
+    assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01)
+
+    problem.update()
+    assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01)
+    assert np.max(problem._u0.x.array) == displacement.value
+
+    displacement.value = 0.02
+    n, converged = solver.solve(u)
+    assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / (youngs_modulus * 0.02)

src/fenics_constitutive/solver.py

@@ -65,7 +65,7 @@ def build_history(
 class IncrSmallStrainProblem(df.fem.petsc.NonlinearProblem):
     def __init__(
         self,
-        laws: list[tuple[IncrSmallStrainModel, np.ndarray]],
+        laws: list[tuple[IncrSmallStrainModel, np.ndarray]] | IncrSmallStrainModel,
         u: df.fem.Function,
         bcs: list[df.fem.DirichletBCMetaClass],
         q_degree: int = 1,
@@ -76,6 +76,8 @@ def __init__(
         map_c = mesh.topology.index_map(mesh.topology.dim)
         num_cells = map_c.size_local + map_c.num_ghosts
         cells = np.arange(0, num_cells, dtype=np.int32)
+        if isinstance(laws, IncrSmallStrainModel):
+            laws = [(laws, cells)]
 
         constraint = laws[0][0].constraint
         assert all(
@@ -99,13 +101,20 @@ def __init__(
             quad_scheme="default",
             shape=(constraint.stress_strain_dim(), constraint.stress_strain_dim()),
         )
+        Q_grad_u_e = ufl.TensorElement(
+            "Quadrature",
+            mesh.ufl_cell(),
+            q_degree,
+            quad_scheme="default",
+            shape=(gdim, gdim),
+        )
         QV = df.fem.FunctionSpace(mesh, QVe)
         QT = df.fem.FunctionSpace(mesh, QTe)
 
         self.laws: list[tuple[IncrSmallStrainModel, np.ndarray]] = []
         self.submesh_maps: list[SubSpaceMap] = []
 
-        self._del_strain = []
+        self._del_grad_u = []
         self._stress = []
         self._history_0 = []
         self._history_1 = []
@@ -123,8 +132,11 @@ def __init__(
                         cells, QV, return_subspace=True
                     )
                     self.submesh_maps.append(subspace_map)
-                    self._del_strain.append(df.fem.Function(QV_subspace))
                     self._stress.append(df.fem.Function(QV_subspace))
+
+                    #subspace for grad u
+                    Q_grad_u_subspace = df.fem.FunctionSpace(submesh, Q_grad_u_e)
+                    self._del_grad_u.append(df.fem.Function(Q_grad_u_subspace))
 
                     #subspace for tanget
                     QT_subspace = df.fem.FunctionSpace(submesh, QTe)
@@ -143,7 +155,9 @@ def __init__(
                 law, cells = laws[0]
                 self.laws.append((law, cells))
 
-                self._del_strain.append(df.fem.Function(QV))
+                #subspace for grad u
+                Q_grad_u_space = df.fem.FunctionSpace(mesh, Q_grad_u_e)
+                self._del_grad_u.append(df.fem.Function(Q_grad_u_space))
 
                 #Spaces for history
                 history_0 = build_history(law, mesh, q_degree)
@@ -184,8 +198,8 @@ def __init__(
         basix_celltype = getattr(basix.CellType, mesh.topology.cell_type.name)
         self.q_points, _ = basix.make_quadrature(basix_celltype, q_degree)
 
-        self.del_strain_expr = df.fem.Expression(
-            ufl_mandel_strain(self._u - self._u0, constraint), self.q_points
+        self.del_grad_u_expr = df.fem.Expression(
+            ufl.nabla_grad(self._u - self._u0), self.q_points
         )
 
         # ### Note on JIT compilation of UFL forms
@@ -225,20 +239,20 @@ def form(self, x: PETSc.Vec):
         super().form(x)
 
         assert (
-            x is self._u.vector
+            x.array.data == self._u.vector.array.data
         ), "The solution vector must be the same as the one passed to the MechanicsProblem"
         if len(self.laws) > 1:
             for k, (law, cells) in enumerate(self.laws):
                 with df.common.Timer("strain_evaluation"):
                     #TODO: test this!!
-                    self.del_strain_expr.eval(cells, self._del_strain[k].x.array.reshape(cells.size, -1))
+                    self.del_grad_u_expr.eval(cells, self._del_grad_u[k].x.array.reshape(cells.size, -1))
 
                 with df.common.Timer("stress_evaluation"):
                     self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k])
                     self._history_1[k].x.array[:] = self._history_0[k].x.array
                     law.evaluate(
                         self._time,
-                        self._del_strain[k].x.array,
+                        self._del_grad_u[k].x.array,
                         self._stress[k].x.array,
                         self._tangent[k].x.array,
                         self._history_1[k].x.array,  # history,
@@ -250,17 +264,18 @@ def form(self, x: PETSc.Vec):
         else:
             law, cells = self.laws[0]
             with df.common.Timer("strain_evaluation"):
-                self.del_strain_expr.eval(cells, self._del_strain[0].x.array)
+                self.del_grad_u_expr.eval(cells, self._del_grad_u[0].x.array.reshape(cells.size, -1))
 
             with df.common.Timer("stress_evaluation"):
                 self.stress_1.x.array[:] = self.stress_0.x.array
-                self._history_1[0].x.array[:] = self._history_0[0].x.array
+                if law.history_dim is not None:
+                    self._history_1[0].x.array[:] = self._history_0[0].x.array
                 law.evaluate(
                     self._time,
-                    self._del_strain[0].x.array,
+                    self._del_grad_u[0].x.array,
                     self.stress_1.x.array,
                     self.tangent.x.array,
-                    self._history_1[0].x.array,  # history,
+                    self._history_1[0].x.array if law.history_dim is not None else None,  # history,
                 )
 
         self.stress_1.x.scatter_forward()

src/fenics_constitutive/stress_strain.py

@@ -20,7 +20,10 @@ def ufl_mandel_strain(
         Vector-valued UFL expression of the mandel strain.
     """
     strain_dim = constraint.stress_strain_dim()
-    assert u.ufl_shape == (constraint.geometric_dim(),)
+
+    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 Constraint.UNIAXIAL_STRAIN:
             return ufl.nabla_grad(u)

tests/test_maps.py

@@ -3,7 +3,7 @@
 import ufl
 from mpi4py import MPI
 
-from fenics_constitutive import SubSpaceMap, _build_view, build_subspace_map
+from fenics_constitutive import SubSpaceMap, build_subspace_map
 
 
 def test_subspace_vector_map_vector_equals_tensor_map():
@@ -26,8 +26,8 @@ def test_subspace_vector_map_vector_equals_tensor_map():
     QT_space = df.fem.FunctionSpace(mesh, QT)
 
     Q_map, _ = build_subspace_map(cells, Q_space)
-    QV_map, _ = _build_view(cells, QV_space)
-    QT_map, _ = _build_view(cells, QT_space)
+    QV_map, _ = build_subspace_map(cells, QV_space)
+    QT_map, _ = build_subspace_map(cells, QT_space)
 
     assert np.all(Q_map.parent == QV_map.parent)
     assert np.all(Q_map.child == QV_map.child)