Dolfinx 0.9.0

environment.yml

@@ -1,10 +1,9 @@
 name: fenics-constitutive
 channels:
   - conda-forge
-  - defaults
 dependencies:
   - python>=3.10
-  - fenics-dolfinx=0.6.0
+  - fenics-dolfinx>=0.9.0
   - numpy
   - pip
   - pytest

examples/linear_elasticity/test_elasticity.py

@@ -13,15 +13,16 @@
     IncrSmallStrainProblem,
     PlaneStrainFrom3D,
     UniaxialStrainFrom3D,
+    norm,
 )
 
 youngs_modulus = 42.0
 poissons_ratio = 0.3
 
 
 def test_uniaxial_stress():
-    mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 4)
-    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
+    V = df.fem.functionspace(mesh, ("CG", 1))
     u = df.fem.Function(V)
     law = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -51,25 +52,25 @@ def right_boundary(x):
     n, converged = solver.solve(u)
 
     # Compare the result with the analytical solution
-    assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (
-        youngs_modulus * 0.01
-    )
+    diff = problem.stress_1 - ufl.as_vector([youngs_modulus * 0.01,])
+    assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.01
 
     problem.update()
     # Check that the stress is updated correctly
-    assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (
-        youngs_modulus * 0.01
-    )
+    diff = problem.stress_0 - ufl.as_vector([youngs_modulus * 0.01,])
+    assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.01
+
     # Check that the displacement is updated correctly
-    assert np.max(problem._u0.x.array) == displacement.value
+    max_u = MPI.COMM_WORLD.allreduce(np.max(problem._u0.x.array), MPI.MAX)
+    assert max_u == displacement.value
 
     displacement.value = 0.02
     n, converged = solver.solve(u)
 
     # Compare the result of the updated problem with new BC with the analytical solution
-    assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / (
-        youngs_modulus * 0.02
-    )
+    diff = problem.stress_1 - ufl.as_vector([youngs_modulus * 0.02,])
+    assert norm(diff, problem.dxm) < 1e-10 / youngs_modulus * 0.02
+
 
 
 @pytest.mark.parametrize(
@@ -83,22 +84,23 @@ def right_boundary(x):
 )
 def test_uniaxial_stress_two_laws(factor: float):
     mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
-    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1))
     u = df.fem.Function(V)
+    cells_local = mesh.topology.index_map(mesh.topology.dim).global_to_local(np.arange(2, dtype=np.int32))
     laws = [
         (
             LinearElasticityModel(
                 parameters={"E": youngs_modulus, "nu": poissons_ratio},
                 constraint=Constraint.UNIAXIAL_STRESS,
             ),
-            np.array([0], dtype=np.int32),
+            cells_local[0:1],
         ),
         (
             LinearElasticityModel(
                 parameters={"E": factor * youngs_modulus, "nu": poissons_ratio},
                 constraint=Constraint.UNIAXIAL_STRESS,
             ),
-            np.array([1], dtype=np.int32),
+            cells_local[1:2],
         ),
     ]
 
@@ -138,7 +140,7 @@ def right_boundary(x):
 
 def test_uniaxial_strain():
     mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
-    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1))
     u = df.fem.Function(V)
     law = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -173,10 +175,8 @@ def right_boundary(x):
         * (1.0 - poissons_ratio)
         / ((1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio))
     ) * displacement.value
-
-    assert abs(problem.stress_0.x.array[0] - analytical_stress) < 1e-10 / (
-        analytical_stress
-    )
+    diff = problem.stress_0 - ufl.as_vector([analytical_stress,])
+    assert norm(diff, problem.dxm) < 1e-10 / analytical_stress
 
     # test the converter from 3D model to uniaxial strain model
     law_3d = LinearElasticityModel(
@@ -196,32 +196,28 @@ def right_boundary(x):
     problem_3d.update()
 
     # test that sigma_11 is the same as the analytical solution
-    assert abs(problem_3d.stress_0.x.array[0] - analytical_stress) < 1e-10 / (
-        analytical_stress
-    )
+    diff = problem_3d.stress_0 - ufl.as_vector([analytical_stress,])
+    assert norm(diff, problem_3d.dxm) < 1e-10 / analytical_stress
+
     # test that the stresses of the problem with uniaxial strain model
     # are the same as the stresses of the problem with the converted 3D model
-    assert (
-        np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array)
-        / np.linalg.norm(problem.stress_0.x.array)
-        < 1e-14
-    )
+    diff =  problem_3d.stress_0 - problem.stress_0
+    assert norm(diff, problem_3d.dxm) < 1e-10 / norm(problem.stress_0, problem.dxm)
 
     # test that the shear stresses are zero. Since this is uniaxial strain, the
     # stress can have nonzero \sigma_22 and \sigma_33 components
     assert np.linalg.norm(wrapped_1d_law.stress_3d[3:6]) < 1e-14
     # test that the displacement is the same in both versions
