Merge pull request #88 from finsberg/remove-utils

Remove utils.py in test folder and use scifem more

Author: finsberg

src/pulse/kinematics.py

@@ -81,6 +81,10 @@ def IsochoricDeformationGradient(F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
     return pow(J, -1.0 / float(dim)) * F
 
 
+def IsochoricDeformationGradient_from_u(u) -> ufl.core.expr.Expr:
+    return IsochoricDeformationGradient(DeformationGradient(u))
+
+
 def Jacobian(F: ufl.core.expr.Expr) -> ufl.core.expr.Expr:
     r"""Determinant of the deformation gradient
 

src/pulse/utils.py

@@ -1,52 +1,8 @@
 import dolfinx
 import numpy as np
 import numpy.typing as npt
-
-
-def vertex_to_dofmap(V: dolfinx.fem.FunctionSpace) -> npt.NDArray[np.int32]:
-    """Compute a mapping from vertices to dofs in a function space
-
-    Parameters
-    ----------
-    V : dolfinx.fem.FunctionSpace
-        The function space
-
-    Returns
-    -------
-    npt.NDArray[np.int32]
-        The mapping from vertices to dofs
-    """
-    mesh = V.mesh
-    num_vertices_per_cell = dolfinx.cpp.mesh.cell_num_entities(mesh.topology.cell_type, 0)
-
-    dof_layout = np.empty((num_vertices_per_cell,), dtype=np.int32)
-    for i in range(num_vertices_per_cell):
-        var = V.dofmap.dof_layout.entity_dofs(0, i)
-        assert len(var) == 1
-        dof_layout[i] = var[0]
-
-    num_vertices = mesh.topology.index_map(0).size_local + mesh.topology.index_map(0).num_ghosts
-
-    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
-    assert (c_to_v.offsets[1:] - c_to_v.offsets[:-1] == c_to_v.offsets[1]).all(), (
-        "Single cell type supported"
-    )
-
-    vertex_to_dof_map = np.empty(num_vertices, dtype=np.int32)
-    vertex_to_dof_map[c_to_v.array] = V.dofmap.list[:, dof_layout].reshape(-1)
-    return vertex_to_dof_map
-
-
-def unroll_dofmap(dofs: npt.NDArray[np.int32], bs: int) -> npt.NDArray[np.int32]:
-    """
-    Given a two-dimensional dofmap of size `(num_cells, num_dofs_per_cell)`
-    Expand the dofmap by its block size such that the resulting array
-    is of size `(num_cells, bs*num_dofs_per_cell)`
-    """
-    num_cells, num_dofs_per_cell = dofs.shape
-    unrolled_dofmap = np.repeat(dofs, bs).reshape(num_cells, num_dofs_per_cell * bs) * bs
-    unrolled_dofmap += np.tile(np.arange(bs), num_dofs_per_cell)
-    return unrolled_dofmap
+import scifem
+import ufl
 
 
 def evaluate_at_vertex_tag(
@@ -72,9 +28,9 @@ def evaluate_at_vertex_tag(
     npt.NDArray[np.int32]
         The values of `u` at the vertices tagged with `tag`
     """
-    v2d = vertex_to_dofmap(u.function_space)
+    v2d = scifem.vertex_to_dofmap(u.function_space)
     block_index = v2d[vt.find(tag)]
-    dofs = unroll_dofmap(block_index.reshape(-1, 1), u.function_space.dofmap.bs)
+    dofs = scifem.utils.unroll_dofmap(block_index.reshape(-1, 1), u.function_space.dofmap.bs)
     return u.x.array[dofs]
 
 
@@ -115,3 +71,32 @@ def gather_broadcast_array(comm, local_array):
     else:
         global_array = np.array([])
     return comm.bcast(global_array, root=0)
+
+
+def matrix_is_zero(A: ufl.core.expr.Expr) -> bool:
+    n = ufl.domain.find_geometric_dimension(A)
+    for i in range(n):
+        for j in range(n):
+            value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(A[i, j] * ufl.dx))
+            print(i, j, value)
+            is_zero = np.isclose(value, 0)
+            if not is_zero:
+                return False
+    return True
+
+
+def float2object(
+    f: float,
+    obj_str: str,
+    mesh: dolfinx.mesh.Mesh,
+    V: dolfinx.fem.FunctionSpace,
+):
+    if obj_str == "float":
+        return f
+    if obj_str == "Constant":
+        return dolfinx.fem.Constant(mesh, f)
+    if obj_str == "Function":
+        v = dolfinx.fem.Function(V)
+        v.x.array[:] = f
+        return v
+    raise ValueError(f"Invalid object string {obj_str!r}")

