diff --git a/src/fenicsxconcrete/experimental_setup/am_multiple_layers.py b/src/fenicsxconcrete/experimental_setup/am_multiple_layers.py new file mode 100644 index 0000000..213daeb --- /dev/null +++ b/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") diff --git a/tests/experimental_setup/test_experimental_setups.py b/tests/experimental_setup/test_experimental_setups.py index cdee161..2717110 100644 --- a/tests/experimental_setup/test_experimental_setups.py +++ b/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() diff --git a/tests/finite_element_problem/test_am_layers.py b/tests/finite_element_problem/test_am_layers.py new file mode 100644 index 0000000..db174f0 --- /dev/null +++ b/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)