+    diff = problem_3d._u - problem._u
     assert (
-        np.linalg.norm(problem_3d._u.x.array - problem._u.x.array)
-        / np.linalg.norm(problem._u.x.array)
-        < 1e-14
+        norm(diff, problem_3d.dxm) < 1e-14 / norm(problem._u, problem.dxm)
     )
 
 
 def test_plane_strain():
     # sanity check if out of plane stress is NOT zero
     mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1,(2,)))
     u = df.fem.Function(V)
     law = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -249,14 +245,9 @@ def right_boundary(x):
     solver = NewtonSolver(MPI.COMM_WORLD, problem)
     n, converged = solver.solve(u)
     problem.update()
-    assert (
-        np.linalg.norm(
-            problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[
-                :, 2
-            ]
-        )
-        > 1e-2
-    )
+    # test that the stress is nonzero in 33 direction
+    assert norm(problem.stress_0[2], problem.dxm) > 1e-2
+
     # test the model conversion from 3D to plane strain
     law_3d = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -274,32 +265,22 @@ def right_boundary(x):
     n, converged = solver_3d.solve(u_3d)
     problem_3d.update()
     # test that the stress is nonzero in 33 direction
-    assert (
-        np.linalg.norm(
-            problem_3d.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[
-                :, 2
-            ]
-        )
-        > 1e-2
-    )
+    assert norm(problem_3d.stress_0[2], problem.dxm) > 1e-2
+
     # test that the displacement is the same in both versions
+    diff = problem_3d._u - problem._u
     assert (
-        np.linalg.norm(problem_3d._u.x.array - problem._u.x.array)
-        / np.linalg.norm(problem._u.x.array)
-        < 1e-14
+        norm(diff, problem.dxm)/ norm(problem._u, problem.dxm) < 1e-14 
     )
     # test that the stresses are the same in both versions
-    assert (
-        np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array)
-        / np.linalg.norm(problem.stress_0.x.array)
-        < 1e-14
-    )
+    diff = problem_3d.stress_0 - problem.stress_0
+    assert norm(diff, problem.dxm) / norm(problem.stress_0, problem.dxm)< 1e-10 
 
 
 def test_plane_stress():
     # just a sanity check if out of plane stress is really zero
     mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1,(2,)))
     u = df.fem.Function(V)
     law = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -328,20 +309,15 @@ def right_boundary(x):
     solver = NewtonSolver(MPI.COMM_WORLD, problem)
     n, converged = solver.solve(u)
     problem.update()
-    assert (
-        np.linalg.norm(
-            problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[
-                :, 2
-            ]
-        )
-        < 1e-10
-    )
+    # test that the out of plane stress is zero
+    assert norm(problem.stress_0[2], problem.dxm) < 1e-10
+
 
 
 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))
+    V = df.fem.functionspace(mesh, ("CG", 1,(3,)))
     u = df.fem.Function(V)
     law = LinearElasticityModel(
         parameters={"E": youngs_modulus, "nu": poissons_ratio},
@@ -395,10 +371,13 @@ def sigma(v):
     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
-    )
+    diff = u_fenics - u
+    assert norm(diff, problem.dxm) < 1e-8 / norm(u, problem.dxm)
 
 
 if __name__ == "__main__":
     test_uniaxial_stress()
+    test_uniaxial_strain()
+    test_plane_stress()
+    test_plane_strain()
+    test_3d()

examples/mises_plasticity/test_plasticity.py

@@ -10,8 +10,9 @@
 
 
 def test_uniaxial_strain_3d():
+    #no MPI, one cell only
     mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1)
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1, (3,)))
     u = df.fem.Function(V)
     matparam = {
         "p_ka": 175000,
@@ -49,8 +50,8 @@ def right(x):
     )
 
     dirichlet = [fix_ux_left, move_ux_right]
-    #
-    problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=2)
+
+    problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=1)
 
     solver = NewtonSolver(MPI.COMM_WORLD, problem)
 
