52 improve coverage two small things to fix (#56)

* improve cylinder tests

* imporve point at function

* delete temp mesh files

* remove comment

* rename cylinder test file

* [compression cylinder] use tempfile and gmshio instead of meshio

* attempt to make style check pass

---------

Co-authored-by: Philipp Diercks <[email protected]>

pyproject.toml

@@ -40,3 +40,9 @@ line-length = 119
 
 [tool.isort]
 profile = "black"
+multi_line_output = 3
+include_trailing_comma = true
+force_grid_wrap = 0
+use_parentheses = true
+ensure_newline_before_comments = true
+line_length = 119

src/fenicsxconcrete/experimental_setup/compression_cylinder.py

@@ -1,12 +1,12 @@
-import os
+import tempfile
 from collections.abc import Callable
 
 import dolfinx as df
 import gmsh
-import meshio
 import numpy as np
 import pint
 import ufl
+from dolfinx.io import gmshio
 from mpi4py import MPI
 
 from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions
@@ -31,12 +31,6 @@ def generate_cylinder_mesh(radius: float, height: float, mesh_density: float, el
     mesh : cylinder mesh for dolfin
     """
 
-    # file names and location
-    folder_name = "mesh"
-    file_name = "cylinder"
-    msh_file = f"{folder_name}/{file_name}.msh"
-    xdmf_file = f"{folder_name}/{file_name}.xdmf"
-
     # start gmsh
     gmsh.initialize()
     gmsh.option.setNumber("General.Verbosity", 3)  # only print warnings etc
@@ -46,10 +40,10 @@ def generate_cylinder_mesh(radius: float, height: float, mesh_density: float, el
     # syntax: add_cylinder(x,y,z,dx,dy,dz,radius,angle in radian)
     membrane = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, height, radius, angle=2 * np.pi)
     gmsh.model.occ.synchronize()
-    dim = 3
+    gdim = 3
     # only physical groups get exported
     # syntax: add_physical_group(dim , list of 3d objects, tag)
-    gmsh.model.addPhysicalGroup(dim, [membrane], 1)
+    gmsh.model.addPhysicalGroup(gdim, [membrane], 1)
 
     # meshing
     characteristic_length = height / mesh_density
@@ -59,23 +53,20 @@ def generate_cylinder_mesh(radius: float, height: float, mesh_density: float, el
     # setting the order of the elements
     gmsh.option.setNumber("Mesh.ElementOrder", element_degree)
     gmsh.model.mesh.setOrder(element_degree)
-    gmsh.model.mesh.generate(dim)
-
-    # save it to disk as msh in folder
-    if not os.path.exists(folder_name):  # creat mesh folder if it does nto exists
-        os.mkdir(folder_name)
+    gmsh.model.mesh.generate(gdim)
 
-    # write file
-    gmsh.write(msh_file)
+    # write to tmp file
+    msh_file = tempfile.NamedTemporaryFile(suffix=".msh")
+    gmsh.write(msh_file.name)
     gmsh.finalize()
 
-    # convert msh to xdmf
-    meshio_mesh = meshio.read(msh_file)
-    meshio.write(xdmf_file, meshio_mesh)
+    # reads in the mesh on a single process
+    # and then distributes the cells over available ranks
+    # returns mesh, cell_tags, facet_tags
+    mesh, _, _ = gmshio.read_from_msh(msh_file.name, MPI.COMM_WORLD, gdim=gdim)
 
-    # read xdmf as dolfin mesh
-    with df.io.XDMFFile(MPI.COMM_WORLD, xdmf_file, "r") as mesh_file:
-        mesh = mesh_file.read_mesh(name="Grid")
+    # tmp file is deleted when closed
+    msh_file.close()
 
     return mesh
 
@@ -240,9 +231,6 @@ def create_displacement_boundary(self, V: df.fem.FunctionSpace) -> list[df.fem.b
                 bc_generator.add_dirichlet_bc(np.float64(0.0), point_at(x_min_boundary_point), 1, "geometrical", 1)
                 bc_generator.add_dirichlet_bc(np.float64(0.0), point_at(x_max_boundary_point), 1, "geometrical", 1)
                 bc_generator.add_dirichlet_bc(np.float64(0.0), point_at(y_boundary_point), 0, "geometrical", 0)
-
-            else:
-                raise ValueError(f"dim setting: {self.p['dim']}, not implemented for cylinder bc setup: free")
         else:
             raise ValueError(f"Wrong boundary setting: {self.p['bc_setting']}, for cylinder setup")
 

tests/boundary_conditions/test_bcs_get_boundary_dofs.py

@@ -4,10 +4,7 @@
 import numpy as np
 from mpi4py import MPI
 
-from fenicsxconcrete.boundary_conditions.bcs import (
-    BoundaryConditions,
-    get_boundary_dofs,
-)
+from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions, get_boundary_dofs
 from fenicsxconcrete.boundary_conditions.boundary import plane_at
 
 

tests/finite_element_problem/test_linear_cylinder.py

@@ -0,0 +1,93 @@
+import numpy as np
+import pint
+import pytest
+
+from fenicsxconcrete.experimental_setup.compression_cylinder import CompressionCylinder
+from fenicsxconcrete.finite_element_problem.linear_elasticity import LinearElasticity
+from fenicsxconcrete.helper import Parameters
+from fenicsxconcrete.sensor_definition.base_sensor import Sensor
+from fenicsxconcrete.sensor_definition.other_sensor import ReactionForceSensorBottom
+from fenicsxconcrete.unit_registry import ureg
+
+
+def simple_setup(p: Parameters, displacement: float, sensor: Sensor, bc_setting: pint.Quantity) -> None:
+    parameters = Parameters()  # using the current default values
+
+    parameters["log_level"] = "WARNING" * ureg("")
+    parameters["bc_setting"] = bc_setting
+    parameters["mesh_density"] = 10 * ureg("")
+    parameters = parameters + p
+
+    experiment = CompressionCylinder(parameters)
+
+    problem = LinearElasticity(experiment, parameters)
+    problem.add_sensor(sensor)
+
+    problem.experiment.apply_displ_load(displacement)
+
+    problem.solve()  # solving this
+
+    # last measurement
+    return problem.sensors[sensor.name].data[-1]
+
+
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("degree", [1, 2])
+@pytest.mark.parametrize("bc_setting", ["fixed", "free"])
+def test_force_response(bc_setting: int, degree: int, dim: str) -> None:
+    p = Parameters()  # using the current default values
+
+    p["E"] = 1023 * ureg("MPa")
+    p["nu"] = 0.0 * ureg("")
+    p["radius"] = 0.006 * ureg("m")
+    p["height"] = 0.012 * ureg("m")
+    displacement = -0.003 * ureg("m")
+    p["dim"] = dim * ureg("")
+    p["bc_setting"] = bc_setting * ureg("")
+    p["degree"] = degree * ureg("")
+
+    sensor = ReactionForceSensorBottom()
+    measured = simple_setup(p, displacement, sensor, p["bc_setting"])
+
+    result = None
+    if dim == 2:
+        result = p["E"] * p["radius"] * 2 * displacement / p["height"]
+    elif dim == 3:
+        result = p["E"] * np.pi * p["radius"] ** 2 * displacement / p["height"]
+
+    assert measured == pytest.approx(result.magnitude, 0.01)
+
+
+@pytest.mark.parametrize("bc_setting", ["fixed", "free"])
+def test_errors_dimensions(bc_setting: str) -> None:
+    p = Parameters()  # using the current default values
+    p["E"] = 1023 * ureg("MPa")
+    p["nu"] = 0.0 * ureg("")
+    p["radius"] = 0.006 * ureg("m")
+    p["height"] = 0.012 * ureg("m")
+    displacement = -0.003 * ureg("m")
+    p["bc_setting"] = bc_setting * ureg("")
+    p["degree"] = 2 * ureg("")
+
+    sensor = ReactionForceSensorBottom()
+
+    with pytest.raises(ValueError):
+        p["dim"] = 4 * ureg("")
+        measured = simple_setup(p, displacement, sensor, p["bc_setting"])
+
+
+def test_errors_bc_setting() -> None:
+    p = Parameters()  # using the current default values
+    p["E"] = 1023 * ureg("MPa")
+    p["nu"] = 0.0 * ureg("")
+    p["radius"] = 0.006 * ureg("m")
+    p["height"] = 0.012 * ureg("m")
+    displacement = -0.003 * ureg("m")
+    p["dim"] = 3 * ureg("")
+    p["degree"] = 2 * ureg("")
+
+    sensor = ReactionForceSensorBottom()
+
+    with pytest.raises(ValueError):
+        p["bc_setting"] = "wrong" * ureg("")
+        measured = simple_setup(p, displacement, sensor, p["bc_setting"])