54 add am experimental set up (#68)

* add am experiment set up

* check coverage

* add pseudo test for am set-up

* delete unnecessary parts

* add comment to plane_strain / plane_stress parameter

---------

Co-authored-by: aradermacher <[email protected]>

src/fenicsxconcrete/experimental_setup/am_multiple_layers.py

@@ -0,0 +1,152 @@
+import logging
+from collections.abc import Callable
+
+import dolfinx as df
+import numpy as np
+import pint
+import ufl
+from mpi4py import MPI
+from petsc4py.PETSc import ScalarType
+
+from fenicsxconcrete.boundary_conditions.bcs import BoundaryConditions
+from fenicsxconcrete.boundary_conditions.boundary import plane_at
+from fenicsxconcrete.experimental_setup.base_experiment import Experiment
+from fenicsxconcrete.helper import Parameters
+from fenicsxconcrete.unit_registry import ureg
+
+
+class AmMultipleLayers(Experiment):
+    """sets up a simple layered structure
+
+    all layers of the same height are on top of each other, the boundary on the bottom is fixed
+    the mesh includes all (activation via pseudo-density)
+
+    """
+
+    def __init__(self, parameters: dict[str, pint.Quantity]):
+        """defines default parameters, for the rest, see base class"""
+
+        # initialize default parameters for the setup
+        default_p = Parameters()
+        # default_p['dummy'] = 'example' * ureg('')  # example default parameter for this class
+
+        # updating parameters, overriding defaults
+        default_p.update(parameters)
+        self.logger = logging.getLogger(__name__)
+        super().__init__(default_p)
+
+    @staticmethod
+    def default_parameters() -> dict[str, pint.Quantity]:
+        '''set up a set of working values as example"""
+
+        Returns: dictionary with required parameter
+
+        '''
+        # this must de defined in each setup class
+
+        setup_parameters = {}
+        # geometry
+        setup_parameters["dim"] = 2 * ureg("")
+        setup_parameters["num_layers"] = 10 * ureg("")  # number of layers in y
+        setup_parameters["layer_length"] = 0.5 * ureg("m")  # x_dimension
+        setup_parameters["layer_height"] = 0.01 * ureg("m")  # Dy dimension
+        # only relevant for 3D case [z-dimension]
+        setup_parameters["layer_width"] = 0.05 * ureg("m")
+
+        # mesh
+        setup_parameters["num_elements_layer_length"] = 10 * ureg("")
+        setup_parameters["num_elements_layer_height"] = 1 * ureg("")
+        # only relevant for 3D case
+        setup_parameters["num_elements_layer_width"] = 2 * ureg("")
+
+        # only relevant for 2D case
+        # setup_parameters["stress_case"] = "plane_strain" # not yet implemented
+
+        return setup_parameters
+
+    def setup(self) -> None:
+        """define the mesh for 2D and 3D"""
+        self.logger.debug("setup mesh for %s", self.p["dim"])
+        print(self.p["dim"])
+
+        if self.p["dim"] == 2:
+            self.mesh = df.mesh.create_rectangle(
+                comm=MPI.COMM_WORLD,
+                points=[(0.0, 0.0), (self.p["layer_length"], self.p["num_layer"] * self.p["layer_height"])],
+                n=(self.p["num_elements_layer_length"], self.p["num_layer"] * self.p["num_elements_layer_height"]),
+                cell_type=df.mesh.CellType.quadrilateral,
+            )
+        elif self.p["dim"] == 3:
+            self.mesh = df.mesh.create_box(
+                comm=MPI.COMM_WORLD,
+                points=[
+                    (0.0, 0.0, 0.0),
+                    (self.p["layer_length"], self.p["layer_width"], self.p["num_layer"] * self.p["layer_height"]),
+                ],
+                n=[
+                    self.p["num_elements_layer_length"],
+                    self.p["num_elements_layer_width"],
+                    self.p["num_layer"] * self.p["num_elements_layer_height"],
+                ],
+                cell_type=df.mesh.CellType.hexahedron,
+            )
+        else:
+            raise ValueError(f'wrong dimension: {self.p["dim"]} is not implemented for problem setup')
+
+    def create_displacement_boundary(self, V: df.fem.FunctionSpace) -> list[df.fem.bcs.DirichletBCMetaClass]:
+        """define displacement boundary as fixed at bottom
+
+        Args:
+            V: function space
+
+        Returns: list of dirichlet boundary conditions
+
+        """
+        #
+        bc_generator = BoundaryConditions(self.mesh, V)
+
+        if self.p["dim"] == 2:
+            # fix dofs at bottom
+            bc_generator.add_dirichlet_bc(
+                np.array([0.0, 0.0], dtype=ScalarType),
+                boundary=self.boundary_bottom(),
+                method="geometrical",
+            )
+
+        elif self.p["dim"] == 3:
+            # fix dofs at bottom
+            bc_generator.add_dirichlet_bc(
+                np.array([0.0, 0.0, 0.0], dtype=ScalarType),
+                boundary=self.boundary_bottom(),
+                method="geometrical",
+            )
+
+        return bc_generator.bcs
+
+    def create_body_force(self, v: ufl.argument.Argument) -> ufl.form.Form:
+        """apply body force
+
+        Args:
+            v: test function
+
+        Returns: form for body load
+
+        """
+        force_vector = np.zeros(self.p["dim"])
+        force_vector[-1] = -self.p["rho"] * self.p["g"]
+
+        f = df.fem.Constant(self.mesh, ScalarType(force_vector))
+        L = ufl.dot(f, v) * ufl.dx
+
+        return L
+
+    def boundary_bottom(self) -> Callable:
+        """specify boundary: plane at bottom
+
+        Returns: fct defining if dof is at boundary
+
+        """
+        if self.p["dim"] == 2:
+            return plane_at(0.0, "y")
+        elif self.p["dim"] == 3:
+            return plane_at(0.0, "z")

tests/experimental_setup/test_experimental_setups.py

@@ -2,6 +2,7 @@
 
 import pytest
 
+from fenicsxconcrete.experimental_setup.am_multiple_layers import AmMultipleLayers
 from fenicsxconcrete.experimental_setup.base_experiment import Experiment
 from fenicsxconcrete.experimental_setup.cantilever_beam import CantileverBeam
 from fenicsxconcrete.experimental_setup.compression_cylinder import CompressionCylinder
@@ -11,7 +12,7 @@
 from fenicsxconcrete.unit_registry import ureg
 
 
-@pytest.mark.parametrize("setup", [CantileverBeam, TensileBeam, SimpleBeam, CompressionCylinder])
+@pytest.mark.parametrize("setup", [CantileverBeam, TensileBeam, SimpleBeam, CompressionCylinder, AmMultipleLayers])
 def test_default_parameters(setup: Experiment) -> None:
     """This function creates experimental setups with the respective default dictionaries
 
@@ -34,7 +35,7 @@ def test_default_parameters(setup: Experiment) -> None:
 
 
 # to imporve coverage, I want to test the error messages
-@pytest.mark.parametrize("setup", [CantileverBeam, TensileBeam, SimpleBeam, CompressionCylinder])
+@pytest.mark.parametrize("setup", [CantileverBeam, TensileBeam, SimpleBeam, CompressionCylinder, AmMultipleLayers])
 def test_default_parameters(setup: Experiment) -> None:
     setup_parameters = setup.default_parameters()
 

tests/finite_element_problem/test_am_layers.py

@@ -0,0 +1,77 @@
+import os
+from pathlib import Path
+
+import pytest
+
+from fenicsxconcrete.experimental_setup.am_multiple_layers import AmMultipleLayers
+from fenicsxconcrete.finite_element_problem.linear_elasticity import LinearElasticity
+from fenicsxconcrete.helper import Parameters
+from fenicsxconcrete.sensor_definition.other_sensor import ReactionForceSensorBottom
+from fenicsxconcrete.unit_registry import ureg
+
+
+def set_test_parameters(dim: int) -> Parameters:
+    """set up a test parameter set
+
+    Args:
+        dim: dimension of problem
+
+    Returns: filled instance of Parameters
+
+    """
+    setup_parameters = Parameters()
+
+    setup_parameters["dim"] = dim * ureg("")
+    # setup_parameters["stress_state"] = "plane_strain"
+    setup_parameters["num_layer"] = 5 * ureg("")  # changed in single layer test!!
+    setup_parameters["layer_height"] = 1 / 100 * ureg("m")
+    setup_parameters["layer_length"] = 50 / 100 * ureg("m")
+    setup_parameters["layer_width"] = 5 / 100 * ureg("m")
+
+    setup_parameters["num_elements_layer_length"] = 10 * ureg("")
+    setup_parameters["num_elements_layer_height"] = 1 * ureg("")
+    setup_parameters["num_elements_layer_width"] = 2 * ureg("")
+
+    setup_parameters["rho"] = 2070.0 * ureg("kg/m^3")
+    setup_parameters["E"] = 0.078e6 * ureg("N/m^2")
+    setup_parameters["nu"] = 0.3 * ureg("")
+    # setup_parameters['g'] = 9.81 # in material_problem.py default value
+
+    return setup_parameters
+
+
+@pytest.mark.parametrize("dimension", [2, 3])
+def test_am_single_layer(dimension: int) -> None:
+    """
+    simple test of am experiment set up with dummy linear elastic material to get coverage
+    Note: to be changed if AM MaterilapProblem is implemented"""
+
+    # defining parameters
+    setup_parameters = set_test_parameters(dimension)
+
+    # setting up the problem
+    experiment = AmMultipleLayers(setup_parameters)
+    problem = LinearElasticity(experiment, setup_parameters)
+    problem.add_sensor(ReactionForceSensorBottom())
+
+    # solving and plotting
+    problem.solve()
+    problem.pv_plot()
+
+    # check sensor output
+    force_bottom = problem.sensors["ReactionForceSensorBottom"].data
+
+    dead_load = (
+        problem.parameters["g"]
+        * problem.parameters["rho"]
+        * problem.parameters["layer_length"]
+        * problem.parameters["num_layer"]
+        * problem.parameters["layer_height"]
+    )
+    if dimension == 2:
+        dead_load *= 1 * ureg("m")
+    elif dimension == 3:
+        dead_load *= setup_parameters["layer_width"]
+
+    # dead load of full structure
+    assert sum(force_bottom) == pytest.approx(-dead_load.magnitude)