@@ -62,15 +63,15 @@ def right(x):
     load = [0.0]
 
     for inc, time in enumerate(load_steps):
-        print("Load Increment:", inc)
+        #print("Load Increment:", inc)
 
         current_disp = time * max_disp
         scalar_x.value = current_disp
 
         niter, converged = solver.solve(u)
         problem.update()
 
-        print(f"Converged: {converged} in {niter} iterations.")
+        #print(f"Converged: {converged} in {niter} iterations.")
         iterations = np.append(iterations, niter)
 
         stress_values = []
@@ -106,7 +107,7 @@ def right(x):
 
 def test_uniaxial_cyclic_strain_3d():
     mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 1, 1, 1)
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1, (3,)))
     u = df.fem.Function(V)
     matparam = {
         "p_ka": 175000,
@@ -145,7 +146,7 @@ def right(x):
 
     dirichlet = [fix_ux_left, move_ux_right]
     #
-    problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=2)
+    problem = IncrSmallStrainProblem(law, u, dirichlet, q_degree=1)
 
     solver = NewtonSolver(MPI.COMM_WORLD, problem)
 
@@ -157,15 +158,15 @@ def right(x):
     load = [0.0]
 
     for inc, time in enumerate(load_steps):
-        print("Load Increment:", inc)
+        #print("Load Increment:", inc)
 
         current_disp = np.sin(time) * max_disp
         scalar_x.value = current_disp
 
         niter, converged = solver.solve(u)
         problem.update()
 
-        print(f"Converged: {converged} in {niter} iterations.")
+        #print(f"Converged: {converged} in {niter} iterations.")
         iterations = np.append(iterations, niter)
 
         stress_values = []

examples/visco_elasticity/test_viscoelasticity.py

@@ -24,7 +24,7 @@ def test_relaxation_uniaxial_stress(mat: IncrSmallStrainModel):
     """stress relaxation under uniaxial tension test for 1D, displacement controlled"""
 
     mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
-    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1))
     u = df.fem.Function(V)
     law = mat(
         parameters={"E0": youngs_modulus, "E1": visco_modulus, "tau": relaxation_time},
@@ -89,8 +89,8 @@ def right_boundary(x):
         strain.append(problem._history_1[0]["strain"].x.array[-1])
         viscostrain.append(problem._history_1[0]["strain_visco"].x.array[-1])
 
-    print("0", time[0], disp[0], stress[0], strain[0], viscostrain[0])
-    print("end", time[-1], disp[-1], stress[-1], strain[-1], viscostrain[-1])
+    #print("0", time[0], disp[0], stress[0], strain[0], viscostrain[0])
+    #print("end", time[-1], disp[-1], stress[-1], strain[-1], viscostrain[-1])
     # analytic solution
     if isinstance(law, SpringKelvinModel):
         # analytic solution for 1D Kelvin model
@@ -158,7 +158,7 @@ def z_boundary(x):
     else:
         raise ValueError(f"Dimension {dim} not supported")
 
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1, (dim,)))
     u = df.fem.Function(V)
 
     def left_boundary(x):
@@ -288,7 +288,7 @@ def test_kelvin_vs_maxwell():
     """1D test with uniaxial stress, compare Kelvin and Maxwell model."""
 
     mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2)
-    V = df.fem.FunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1))
     u = df.fem.Function(V)
 
     # Kelvin material
@@ -402,7 +402,7 @@ def z_boundary(x):
     else:
         raise ValueError(f"Dimension {dim} not supported")
 
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1,(dim,)))
     u = df.fem.Function(V)
 
     def left_boundary(x):
@@ -525,7 +525,7 @@ def y_boundary(x):
 
 def create_meshtags(
     domain: df.mesh.Mesh, entity_dim: int, markers: dict[str, tuple[int, Callable]]
-) -> tuple[df.mesh.MeshTagsMetaClass, dict[str, int]]:
+) -> tuple[df.mesh.MeshTags, dict[str, int]]:
     """Creates meshtags for the given markers.
     This code is part of the FEniCSx tutorial
     by Jørgen S. Dokken.
@@ -605,7 +605,7 @@ def z1_boundary(x):
         )
         fixed_vector = np.array([0.0, 0.0, 0.0])
 
-    V = df.fem.VectorFunctionSpace(mesh, ("CG", 1))
+    V = df.fem.functionspace(mesh, ("CG", 1,(dim,)))
     u = df.fem.Function(V)
 
     # boundaries