3d test case

examples/linear_elasticity/test_elasticity.py

@@ -3,6 +3,7 @@
 import dolfinx as df
 import numpy as np
 import pytest
+import ufl
 from dolfinx.nls.petsc import NewtonSolver
 from linear_elasticity_model import LinearElasticityModel
 from mpi4py import MPI
@@ -162,3 +163,61 @@ def right_boundary(x):
     assert abs(problem.stress_0.x.array[0] - analytical_stress) < 1e-10 / (
         analytical_stress
     )
+
+
+def test_3d():
+    # test the 3d case against a pure fenics implementation
+    mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 2, 2, 2)
+    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    u = df.fem.Function(V)
+    law = LinearElasticityModel(
+        parameters={"E": youngs_modulus, "nu": poissons_ratio},
+        constraint=Constraint.FULL,
+    )
+
+    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_2d, np.ndarray([0.01,0.0]))
+    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(np.array([0.0,0.0, 0.0]), dofs_left, V)
+    bc_right = df.fem.dirichletbc(np.array([0.01, 0.0, 0.0]), dofs_right, V)
+
+    problem = IncrSmallStrainProblem(
+        law,
+        u,
+        [bc_left, bc_right],
+    )
+
+    solver = NewtonSolver(MPI.COMM_WORLD, problem)
+    n, converged = solver.solve(u)
+    problem.update()
+
+    v_, u_ = ufl.TestFunction(V), ufl.TrialFunction(V)
+
+    def eps(v):
+        return ufl.sym(ufl.grad(v))
+
+    def sigma(v):
+        eps = ufl.sym(ufl.grad(v))
+        lam, mu = (
+            youngs_modulus
+            * poissons_ratio
+            / ((1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio)),
+            youngs_modulus / (2.0 * (1.0 + poissons_ratio)),
+        )
+        return 2.0 * mu * eps + lam * ufl.tr(eps) * ufl.Identity(len(v))
+    zero = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0]))
+    a = ufl.inner(ufl.grad(v_), sigma(u_)) * ufl.dx
+    L = ufl.dot(zero ,v_) * ufl.dx
+
+    u_fenics = u.copy()
+    problem_fenics = df.fem.petsc.LinearProblem(a, L, [bc_left, bc_right], u=u_fenics)
+    problem_fenics.solve()
+
+    # Check that the solution is the same
+    assert np.linalg.norm(u_fenics.x.array - u.x.array) < 1e-8/np.linalg.norm(u_fenics.x.array)

src/fenics_constitutive/solver.py

@@ -15,12 +15,12 @@ def build_history(
     law: IncrSmallStrainModel, mesh: df.mesh.Mesh, q_degree: int
 ) -> df.fem.Function | dict[str, df.fem.Function] | None:
     """Build the history space and function(s) for the given law.
-    
+
     Args:
         law: The constitutive law.
         mesh: Either the full mesh for a homogenous domain or the submesh.
         q_degree: The quadrature degree.
-    
+
     Returns:
         The history function(s) for the given law.
 
@@ -85,7 +85,9 @@ def __init__(
         ), "All laws must have the same constraint"
 
         gdim = mesh.ufl_cell().geometric_dimension()
-        assert constraint.geometric_dim() == gdim, "Geometric dimension mismatch between mesh and laws"
+        assert (
+            constraint.geometric_dim() == gdim
+        ), "Geometric dimension mismatch between mesh and laws"
 
         QVe = ufl.VectorElement(
             "Quadrature",
@@ -133,16 +135,16 @@ def __init__(
                     )
                     self.submesh_maps.append(subspace_map)
                     self._stress.append(df.fem.Function(QV_subspace))
-                    
-                    #subspace for grad u
+
+                    # 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
+                    # subspace for tanget
                     QT_subspace = df.fem.FunctionSpace(submesh, QTe)
                     self._tangent.append(df.fem.Function(QT_subspace))
 
-                    #subspaces for history
+                    # subspaces for history
                     history_0 = build_history(law, submesh, q_degree)
                     history_1 = (
                         {key: fn.copy() for key, fn in history_0.items()}
@@ -154,12 +156,12 @@ def __init__(
             else:
                 law, cells = laws[0]
                 self.laws.append((law, cells))
-                
-                #subspace for grad u
+
+                # 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
+                # Spaces for history
                 history_0 = build_history(law, mesh, q_degree)
                 history_1 = (
                     {key: fn.copy() for key, fn in history_0.items()}
@@ -222,7 +224,9 @@ def a(self) -> df.fem.FormMetaClass:
                 self._u,
                 self._bcs,
                 self.dR_form,
-                form_compiler_options=self._form_compiler_options if self._form_compiler_options is not None else {},
+                form_compiler_options=self._form_compiler_options
+                if self._form_compiler_options is not None
+                else {},
                 jit_options=self._jit_options if self._jit_options is not None else {},
             )
 
@@ -244,8 +248,10 @@ def form(self, x: PETSc.Vec):
         if len(self.laws) > 1:
             for k, (law, cells) in enumerate(self.laws):
                 with df.common.Timer("strain_evaluation"):
-                    #TODO: test this!!
-                    self.del_grad_u_expr.eval(cells, self._del_grad_u[k].x.array.reshape(cells.size, -1))
+                    # TODO: test this!!
+                    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])
@@ -256,7 +262,9 @@ def form(self, x: PETSc.Vec):
                         self._del_grad_u[k].x.array,
                         self._stress[k].x.array,
                         self._tangent[k].x.array,
-                        self._history_1[k].x.array if law.history_dim is not None else None,
+                        self._history_1[k].x.array
+                        if law.history_dim is not None
+                        else None,
                     )
 
                 with df.common.Timer("stress-local-to-global"):
@@ -265,7 +273,9 @@ def form(self, x: PETSc.Vec):
         else:
             law, cells = self.laws[0]
             with df.common.Timer("strain_evaluation"):
-                self.del_grad_u_expr.eval(cells, self._del_grad_u[0].x.array.reshape(cells.size, -1))
+                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
@@ -276,7 +286,9 @@ def form(self, x: PETSc.Vec):
                     self._del_grad_u[0].x.array,
                     self.stress_1.x.array,
                     self.tangent.x.array,
-                    self._history_1[0].x.array if law.history_dim is not None else None,  # history,
+                    self._history_1[0].x.array
+                    if law.history_dim is not None
+                    else None,  # history,
                 )
 
         self.stress_1.x.scatter_forward()
@@ -290,7 +302,7 @@ def update(self):
 
         self.stress_0.x.array[:] = self.stress_1.x.array
         self.stress_0.x.scatter_forward()
-        
+
         for k, (law, _) in enumerate(self.laws):
             match law.history_dim:
                 case int():
@@ -300,5 +312,7 @@ def update(self):
                     pass
                 case dict():
                     for key in law.history_dim:
-                        self._history_0[k][key].x.array[:] = self._history_1[k][key].x.array
+                        self._history_0[k][key].x.array[:] = self._history_1[k][
+                            key
+                        ].x.array
                         self._history_0[k][key].x.scatter_forward()