tests/test_invariants.py

@@ -2,7 +2,6 @@
 import numpy as np
 import pytest
 import ufl
-import utils
 
 from pulse import invariants, kinematics
 
@@ -11,7 +10,7 @@
     "cls, expected",
     (
         (kinematics.DeformationGradient, 4 * 3),
-        (utils.IsochoricDeformationGradient, pow(8, -2 / 3) * 4 * 3),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -2 / 3) * 4 * 3),
     ),
 )
 def test_I1(cls, expected, u) -> None:
@@ -30,7 +29,7 @@ def test_I1(cls, expected, u) -> None:
     (
         (kinematics.DeformationGradient, 0.5 * ((4 * 3) ** 2 - (4**2) * 3)),
         (
-            utils.IsochoricDeformationGradient,
+            kinematics.IsochoricDeformationGradient_from_u,
             0.5 * (((pow(8, -2 / 3) * 4 * 3) ** 2) - ((pow(8, -2 / 3) * 4) ** 2) * 3),
         ),
     ),
@@ -49,7 +48,7 @@ def test_I2(cls, expected, u) -> None:
     "cls, expected",
     (
         (kinematics.DeformationGradient, 4**3),
-        (utils.IsochoricDeformationGradient, (pow(8, -2 / 3) * 4) ** 3),
+        (kinematics.IsochoricDeformationGradient_from_u, (pow(8, -2 / 3) * 4) ** 3),
     ),
 )
 def test_I3(cls, expected, u) -> None:
@@ -67,7 +66,7 @@ def test_I3(cls, expected, u) -> None:
     "cls, expected",
     (
         (kinematics.DeformationGradient, 4),
-        (utils.IsochoricDeformationGradient, pow(8, -2 / 3) * 4),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -2 / 3) * 4),
     ),
 )
 def test_I4(cls, expected, u, mesh) -> None:
@@ -86,7 +85,7 @@ def test_I4(cls, expected, u, mesh) -> None:
     "cls, expected",
     (
         (kinematics.DeformationGradient, 4**2),
-        (utils.IsochoricDeformationGradient, (pow(8, -2 / 3) * 4) ** 2),
+        (kinematics.IsochoricDeformationGradient_from_u, (pow(8, -2 / 3) * 4) ** 2),
     ),
 )
 def test_I5(cls, expected, u, mesh) -> None:
@@ -103,7 +102,7 @@ def test_I5(cls, expected, u, mesh) -> None:
 
 @pytest.mark.parametrize(
     "cls, expected",
-    ((kinematics.DeformationGradient, 0), (utils.IsochoricDeformationGradient, 0)),
+    ((kinematics.DeformationGradient, 0), (kinematics.IsochoricDeformationGradient_from_u, 0)),
 )
 def test_I8(cls, expected, u, mesh) -> None:
     u.interpolate(lambda x: x)

tests/test_kinematics.py

@@ -1,8 +1,7 @@
 import pytest
 import ufl
-import utils
 
-from pulse import kinematics
+from pulse import kinematics, utils
 
 
 def test_SecondOrderIdentity(u) -> None:
@@ -13,7 +12,7 @@ def test_SecondOrderIdentity(u) -> None:
     "cls, factor",
     (
         (kinematics.DeformationGradient, 2),
-        (utils.IsochoricDeformationGradient, pow(8, -1 / 3) * 2),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -1 / 3) * 2),
     ),
 )
 def test_DeformationGradient(cls, factor, u) -> None:
@@ -47,7 +46,7 @@ def test_DeformationGradient(cls, factor, u) -> None:
     "cls, factor",
     (
         (kinematics.DeformationGradient, 4),
-        (utils.IsochoricDeformationGradient, pow(8, -2 / 3) * 4),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -2 / 3) * 4),
     ),
 )
 def test_RightCauchyGreen(cls, factor, u) -> None:
@@ -63,7 +62,7 @@ def test_RightCauchyGreen(cls, factor, u) -> None:
     "cls, factor",
     (
         (kinematics.DeformationGradient, 4),
-        (utils.IsochoricDeformationGradient, pow(8, -2 / 3) * 4),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -2 / 3) * 4),
     ),
 )
 def test_LeftCauchyGreen(cls, factor, u) -> None:
@@ -79,7 +78,7 @@ def test_LeftCauchyGreen(cls, factor, u) -> None:
     "cls, factor",
     (
         (kinematics.DeformationGradient, 1),
-        (utils.IsochoricDeformationGradient, pow(8, -1 / 3) * 2 - 1),
+        (kinematics.IsochoricDeformationGradient_from_u, pow(8, -1 / 3) * 2 - 1),
     ),
 )
 def test_EngineeringStrain(cls, factor, u) -> None:
@@ -95,7 +94,7 @@ def test_EngineeringStrain(cls, factor, u) -> None:
     "cls, factor",
     (
         (kinematics.DeformationGradient, 1.5),
-        (utils.IsochoricDeformationGradient, 0.5 * (pow(8, -2 / 3) * 4 - 1)),
+        (kinematics.IsochoricDeformationGradient_from_u, 0.5 * (pow(8, -2 / 3) * 4 - 1)),
     ),
 )
 def test_GreenLagrangeStrain(cls, factor, u) -> None:

tests/test_material_models.py

@@ -3,7 +3,6 @@
 import dolfinx
 import pytest
 import ufl
-import utils
 
 import pulse
 
@@ -12,7 +11,7 @@
 def test_linear_elastic_model(obj_str, mesh, P1, u) -> None:
     E = 2.0
     _nu = 0.2
-    nu = utils.float2object(f=_nu, obj_str=obj_str, mesh=mesh, V=P1)
+    nu = pulse.utils.float2object(f=_nu, obj_str=obj_str, mesh=mesh, V=P1)
     model = pulse.LinearElastic(E=pulse.Variable(E, "Pa"), nu=nu)
 
     u.interpolate(lambda x: x)
@@ -23,14 +22,14 @@ def test_linear_elastic_model(obj_str, mesh, P1, u) -> None:
     sigma = model.sigma(F)
     I = ufl.Identity(3)
     zero = sigma - (E / (1 + _nu)) * (1 + (_nu / (1 - 2 * _nu)) * 3) * I
-    assert utils.matrix_is_zero(zero)
+    assert pulse.utils.matrix_is_zero(zero)
 
 
 @pytest.mark.parametrize("obj_str", ("float", "Constant", "Function"))
 def test_linear_elastic_model_with_invalid_range(obj_str, mesh, P1) -> None:
     E = pulse.Variable(2.0, "kPa")
     _nu = 0.5
-    nu = utils.float2object(f=_nu, obj_str=obj_str, mesh=mesh, V=P1)
+    nu = pulse.utils.float2object(f=_nu, obj_str=obj_str, mesh=mesh, V=P1)
 
     with pytest.raises(pulse.exceptions.InvalidRangeError):
         pulse.LinearElastic(E=E, nu=nu)