diff --git a/requirements.txt b/requirements.txt index 184da3639..73c519b40 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,2 +1,2 @@ numpy==1.24.2 -gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@main +gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@vakantie_branch diff --git a/requirements_dev.txt b/requirements_dev.txt index d6af7035a..4bf5a2509 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -1,5 +1,5 @@ numpy==1.24.2 -gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@main +gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@vakantie_branch pytest==7.2.2 pytest-cov==4.0.0 tox==4.4.11 diff --git a/run_stem/demo_create_gmsh_mesh.py b/run_stem/demo_create_gmsh_mesh.py index 54934ca9f..9fb7043c9 100644 --- a/run_stem/demo_create_gmsh_mesh.py +++ b/run_stem/demo_create_gmsh_mesh.py @@ -24,7 +24,6 @@ mesh_output_dir = "./" - gmsh_io = GmshIO() gmsh_io.generate_gmsh_mesh(input_points, extrusion_length, element_size, dims, name_label, mesh_name, mesh_output_dir, diff --git a/setup.cfg b/setup.cfg index 75bd09b72..e66946d5b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -17,7 +17,7 @@ packages = include_package_data = True install_requires = numpy>=1.24 - gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@main + gmsh_utils @ git+https://github.com/StemVibrations/gmsh_utils@vakantie_branch python_requires = >=3.8 [options.extras_require] diff --git a/stem/IO/kratos_additional_processes_io.py b/stem/IO/kratos_additional_processes_io.py index db4fdb66a..759f9419b 100644 --- a/stem/IO/kratos_additional_processes_io.py +++ b/stem/IO/kratos_additional_processes_io.py @@ -42,11 +42,12 @@ def __create_excavation_dict( "python_module": "apply_excavation_process", "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", "process_name": "ApplyExcavationProcess", - "Parameters": parameters.__dict__, + "Parameters": {}, } boundary_dict["Parameters"]["model_part_name"] = f"{self.domain}.{part_name}" boundary_dict["Parameters"]["variable_name"] = "EXCAVATION" + boundary_dict["Parameters"]["deactivate_soil_part"] = parameters.deactivate_body_model_part return boundary_dict diff --git a/stem/IO/kratos_output_io.py b/stem/IO/kratos_output_io.py index b3bd525df..998563c0c 100644 --- a/stem/IO/kratos_output_io.py +++ b/stem/IO/kratos_output_io.py @@ -208,6 +208,7 @@ def __create_json_output_dict( "gauss_points_output_variables": [ op.name for op in output_parameters.gauss_point_results ], + "time_frequency": output_parameters.time_frequency, }, } return output_dict diff --git a/stem/IO/kratos_water_boundaries_io.py b/stem/IO/kratos_water_boundaries_io.py new file mode 100644 index 000000000..5bd94a5a9 --- /dev/null +++ b/stem/IO/kratos_water_boundaries_io.py @@ -0,0 +1,171 @@ +from typing import Dict, Any + +from stem.water_boundaries import WaterBoundary, PhreaticMultiLineBoundary, InterpolateLineBoundary, \ + WaterBoundaryParameters, PhreaticLine + + +class KratosWaterBoundariesIO: + """ + Class to create the water boundary process dictionary for the ProjectParameters.json file in Kratos + + + """ + + def __init__(self, domain: str): + """ + Constructor of KratosWaterBoundariesIO class + + Args: + domain: Name of the Kratos domain + + + """ + self.domain = domain + + def __create_phreatic_line_dict(self, name: str, type: str, water_boundary: PhreaticLine) -> Dict[str, Any]: + """ + Creates a dictionary containing the phreatic line parameters + + + Args: + - name: Name of the boundary + - type: Type of the boundary + - water_boundary: Phreatic line boundary object + + Returns: + - Dict[str, Any]: dictionary containing the phreatic line parameters + + """ + boundary_dict_phreatic_line: Dict[str, Any] = { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": f"{self.domain}.{name}", + "variable_name": "WATER_PRESSURE", + "is_fixed": water_boundary.is_fixed, + "table": [0, 0], + "fluid_pressure_type": type, + "gravity_direction": water_boundary.gravity_direction, + "out_of_plane_direction": water_boundary.out_of_plane_direction, + "specific_weight": water_boundary.specific_weight, + "first_reference_coordinate": water_boundary.first_reference_coordinate, + "second_reference_coordinate": water_boundary.second_reference_coordinate, + "value": water_boundary.value, + } + } + return boundary_dict_phreatic_line + + def __create_phreatic_multi_line_dict(self, name: str, type: str, water_boundary: PhreaticMultiLineBoundary) -> \ + Dict[str, Any]: + """ + Creates a dictionary containing the phreatic multi line parameters + + Args: + - name: Name of the boundary + - type: Type of the boundary + - water_boundary: Multi line phreatic line boundary object + + Returns: + - Dict[str, Any]: dictionary containing the phreatic line parameters + + """ + parameters: Dict[str, Any] = { + "model_part_name": f"{self.domain}.{name}", + "variable_name": "WATER_PRESSURE", + "table": [0, 0, 0], + "value": water_boundary.water_pressure, + "is_fixed": water_boundary.is_fixed, + "gravity_direction": water_boundary.gravity_direction, + "out_of_plane_direction": water_boundary.out_of_plane_direction, + "fluid_pressure_type": type, + "specific_weight": water_boundary.specific_weight, + "x_coordinates": water_boundary.x_coordinates, + "y_coordinates": water_boundary.y_coordinates, + "z_coordinates": water_boundary.z_coordinates, + } + boundary_dict_multi_line: Dict[str, Any] = { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": parameters, + } + return boundary_dict_multi_line + + def __create_interpolation_line_dict(self, name: str, type: str, water_boundary: InterpolateLineBoundary) -> Dict[ + str, Any]: + """ + Creates a dictionary containing the interpolation line parameters + + Args: + - name: Name of the boundary + - type: Type of the boundary + - water_boundary: Interpolation line boundary object + + Returns: + - Dict[str, Any]: dictionary containing the phreatic line parameters + + """ + + boundary_dict_interpolate: Dict[str, Any] = { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": f"{self.domain}.{name}", + "variable_name": "WATER_PRESSURE", + "is_fixed": water_boundary.is_fixed, + "table": 0, + "fluid_pressure_type": type, + "gravity_direction": water_boundary.gravity_direction, + "out_of_plane_direction": water_boundary.out_of_plane_direction, + } + } + return boundary_dict_interpolate + + def __create_water_boundary_dict(self, name: str, type: str, water_boundary: WaterBoundaryParameters) -> Dict[ + str, Any]: + """ + Creates a dictionary containing the water boundary parameters + + Args: + - name: name of the water boundary + - type: type of the water boundary + - water_boundary: water boundary object + + Returns: + - Dict[str, Any]: dictionary containing the water boundary parameters + + """ + if isinstance(water_boundary, PhreaticMultiLineBoundary): + temp_phreatic_multi_line: PhreaticMultiLineBoundary = water_boundary + boundary_dict_multi_line: Dict[str, Any] = self.__create_phreatic_multi_line_dict(name, type, + temp_phreatic_multi_line) + return boundary_dict_multi_line + elif isinstance(water_boundary, InterpolateLineBoundary): + temp_interpolate_line: InterpolateLineBoundary = water_boundary + boundary_dict_interpolate_line: Dict[str, Any] = self.__create_interpolation_line_dict(name, type, + temp_interpolate_line) + return boundary_dict_interpolate_line + elif isinstance(water_boundary, PhreaticLine): + temp_phreatic_line: PhreaticLine = water_boundary + boundary_dict_phreatic_line: Dict[str, Any] = self.__create_phreatic_line_dict(name, type, + temp_phreatic_line) + return boundary_dict_phreatic_line + else: + raise NotImplementedError("This type of boundary is not implemented") + + def create_water_boundary_dict(self, water_boundary: WaterBoundary) -> Dict[str, Any]: + """ + Creates a dictionary containing the water boundary parameters + + Args: + - water_boundary: water boundary object + + Returns: + - Dict[str, Any]: dictionary containing the water boundary parameters + + """ + return self.__create_water_boundary_dict(water_boundary.name, + water_boundary.water_boundary.type, + water_boundary.water_boundary) diff --git a/stem/additional_processes.py b/stem/additional_processes.py index c43d7432c..26b07cef7 100644 --- a/stem/additional_processes.py +++ b/stem/additional_processes.py @@ -20,7 +20,7 @@ class Excavation(AdditionalProcessesParametersABC): - :class:`AdditionalProcessesParametersABC` Attributes: - - deactivate_soil_part (bool): Deactivate or not the body model part + - deactivate_body_model_part (bool): Deactivate or not the body model part """ - deactivate_soil_part: bool + deactivate_body_model_part: bool diff --git a/stem/geometry.py b/stem/geometry.py index bf083cb43..cf6e7f86d 100644 --- a/stem/geometry.py +++ b/stem/geometry.py @@ -29,7 +29,7 @@ class Point(GeometricalObjectABC): Attributes: - __id (int): A unique identifier for the point. - - coordinates (List[float]): An iterable of floats representing the x, y and z coordinates of the point. + - coordinates (Sequence[float]): A sequence of floats representing the x, y and z coordinates of the point. """ def __init__(self, id: int): """ @@ -39,7 +39,24 @@ def __init__(self, id: int): id (int): The id of the point. """ self.__id: int = id - self.coordinates: List[float] = [] + self.coordinates: Sequence[float] = [] + + @classmethod + def create(cls, coordinates: Sequence[float], id: int): + """ + Creates a point object from a list of coordinates and a point id. + + Args: + - coordinates (Sequence[float]): An iterable of floats representing the x, y and z coordinates of the point. + - id (int): The id of the point. + + Returns: + - :class:`Point`: A point object. + + """ + point = cls(id) + point.coordinates = coordinates + return point @property def id(self) -> int: @@ -73,7 +90,7 @@ class Line(GeometricalObjectABC): Attributes: - id (int): A unique identifier for the line. - - point_ids (List[int]): An Iterable of two integers representing the ids of the points that make up the\ + - point_ids (Sequence[int]): A sequence of two integers representing the ids of the points that make up the\ line. """ @@ -85,7 +102,25 @@ def __init__(self, id: int): id (int): The id of the line. """ self.__id: int = id - self.point_ids: List[int] = [] + self.point_ids: Sequence[int] = [] + + @classmethod + def create(cls, point_ids: Sequence[int], id: int): + """ + Creates a line object from a list of point ids and a line id. + + Args: + - point_ids (Sequence[int]): A sequence of two integers representing the ids of the points that make up the\ + line. + - id (int): The id of the line. + + Returns: + - :class:`Line`: A line object. + + """ + line = cls(id) + line.point_ids = point_ids + return line @property def id(self) -> int: @@ -118,12 +153,12 @@ class Surface(GeometricalObjectABC): Attributes: - __id (int): A unique identifier for the surface. - - line_ids (List[int]): An Iterable of three or more integers representing the ids of the lines that make\ + - line_ids (Sequence[int]): A sequence of three or more integers representing the ids of the lines that make\ up the surface. """ def __init__(self, id: int): self.__id: int = id - self.line_ids: List[int] = [] + self.line_ids: Sequence[int] = [] @property def id(self) -> int: @@ -146,6 +181,24 @@ def id(self, value: int): """ self.__id = value + @classmethod + def create(cls, line_ids: Sequence[int], id: int): + """ + Creates a surface object from a list of line ids and a surface id. + + Args: + - line_ids (Sequence[int]): A sequence of three or more integers representing the ids of the lines that make\ + up the surface. + - id (int): The id of the surface. + + Returns: + - :class:`Surface`: A surface object. + + """ + surface = cls(id) + surface.line_ids = line_ids + return surface + class Volume(GeometricalObjectABC): """ @@ -156,12 +209,12 @@ class Volume(GeometricalObjectABC): Attributes: - __id (int): A unique identifier for the volume. - - surface_ids (List[int]): An Iterable of four or more integers representing the ids of the surfaces that\ + - surface_ids (Sequence[int]): A sequence of four or more integers representing the ids of the surfaces that\ make up the volume. """ def __init__(self, id: int): self.__id: int = id - self.surface_ids: List[int] = [] + self.surface_ids: Sequence[int] = [] @property def id(self) -> int: @@ -184,23 +237,42 @@ def id(self, value: int): """ self.__id = value + @classmethod + def create(cls, surface_ids: Sequence[int], id: int): + """ + Creates a volume object from a list of surface ids and a volume id. + + Args: + - surface_ids (Sequence[int]): A sequence of four or more integers representing the ids of the surfaces that\ + make up the volume. + - id (int): The id of the volume. + + Returns: + - :class:`Volume`: A volume object. + + """ + volume = cls(id) + volume.surface_ids = surface_ids + return volume + class Geometry: """ A class to represent a collection of geometric objects in a zero-, one-, two- or three-dimensional space. Attributes: - - points (Optional[List[:class:`Point`]]): An Iterable of Point objects representing the points in the geometry. - - lines (Optional[List[:class:`Line`]]): An Iterable of Line objects representing the lines in the geometry. - - surfaces (Optional[List[:class:`Surface`]]): An Iterable of Surface objects representing the surfaces in the geometry. - - volumes (Optional[List[:class:`Volume`]]): An Iterable of Volume objects representing the volumes in the geometry. + - points (Dict[int, :class:`Point`]): An dictionary of Point objects representing the points in the geometry. + - lines (Dict[int, :class:`Line`]): A dictionary of Line objects representing the lines in the geometry. + - surfaces (Dict[int, :class:`Surface`]): A dictionary of Surface objects representing the surfaces in the \ + geometry. + - volumes (Dict[int, :class:`Volume`]): A dictionary of Volume objects representing the volumes in the geometry. """ - def __init__(self, points: Optional[List[Point]] = None, lines: Optional[List[Line]] = None, - surfaces: Optional[List[Surface]] = None, volumes: Optional[List[Volume]] = None): - self.points: Optional[List[Point]] = points - self.lines: Optional[List[Line]] = lines - self.surfaces: Optional[List[Surface]] = surfaces - self.volumes: Optional[List[Volume]] = volumes + def __init__(self, points: Dict[int, Point] = {}, lines: Dict[int, Line] = {}, + surfaces: Dict[int, Surface] = {}, volumes: Dict[int, Volume] = {}): + self.points: Dict[int, Point] = points + self.lines: Dict[int, Line] = lines + self.surfaces: Dict[int, Surface] = surfaces + self.volumes: Dict[int, Volume] = volumes @staticmethod def __get_unique_entities_by_ids(entities: Sequence[GeometricalObjectABC]): @@ -211,7 +283,7 @@ def __get_unique_entities_by_ids(entities: Sequence[GeometricalObjectABC]): - entities (Sequence[:class:`GeometricalObjectABC`]): An Sequence of geometrical entities. Returns: - - unique_entities (List[:class:`GeometricalObjectABC`): A list of unique geometrical entities entities. + - Sequence[:class:`GeometricalObjectABC`]: A sequence of unique geometrical entities entities. """ unique_entity_ids = [] @@ -232,13 +304,11 @@ def __set_point(geo_data: Dict[str, Any], point_id: int): - point_id (int): The id of the line to create. Returns: - - point (:class:`Point`): The point object. + - :class:`Point`: The point object. """ # create point - point = Point(point_id) - point.coordinates = geo_data["points"][point.id] - return point + return Point.create(geo_data["points"][point_id],point_id) @staticmethod def __set_line(geo_data: Dict[str,Any], line_id: int): @@ -250,15 +320,15 @@ def __set_line(geo_data: Dict[str,Any], line_id: int): - line_id (int): The id of the line to create. Returns: - - line (:class:`Line`): The line object. + - Tuple[:class:`Line`, Sequence[:class:`Point`]]: The line object and the points that make up the line. """ # Initialise point list points = [] # create line and lower dimensional objects - line = Line(abs(line_id)) - line.point_ids = geo_data["lines"][line.id] + line_id = abs(line_id) + line = Line.create(geo_data["lines"][line_id], line_id) for point_id in line.point_ids: points.append(Geometry.__set_point(geo_data, point_id)) return line, points @@ -273,7 +343,8 @@ def __create_surface(geo_data: Dict[str, Any], surface_id: int): - surface_id (int): The id of the surface to create. Returns: - - surface (:class:`Surface`): The surface object. + - Tuple[:class:`Surface`, Sequence[:class:`Line`], Sequence[:class:`Point`]]: The surface object, \ + the lines that make up the surface and the points that make up the lines. """ # Initialise point and line lists @@ -281,8 +352,8 @@ def __create_surface(geo_data: Dict[str, Any], surface_id: int): lines = [] # create surface and lower dimensional objects - surface = Surface(abs(surface_id)) - surface.line_ids = geo_data["surfaces"][surface.id] + surface_id = abs(surface_id) + surface = Surface.create(geo_data["surfaces"][surface_id], surface_id) for line_id in surface.line_ids: line, line_points = Geometry.__set_line(geo_data, line_id) @@ -291,6 +362,43 @@ def __create_surface(geo_data: Dict[str, Any], surface_id: int): return surface, lines, points + @classmethod + def create_geometry_from_geo_data(cls, geo_data: Dict[str,Any]): + """ + Creates the geometry from gmsh geo_data + + Args: + - geo_data (Dict[str, Any]): A dictionary containing the geometry data as provided by gmsh_utils. + + Returns: + - :class:`Geometry`: The geometry object. + """ + + # initialise geometry dictionaries + points = {} + lines = {} + surfaces = {} + volumes = {} + + # add volumes to geometry + for key, value in geo_data["volumes"].items(): + volumes[key] = Volume.create(value,key) + + # add surfaces to geometry + for key, value in geo_data["surfaces"].items(): + surfaces[key] = Surface.create(value, key) + + # add lines to geometry + for key, value in geo_data["lines"].items(): + lines[key] = Line.create(value,key) + + # add points to geometry + for key, value in geo_data["points"].items(): + points[key] = Point.create(value,key) + + # create the geometry class + return cls(points, lines, surfaces, volumes) + @classmethod def create_geometry_from_gmsh_group(cls, geo_data: Dict[str, Any], group_name: str): """ @@ -301,14 +409,14 @@ def create_geometry_from_gmsh_group(cls, geo_data: Dict[str, Any], group_name: s - group_name (str): The name of the group to create the geometry from. Returns: - - geometry (:class:`Geometry`): A Geometry object containing the geometric objects in the group. + - :class:`Geometry`: A Geometry object containing the geometric objects in the group. """ - # initialize point, line, surface and volume lists - points = [] - lines = [] - surfaces = [] - volumes = [] + # initialize point, line, surface and volume dictionaries + points = {} + lines = {} + surfaces = {} + volumes = {} group_data = geo_data["physical_groups"][group_name] ndim_group = group_data["ndim"] @@ -316,42 +424,41 @@ def create_geometry_from_gmsh_group(cls, geo_data: Dict[str, Any], group_name: s if ndim_group == 0: # create points for id in group_data["geometry_ids"]: - points.append(Geometry.__set_point(geo_data, id)) + points[id] = Geometry.__set_point(geo_data, id) elif ndim_group == 1: # create lines and lower dimensional objects for id in group_data["geometry_ids"]: line, line_points = Geometry.__set_line(geo_data, id) - lines.append(line) - points.extend(line_points) + lines[id] = line + for point in line_points: + points[point.id] = point elif ndim_group == 2: # create surfaces and lower dimensional objects for id in group_data["geometry_ids"]: surface, surface_lines, surface_points = Geometry.__create_surface(geo_data, id) - surfaces.append(surface) - lines.extend(surface_lines) - points.extend(surface_points) + surfaces[id] = surface + for line in surface_lines: + lines[line.id] = line + for point in surface_points: + points[point.id] = point elif ndim_group == 3: # Create volumes and lower dimensional objects for id in group_data["geometry_ids"]: - volume = Volume(id) - volume.surface_ids = geo_data["volumes"][volume.id] + volume = Volume.create(geo_data["volumes"][id], id) # create surfaces and lower dimensional objects which are part of the current volume for surface_id in volume.surface_ids: surface, surface_lines, surface_points = Geometry.__create_surface(geo_data, surface_id) - surfaces.append(surface) - lines.extend(surface_lines) - points.extend(surface_points) - volumes.append(volume) - - # remove duplicates from points, lines, surfaces, volumes - unique_volumes = Geometry.__get_unique_entities_by_ids(volumes) - unique_surfaces = Geometry.__get_unique_entities_by_ids(surfaces) - unique_lines = Geometry.__get_unique_entities_by_ids(lines) - unique_points = Geometry.__get_unique_entities_by_ids(points) - - return cls(unique_points, unique_lines, unique_surfaces, unique_volumes) + surfaces[abs(surface_id)] = surface + for line in surface_lines: + lines[line.id] = line + for point in surface_points: + points[point.id] = point + + volumes[id] = volume + + return cls(points, lines, surfaces, volumes) diff --git a/stem/load.py b/stem/load.py index a013a600f..3d172cfd3 100644 --- a/stem/load.py +++ b/stem/load.py @@ -24,8 +24,8 @@ class PointLoad(LoadParametersABC): - value (List[float]): Entity of the load in the 3 directions [N]. """ - active: List[bool] = field(default_factory=lambda: [True, True, True]) - value: List[float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + active: List[bool] + value: List[float] @dataclass @@ -37,8 +37,8 @@ class LineLoad(LoadParametersABC): - active (List[bool]): Activate/deactivate load for each direction. - value (List[float]): Entity of the load in the 3 directions [N]. """ - active: List[bool] = field(default_factory=lambda: [True, True, True]) - value: List[float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + active: List[bool] + value: List[float] @dataclass @@ -50,8 +50,8 @@ class SurfaceLoad(LoadParametersABC): - active (List[bool]): Activate/deactivate load for each direction. - value (List[float]): Entity of the load in the 3 directions [N]. """ - active: List[bool] = field(default_factory=lambda: [True, True, True]) - value: List[float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + active: List[bool] + value: List[float] @dataclass @@ -72,10 +72,10 @@ class MovingLoad(LoadParametersABC): - offset (float): Offset of the moving load [m]. """ - load: Union[List[float], List[str]] = field(default_factory=lambda: [0.0, 0.0, 0.0]) - direction: List[float] = field(default_factory=lambda: [1, 1, 1]) - velocity: Union[float, str] = 0.0 - origin: List[float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + load: Union[List[float], List[str]] + direction: List[float] + velocity: Union[float, str] + origin: List[float] offset: float = 0.0 @@ -92,5 +92,5 @@ class GravityLoad(LoadParametersABC): - value (List[float]): Entity of the gravity acceleration in the 3 directions [m/s^2]. Should be -9.81 only in the vertical direction """ - active: List[bool] = field(default_factory=lambda: [False, False, False]) - value: List[float] = field(default_factory=lambda: [0.0, 0.0, 0.0]) + active: List[bool] + value: List[float] diff --git a/stem/mesh.py b/stem/mesh.py index 6ea3323b4..329c20b10 100644 --- a/stem/mesh.py +++ b/stem/mesh.py @@ -1,53 +1,109 @@ -from typing import Dict, List, Tuple, Union, Any +from typing import Dict, List, Tuple, Sequence, Union, Any, Optional +from enum import Enum +from dataclasses import dataclass + import numpy as np import numpy.typing as npt from stem.IO.kratos_io import KratosIO -class Node: +class ElementShape(Enum): """ - Class containing information about a node + Enum class for the element shape. TRIANGLE for triangular elements and tetrahedral elements, QUADRILATERAL for + quadrilateral elements and hexahedral elements. - Attributes: - - id (int): node id - - coordinates (np.array): node coordinates + """ + TRIANGLE = "triangle" + QUADRILATURAL = "quadrilateral" + +class MeshSettings: """ - def __init__(self, id, coordinates): - self.id = id - self.coordinates = coordinates + A class to represent the mesh settings. -class Element: + Attributes: + - element_size (float): The element size (default -1, which means that gmsh determines the size). + - element_shape (:class:`stem.model.ElementShape`): The element shape. TRIANGLE for triangular elements and \ + tetrahedral elements, QUADRILATERAL for quadrilateral elements and hexahedral elements. (default TRIANGLE) + - __element_order (int): The element order. 1 for linear elements, 2 for quadratic elements. (default 1) """ - Class containing information about an element + + def __init__(self, element_size: float = -1, element_order: int = 1, + element_shape: ElementShape = ElementShape.TRIANGLE): + """ + Initialize the mesh settings. + + Args: + - element_size (float): The element size (default -1, which means that gmsh determines the size). + - element_order (int): The element order. 1 for linear elements, 2 for quadratic elements. (default 1) + - element_shape (:class:`stem.model.ElementShape`): The element shape. TRIANGLE for triangular elements and \ + tetrahedral elements, QUADRILATERAL for quadrilateral elements and hexahedral elements. (default TRIANGLE) + """ + self.element_size: float = element_size + self.element_shape: ElementShape = element_shape + + if element_order not in [1, 2]: + raise ValueError("The element order must be 1 or 2. Higher order elements are not supported.") + + self.__element_order: int = element_order + + @property + def element_order(self): + """ + Get the element order. + + Returns: + - int: element order + """ + return self.__element_order + + @element_order.setter + def element_order(self, element_order: int): + """ + Set the element order. The element order must be 1 or 2. + + Args: + - element_order (int): element order + + Raises: + - ValueError: If the element order is not 1 or 2. + """ + + if element_order not in [1, 2]: + raise ValueError("The element order must be 1 or 2. Higher order elements are not supported.") + + self.__element_order = element_order + + +class Node: + """ + Class containing information about a node Attributes: - - id (int): element id - - element_type (str): element type - - node_ids (Union[List[int], npt.NDArray[np.int64]]): node ids + - id (int): node id + - coordinates (Sequence[float]): node coordinates """ - def __init__(self, id: int, element_type: str, node_ids: Union[List[int], npt.NDArray[np.int64]]): + def __init__(self, id: int, coordinates: Sequence[float]): self.id: int = id - self.element_type: str = element_type - self.node_ids: Union[List[int], npt.NDArray[np.int64]] = node_ids + self.coordinates: Sequence[float] = coordinates -class Condition: +class Element: """ - Class containing information about a condition + Class containing information about an element Attributes: - - id (int): condition id + - id (int): element id - element_type (str): element type - - node_ids (Union[List[int], npt.NDArray[np.int64]]): node ids + - node_ids (Sequence[int]): node ids """ - def __init__(self, id: int, element_type: str, node_ids: Union[List[int], npt.NDArray[np.int64]]): + def __init__(self, id: int, element_type: str, node_ids: Sequence[int]): self.id: int = id self.element_type: str = element_type - self.node_ids: Union[List[int], npt.NDArray[np.int64]] = node_ids + self.node_ids: Sequence[int] = node_ids class Mesh: @@ -59,24 +115,57 @@ class Mesh: Attributes: - ndim (int): number of dimensions of the mesh - - nodes (np.array or None): node id followed by node coordinates in an array - - elements (np.array or None): element id followed by connectivities in an array - - conditions (np.array or None): condition id followed by connectivities in an array + - nodes (List[Node]): node id followed by node coordinates in an array + - elements (List[Element]): element id followed by connectivities in an array """ def __init__(self, ndim: int): self.ndim: int = ndim - self.nodes = None - self.elements = None - self.conditions = None - + self.nodes: List[Node] = [] + self.elements: List[Element] = [] @classmethod - def read_mesh_from_gmsh(cls, mesh_file_name: str) -> None: - #todo implement this method to read mesh from gmsh file and create a mesh object with the data read from the - # file. - pass + def create_mesh_from_gmsh_group(cls, mesh_data: Dict[str, Any], group_name: str): + """ + Creates a mesh object from gmsh group + + Args: + - mesh_data (Dict[str, Any]): dictionary of mesh data + - group_name (str): name of the group + + Raises: + - ValueError: If the group name is not found in the mesh data + + Returns: + - :class:`Mesh`: mesh object + """ + + if group_name not in mesh_data["physical_groups"]: + raise ValueError(f"Group {group_name} not found in mesh data") + + # create mesh object + group_data = mesh_data["physical_groups"][group_name] + + group_element_ids = group_data["element_ids"] + group_node_ids = group_data["node_ids"] + group_element_type = group_data["element_type"] + + element_type_data = mesh_data["elements"][group_element_type] + + # create element per element id + elements = [Element(element_id, group_element_type, element_type_data[element_id]) + for element_id in group_element_ids] + + # create node per node id + nodes = [Node(node_id, mesh_data["nodes"][node_id]) for node_id in group_node_ids] + + # add nodes and elements to mesh object + mesh = cls(group_data["ndim"]) + mesh.nodes = nodes + mesh.elements = elements + + return mesh def prepare_data_for_kratos(self, mesh_data: Dict[str, Any]) \ -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.int64]]: @@ -104,7 +193,6 @@ def prepare_data_for_kratos(self, mesh_data: Dict[str, Any]) \ return nodes, all_elements - def write_mesh_to_kratos_structure(self, mesh_data: Dict[str, Any], filename: str) -> None: """ Writes mesh data to the structure which can be read by Kratos diff --git a/stem/model.py b/stem/model.py index 4aec9b724..88fc86705 100644 --- a/stem/model.py +++ b/stem/model.py @@ -1,6 +1,21 @@ -from typing import List +from enum import Enum +from dataclasses import dataclass +from typing import List, Sequence, Dict, Any, Optional, Union, get_args + +import numpy as np +import numpy.typing as npty + +from gmsh_utils import gmsh_IO from stem.model_part import ModelPart, BodyModelPart +from stem.soil_material import * +from stem.structural_material import * +from stem.boundary import * +from stem.geometry import Geometry +from stem.mesh import Mesh, MeshSettings +from stem.load import * +from stem.solver import Problem, StressInitialisationType +from stem.utils import Utils class Model: @@ -8,17 +23,428 @@ class Model: A class to represent the main model. Attributes: + - ndim (int): Number of dimensions of the model - project_parameters (dict): A dictionary containing the project parameters. - - solver (Solver): The solver used to solve the problem. - - body_model_parts (list): A list containing the body model parts. - - process_model_parts (list): A list containing the process model parts. + - solver (:class:`stem.solver.Solver`): The solver used to solve the problem. + - geometry (Optional[:class:`stem.geometry.Geometry`]) The geometry of the whole model. + - body_model_parts (List[:class:`stem.model_part.BodyModelPart`]): A list containing the body model parts. + - process_model_parts (List[:class:`stem.model_part.ModelPart`]): A list containing the process model parts. + - extrusion_length (Optional[Sequence[float]]): The extrusion length in x, y and z direction """ - def __init__(self): - - self.project_parameters = None + def __init__(self, ndim: int): + self.ndim: int = ndim + self.project_parameters: Optional[Problem] = None self.solver = None + self.geometry: Optional[Geometry] = None + self.mesh_settings: MeshSettings = MeshSettings() + self.gmsh_io = gmsh_IO.GmshIO() self.body_model_parts: List[BodyModelPart] = [] self.process_model_parts: List[ModelPart] = [] + self.extrusion_length: Optional[Sequence[float]] = None + + def __del__(self): + """ + Destructor of the Model class. Finalizes the gmsh_io instance. + + """ + self.gmsh_io.finalize_gmsh() + + def __get_geometry_from_geo_data(self, geo_data: Dict[str, Any]): + """ + Get the geometry from the geo_data as generated by gmsh_io. + + Args: + - geo_data (Dict[str, Any]): dictionary containing the geometry data as generated by the gmsh_io + + """ + + self.geometry = Geometry.create_geometry_from_geo_data(geo_data) + + def add_all_layers_from_geo_file(self, geo_file_name: str, body_names: Sequence[str]): + """ + Add all physical groups from a geo file to the model. The physical groups with the names in body_names are + added as body model parts, the other physical groups are added as process model parts. + + Args: + - geo_file_name (str): name of the geo file + - body_names (Sequence[str]): names of the physical groups which should be added as body model parts + + """ + + # read the geo file and generate the geo_data dictionary + self.gmsh_io.read_gmsh_geo(geo_file_name) + + # Reset the gmsh instance with the geo data, as read from the geo file + self.gmsh_io.generate_geo_from_geo_data() + + geo_data = self.gmsh_io.geo_data + + # Create geometry and model part for each physical group in the gmsh geo_data + model_part: Union[ModelPart, BodyModelPart] + for group_name in geo_data["physical_groups"].keys(): + + # create model part, if the group name is in the body names, create a body model part, otherwise a process + # model part + if group_name in body_names: + model_part = BodyModelPart(group_name) + else: + model_part = ModelPart(group_name) + + # set the name and geometry of the model part + model_part.get_geometry_from_geo_data(geo_data, group_name) + + # add model part to either body model parts or process model part + if isinstance(model_part, BodyModelPart): + self.body_model_parts.append(model_part) + else: + self.process_model_parts.append(model_part) + + def add_soil_layer_by_coordinates(self, coordinates: Sequence[Sequence[float]], + material_parameters: Union[SoilMaterial, StructuralMaterial], name: str, + ): + """ + Adds a soil layer to the model by giving a sequence of 2D coordinates. In 3D the 2D geometry is extruded in + the direction of the extrusion_length + + Args: + - coordinates (Sequence[Sequence[float]]): The plane coordinates of the soil layer. + - material_parameters (Union[:class:`stem.soil_material.SoilMaterial`, \ + :class:`stem.structural_material.StructuralMaterial`]): The material parameters of the soil layer. + - name (str): The name of the soil layer. + + Raises: + - ValueError: if extrusion_length is not specified. + """ + + # sort coordinates in anti-clockwise order, such that elements in mesh are also in anti-clockwise order + if Utils.are_2d_coordinates_clockwise(coordinates): + coordinates = coordinates[::-1] + + gmsh_input = {name: {"coordinates": coordinates, "ndim": self.ndim}} + # check if extrusion length is specified in 3D + if self.ndim == 3: + if self.extrusion_length is None: + raise ValueError("Extrusion length must be specified for 3D models") + + gmsh_input[name]["extrusion_length"] = self.extrusion_length + + # todo check if this function in gmsh io can be improved + self.gmsh_io.generate_geometry(gmsh_input, "") + + # create body model part + body_model_part = BodyModelPart(name) + body_model_part.material = material_parameters + + # set the geometry of the body model part + body_model_part.get_geometry_from_geo_data(self.gmsh_io.geo_data, name) + + self.body_model_parts.append(body_model_part) + + def add_load_by_coordinates(self, coordinates: Sequence[Sequence[float]], load_parameters: LoadParametersABC, name: str): + """ + Adds a load to the model by giving a sequence of 3D coordinates. For a 2D model, the third coordinate is + ignored. + + Args: + - coordinates (Sequence[Sequence[float]]): The coordinates of the load. + - load_parameters (:class:`stem.load.LoadParametersABC`): The parameters of the load. + - name (str): The name of the load part. + + Raises: + - ValueError: if load_parameters is not of one of the classes PointLoad, MovingLoad, LineLoad + or SurfaceLoad. + """ + + # todo add validation that load is applied on a body model part + + # validation of inputs + self.validate_coordinates(coordinates) + if isinstance(load_parameters, MovingLoad): + self.__validate_moving_load_parameters(coordinates, load_parameters) + + # create input for gmsh + if isinstance(load_parameters, PointLoad): + gmsh_input = {name: {"coordinates": coordinates, "ndim": 0}} + elif isinstance(load_parameters, LineLoad) or isinstance(load_parameters, MovingLoad): + gmsh_input = {name: {"coordinates": coordinates, "ndim": 1}} + elif isinstance(load_parameters, SurfaceLoad): + gmsh_input = {name: {"coordinates": coordinates, "ndim": 2}} + else: + raise ValueError(f'Invalid load_parameters ({load_parameters.__class__.__name__}) object' + f' provided for the load {name}. Expected one of PointLoad, MovingLoad,' + f' LineLoad or SurfaceLoad.') + + self.gmsh_io.generate_geometry(gmsh_input, "") + + # create model part + model_part = ModelPart(name) + model_part.parameters = load_parameters + + # set the geometry of the model part + model_part.get_geometry_from_geo_data(self.gmsh_io.geo_data, name) + + self.process_model_parts.append(model_part) + + @staticmethod + def validate_coordinates(coordinates: Union[Sequence[Sequence[float]], npty.NDArray[np.float64]]): + """ + Validates the coordinates in input. + + Args: + - coordinates (Sequence[Sequence[float]]): The coordinates of the load. + + Raises: + - ValueError: if coordinates is not convertible to a 2D array (i.e. a sequence of sequences) + - ValueError: if the number of elements (number of coordinates) is not 3. + """ + + # if is not an array, make it array! + + if not isinstance(coordinates, np.ndarray): + coordinates = np.array(coordinates) + + if len(coordinates.shape) != 2: + raise ValueError(f"Coordinates are not a sequence of a sequence or a 2D array.") + + if coordinates.shape[1] != 3: + raise ValueError(f"Coordinates should be 3D but {coordinates.shape[1]} coordinates were given.") + + @staticmethod + def __validate_moving_load_parameters(coordinates: Sequence[Sequence[float]], load_parameters: MovingLoad): + """ + Validates the coordinates in input for the moving load and the trajectory (collinearity of the + points and if the origin is between the point). + + Args: + - coordinates (Sequence[Sequence[float]]): The start-end coordinate of the moving load. + - parameters (:class:`stem.load.LoadParametersABC`): The parameters of the load. + + Raises: + - ValueError: if moving load origin is not on trajectory + """ + + # iterate over each line constituting the trajectory + for ix in range(len(coordinates)-1): + + # check origin is collinear to the edges of the line + collinear_check = Utils.is_collinear( + point=load_parameters.origin, start_point=coordinates[ix],end_point=coordinates[ix+1] + ) + # check origin is between the edges of the line (edges included) + is_between_check = Utils.is_point_between_points( + point=load_parameters.origin, start_point=coordinates[ix], end_point=coordinates[ix+1] + ) + # check if point complies + is_on_line = collinear_check and is_between_check + # exit at the first success of the test (point in the line) + if is_on_line: + return + + # none of the lines contain the origin, then raise an error + raise ValueError(f"Origin is not in the trajectory of the moving load.") + + def add_boundary_condition_by_geometry_ids(self, ndim_boundary: int, geometry_ids: Sequence[int], + boundary_parameters: BoundaryParametersABC, name: str): + """ + Add a boundary condition to the model by giving the geometry ids of the boundary condition. + + Args: + - ndim_boundary (int): dimension of the boundary condition + - geometry_ids (Sequence[int]): geometry ids of the boundary condition + - boundary_condition (:class:`stem.boundary_condition.BoundaryCondition`): boundary condition object + - name (str): name of the boundary condition + + """ + + # add physical group to gmsh + self.gmsh_io.add_physical_group(name, ndim_boundary, geometry_ids) + + # create model part + model_part = ModelPart(name) + + # retrieve geometry from gmsh and add to model part + model_part.get_geometry_from_geo_data(self.gmsh_io.geo_data, name) + + # add boundary parameters to model part + model_part.parameters = boundary_parameters + + self.process_model_parts.append(model_part) + + def synchronise_geometry(self): + """ + Synchronise the geometry of all model parts and synchronise the geometry of the whole model. This function + recalculates all ids and connectivities of all geometrical entities. + + """ + + # synchronize gmsh and extract geo data + self.gmsh_io.synchronize_gmsh() + self.gmsh_io.extract_geo_data() + + # collect all model parts + all_model_parts: List[Union[BodyModelPart, ModelPart]] = [] + all_model_parts.extend(self.body_model_parts) + all_model_parts.extend(self.process_model_parts) + + # Get the geometry from the geo_data for each model part + for model_part in all_model_parts: + model_part.get_geometry_from_geo_data(self.gmsh_io.geo_data, model_part.name) + + # get the complete geometry + self.__get_geometry_from_geo_data(self.gmsh_io.geo_data) + + def generate_mesh(self): + """ + Generate the mesh for the whole model. + + """ + + # generate mesh + self.gmsh_io.generate_mesh(self.ndim, element_size=self.mesh_settings.element_size, + order=self.mesh_settings.element_order) + + # collect all model parts + all_model_parts: List[Union[BodyModelPart, ModelPart]] = [] + all_model_parts.extend(self.body_model_parts) + all_model_parts.extend(self.process_model_parts) + + # add the mesh to each model part + for model_part in all_model_parts: + model_part.mesh = Mesh.create_mesh_from_gmsh_group(self.gmsh_io.mesh_data, model_part.name) + + def __validate_model_part_names(self): + """ + Checks if all model parts have a unique name. + + Raises: + - ValueError: If not all model parts have a name. + - ValueError: If not all model part names are unique . + """ + + # collect all model parts + all_model_parts: List[Union[BodyModelPart, ModelPart]] = [] + all_model_parts.extend(self.body_model_parts) + all_model_parts.extend(self.process_model_parts) + + unique_names = [] + for model_part in all_model_parts: + # Check if all model parts have a name + if model_part.name is None: + raise ValueError("All model parts must have a name") + else: + if model_part.name in unique_names: + raise ValueError("All model parts must have a unique name") + unique_names.append(model_part.name) + + def __add_gravity_model_part(self, gravity_load: GravityLoad, ndim: int, geometry_ids: Sequence[int]): + """ + Add a gravity model part to the complete model. + + Args: + - gravity_load (GravityLoad): The gravity load object. + - ndim (int): The number of dimensions of the on which the gravity load should be applied. + - geometry_ids (Sequence[int]): The geometry on which the gravity load should be applied. + + """ + + # set new model part name + model_part_name = f"gravity_load_{ndim}d" + + # create new gravity physical group and model part + self.gmsh_io.add_physical_group(model_part_name, ndim, geometry_ids) + model_part = ModelPart(model_part_name) + + model_part.parameters = gravity_load + + # add gravity load to process model parts + self.process_model_parts.append(model_part) + + def __add_gravity_load(self, gravity_value: float = -9.81, vertical_axis: int = 1): + """ + Add a gravity load to the complete model. + + Args: + - gravity_value (float): The gravity value [m/s^2]. (default -9.81) + - vertical_axis (int): The vertical axis of the model. x=>0, y=>1, z=>2. (default y, 1) + + """ + + # set gravity load at vertical axis + gravity_load_values: List[float] = [0, 0, 0] + gravity_load_values[vertical_axis] = gravity_value + gravity_load = GravityLoad(value=gravity_load_values, active=[True, True, True]) + + # get all body model part names + body_model_part_names = [body_model_part.name for body_model_part in self.body_model_parts] + + # get geometry ids and ndim for each body model part + model_parts_geometry_ids = np.array([self.gmsh_io.geo_data["physical_groups"][name]["geometry_ids"] for name in + body_model_part_names]) + + model_parts_ndim = np.array([self.gmsh_io.geo_data["physical_groups"][name]["ndim"] + for name in body_model_part_names]).ravel() + + # add gravity load as physical group per dimension + body_geometries_1d = model_parts_geometry_ids[model_parts_ndim == 1].ravel() + if len(body_geometries_1d) > 0: + self.__add_gravity_model_part(gravity_load, 1, body_geometries_1d) + + body_geometries_2d = model_parts_geometry_ids[model_parts_ndim == 2].ravel() + if len(body_geometries_2d) > 0: + self.__add_gravity_model_part(gravity_load, 2, body_geometries_2d) + + body_geometries_3d = model_parts_geometry_ids[model_parts_ndim == 3].ravel() + if len(body_geometries_3d) > 0: + self.__add_gravity_model_part(gravity_load, 3, body_geometries_3d) + + self.synchronise_geometry() + + def validate(self): + """ + Validate the model. \ + - Checks if all model parts have a unique name. + + """ + + self.__validate_model_part_names() + + def __setup_stress_initialisation(self): + """ + Set up the stress initialisation. For K0 procedure and gravity loading, a gravity load is added to the model. + + Raises: + - ValueError: If the project parameters are not set. + + """ + + if self.project_parameters is None: + raise ValueError("Project parameters must be set before setting up the stress initialisation") + + # add gravity load if K0 procedure or gravity loading is used + if (self.project_parameters.settings.stress_initialisation_type == + StressInitialisationType.K0_PROCEDURE) or \ + (self.project_parameters.settings.stress_initialisation_type == + StressInitialisationType.GRAVITY_LOADING): + + self.__add_gravity_load() + + def post_setup(self): + """ + Post setup of the model. \ + - Synchronise the geometry. \ + - Generate the mesh. \ + - Validate the model. \ + - Set up the stress initialisation. + + """ + + self.synchronise_geometry() + self.generate_mesh() + self.validate() + + self.__setup_stress_initialisation() + + diff --git a/stem/model_part.py b/stem/model_part.py index 8ba1d5402..f01a48298 100644 --- a/stem/model_part.py +++ b/stem/model_part.py @@ -1,9 +1,14 @@ from typing import Optional, Union, Dict, Any +from stem.additional_processes import AdditionalProcessesParametersABC +from stem.boundary import BoundaryParametersABC +from stem.load import LoadParametersABC from stem.soil_material import SoilMaterial from stem.structural_material import StructuralMaterial -from stem.geometry import Geometry, Volume, Surface, Line, Point +from stem.geometry import Geometry +from stem.mesh import Mesh + class ModelPart: """ @@ -11,21 +16,38 @@ class ModelPart: like excavation. Attributes: - - name (str): name of the model part - - nodes (np.array or None): node id followed by node coordinates in an array - - elements (np.array or None): element id followed by connectivities in an array - - conditions (np.array or None): condition id followed by connectivities in an array + - __name (str): name of the model part - geometry (Optional[:class:`stem.geometry.Geometry`]): geometry of the model part - - parameters (dict): dictionary containing the model part parameters + - parameters (Optional[Union[:class:`stem.load.LoadParametersABC`, \ + :class:`stem.boundary.BoundaryParametersABC, \ + :class:`stem.additional_processes.AdditionalProcessesParametersABC`]]): process parameters containing the \ + model part parameters. + - mesh (Optional[:class:`stem.mesh.Mesh`]): mesh of the model part """ - def __init__(self): - self.name = None - self.nodes = None - self.elements = None - self.conditions = None + def __init__(self, name: str): + """ + Initialize the model part + Args: + - name (str): name of the model part + """ + self.__name: str = name self.geometry: Optional[Geometry] = None - self.parameters = {} + self.parameters: Optional[ + Union[LoadParametersABC, BoundaryParametersABC,AdditionalProcessesParametersABC] + ] = None + self.mesh: Optional[Mesh] = None + + @property + def name(self): + """ + Get the name of the model part + + Returns: + - str: name of the model part + + """ + return self.__name def get_geometry_from_geo_data(self, geo_data: Dict[str, Any], name: str): """ @@ -43,20 +65,25 @@ class BodyModelPart(ModelPart): """ This class contains model parts which are part of the body, e.g. a soil layer or track components. - Inheritance: + Inheritance: - :class:`ModelPart` Attributes: - - name (str): name of the model part - - nodes (np.array or None): node id followed by node coordinates in an array - - elements (np.array or None): element id followed by connectivities in an array - - conditions (np.array or None): condition id followed by connectivities in an array - - parameters (dict): dictionary containing the model part parameters + - __name (str): name of the model part + - geometry (Optional[:class:`stem.geometry.Geometry`]): geometry of the model part + - mesh (Optional[:class:`stem.mesh.Mesh`]): mesh of the model part + - parameters (Dict[str, Any]): dictionary containing the model part parameters - material (Union[:class:`stem.soil_material.SoilMaterial`, \ :class:`stem.structural_material.StructuralMaterial`]): material of the model part """ - def __init__(self): - super().__init__() + def __init__(self, name: str): + """ + Initialize the body model part + + Args: + - name (str): name of the body model part + """ + super().__init__(name) self.material: Optional[Union[SoilMaterial, StructuralMaterial]] = None diff --git a/stem/utils.py b/stem/utils.py new file mode 100644 index 000000000..e24df4335 --- /dev/null +++ b/stem/utils.py @@ -0,0 +1,114 @@ +from typing import Sequence + +import numpy as np + + +class Utils: + """ + Class containing utility methods. + + """ + + @staticmethod + def are_2d_coordinates_clockwise(coordinates: Sequence[Sequence[float]]): + """ + Checks if the 2D coordinates are given in clockwise order. If the signed area is positive, the coordinates + are given in clockwise order. + + Args: + - coordinates (Sequence[Sequence[float]]): coordinates of the points of a surface + + Returns: + - bool: True if the coordinates are given in clockwise order, False otherwise. + """ + + # calculate signed area of polygon + signed_area = 0.0 + for i in range(len(coordinates) - 1): + signed_area += (coordinates[i + 1][0] - coordinates[i][0]) * (coordinates[i + 1][1] + coordinates[i][1]) + + signed_area += (coordinates[0][0] - coordinates[-1][0]) * (coordinates[0][1] + coordinates[-1][1]) + + # if signed area is positive, the coordinates are given in clockwise order + return signed_area > 0.0 + + @staticmethod + def check_dimensions(points:Sequence[Sequence[float]]): + """ + + Check if points have the same dimensions (2D or 3D). + + Args: + - points: (Sequence[Sequence[float]]): points to be tested + + Raises: + - ValueError: when the points have different dimensions. + - ValueError: when the dimension is not either 2 or 3D. + """ + + lengths = [len(point) for point in points] + if len(np.unique(lengths)) != 1: + raise ValueError("Mismatch in dimension of given points!") + + if any([ll not in [2, 3] for ll in lengths]): + raise ValueError("Dimension of the points should be 2D or 3D.") + + @staticmethod + def is_collinear(point: Sequence[float], start_point: Sequence[float], end_point: Sequence[float], + a_tol: float = 1e-06): + """ + Check if point is aligned with the other two on a line. Points must have the same dimension (2D or 3D) + + Args: + - point (Sequence[float]): point coordinates to be tested + - start_point (Sequence[float]): coordinates of first point of a line + - end_point (Sequence[float]): coordinates of second point of a line + - a_tol (float): absolute tolerance to check collinearity (default 1e-6) + + Raises: + - ValueError: when there is a dimension mismatch in the point dimensions. + + Returns: + - bool: whether the point is aligned or not + """ + + # check dimensions of points for validation + Utils.check_dimensions([point, start_point, end_point]) + + vec_1 = np.asarray(point) - np.asarray(start_point) + vec_2 = np.asarray(end_point) - np.asarray(start_point) + + # cross product of the two vector + cross_product = np.cross(vec_1, vec_2) + # It should be smaller than tolerance for points to be aligned + return np.sum(np.abs(cross_product)) < a_tol + + @staticmethod + def is_point_between_points(point:Sequence[float], start_point:Sequence[float], end_point:Sequence[float]): + """ + Check if point is between the other two. Points must have the same dimension (2D or 3D). + + Args: + - point (Sequence[float]): point coordinates to be tested + - start_point (Sequence[float]): first extreme coordinates of the line + - end_point (Sequence[float]): second extreme coordinates of the line + + Raises: + - ValueError: when there is a dimension mismatch in the point dimensions. + + Returns: + - bool: whether the point is between the other two or not + """ + + # check dimensions of points for validation + Utils.check_dimensions([point, start_point, end_point]) + + # Calculate vectors between the points + vec_1 = np.asarray(point) - np.asarray(start_point) + vec_2 = np.asarray(end_point) - np.asarray(start_point) + + # Calculate the scalar projections of vector1 onto vector2 + scalar_projection = sum(v1 * v2 for v1, v2 in zip(vec_1, vec_2)) / sum(v ** 2 for v in vec_2) + + # Check if the scalar projection is between 0 and 1 (inclusive) + return 0 <= scalar_projection <= 1 diff --git a/stem/water_boundaries.py b/stem/water_boundaries.py new file mode 100644 index 000000000..36f035edb --- /dev/null +++ b/stem/water_boundaries.py @@ -0,0 +1,126 @@ +from typing import List, Union +from dataclasses import dataclass, field +from abc import ABC + + + +@dataclass +class WaterBoundaryParameters(ABC): + """ + Abstract base class for load water boundary parameters + + Args: + - surfaces_assigment (List[str]): List of surfaces to which the water boundary is assigned. + - is_fixed (bool): True if the water boundary is fixed, False otherwise. + - gravity_direction (int): Direction of the gravity vector. + - out_of_plane_direction (int): Direction of the out of plane vector. + + + """ + surfaces_assigment: List[str] + is_fixed: bool + gravity_direction: int + out_of_plane_direction: int + + +@dataclass +class PhreaticMultiLineBoundary(WaterBoundaryParameters): + """ + Class containing the load parameters for a phreatic line boundary condition + + Args: + - x_coordinates (List[float]): X coordinates of the phreatic line [m]. + - y_coordinates (List[float]): Y coordinates of the phreatic line [m]. + - z_coordinates (List[float]): Z coordinates of the phreatic line [m]. + - specific_weight (float): Specific weight of the water. + - water_pressure (float): Water pressure. + + + """ + x_coordinates: List[float] + y_coordinates: List[float] + specific_weight: float + water_pressure: float + z_coordinates: List[float] = field(default_factory=lambda: [0.0]) + + def __post_init__(self): + """ + Post initialization method of the class. It checks that the coordinates are of the same length. + + """ + + # Check that the coordinates are of the same length + if len(self.x_coordinates) != len(self.y_coordinates): + raise ValueError("The x and y coordinates must be of the same length") + # check if coordinate z is defined + if len(self.z_coordinates) > 1: + if len(self.x_coordinates) != len(self.z_coordinates): + raise ValueError("The x/y and z coordinates must be of the same length") + else: + # define default z coordinates + self.z_coordinates = [0.0] * len(self.x_coordinates) + + + @property + def type(self): + return "Phreatic_Multi_Line" + +@dataclass +class InterpolateLineBoundary(WaterBoundaryParameters): + """ + Class containing the boundary parameters for a interpolate line boundary condition. + + """ + + @property + def type(self): + return "Interpolate_Line" + + +@dataclass +class PhreaticLine(WaterBoundaryParameters): + """ + Class containing the boundary parameters for phreatic line boundary condition. This condition is should only contain + two points. + + Args: + - first_reference_coordinate (List[float]): First reference coordinate of the phreatic line [m]. + - second_reference_coordinate (List[float]): Second reference coordinate of the phreatic line [m]. + - specific_weight (float): Specific weight of the water . + - value (float): Value of the water pressure . + + + """ + first_reference_coordinate: List[float] + second_reference_coordinate: List[float] + specific_weight: float + value: float + + @property + def type(self): + return "Phreatic_Line" + + +class WaterBoundary: + """ + Class containing water boundary information acting on a body part + + Args: + - water_boundary (:class:`WaterBoundaryParameters`): Water boundary parameters + - type (str): Type of water boundary + + """ + + def __init__(self, water_boundary_parameters: Union[InterpolateLineBoundary, PhreaticMultiLineBoundary, PhreaticLine], name: str): + """ + Constructor of the class + + Args: + - water_boundary (:class:`WaterBoundaryParameters`): Water boundary parameters + + """ + + self.water_boundary: Union[InterpolateLineBoundary, PhreaticMultiLineBoundary, PhreaticLine] = water_boundary_parameters + self.type: str = self.water_boundary.type + self.name: str = name + diff --git a/tests/test_data/expected_geometry_after_sync_3D.pickle b/tests/test_data/expected_geometry_after_sync_3D.pickle new file mode 100644 index 000000000..5a8c99e73 Binary files /dev/null and b/tests/test_data/expected_geometry_after_sync_3D.pickle differ diff --git a/tests/test_data/expected_water_lines.json b/tests/test_data/expected_water_lines.json new file mode 100644 index 000000000..847433de9 --- /dev/null +++ b/tests/test_data/expected_water_lines.json @@ -0,0 +1,80 @@ +{ + "test": [ + { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.water_soils_1", + "variable_name": "WATER_PRESSURE", + "is_fixed": true, + "value": 0.0, + "table": [ + 0, + 0, + 0 + ], + "fluid_pressure_type": "Phreatic_Multi_Line", + "gravity_direction": 1, + "out_of_plane_direction": 2, + "x_coordinates": [ + -40.0, + -11.4, + 0.0, + 9.0, + 21.5, + 95.0 + ], + "y_coordinates": [ + 0.44, + 0.44, + 3.0, + 3.0, + -0.5, + -0.5 + ], + "z_coordinates": [ + 0, + 0, + 0, + 0, + 0, + 0 + ], + "specific_weight": 9.81 + } + }, + { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.water_soils_2", + "variable_name": "WATER_PRESSURE", + "is_fixed": true, + "table": 0, + "fluid_pressure_type": "Interpolate_Line", + "gravity_direction": 1, + "out_of_plane_direction": 2 + } + }, + { + "python_module": "apply_scalar_constraint_table_process", + "kratos_module": "KratosMultiphysics.GeoMechanicsApplication", + "process_name": "ApplyScalarConstraintTableProcess", + "Parameters": { + "model_part_name": "PorousDomain.water_soils_3", + "variable_name": "WATER_PRESSURE", + "is_fixed": true, + "value": 0.0, + "table": [0, 0], + "fluid_pressure_type": "Phreatic_Line", + "gravity_direction": 1, + "out_of_plane_direction": 2, + "first_reference_coordinate" : [0.0,1.0,0.0], + "second_reference_coordinate": [1.0,0.5,0.0], + "specific_weight": 10000.0 + } + } + ] +} \ No newline at end of file diff --git a/tests/test_data/gmsh_utils_two_blocks_2D.geo b/tests/test_data/gmsh_utils_two_blocks_2D.geo new file mode 100644 index 000000000..0a6847470 --- /dev/null +++ b/tests/test_data/gmsh_utils_two_blocks_2D.geo @@ -0,0 +1,35 @@ +// Gmsh project: created with gmsh-3.0.6-Windows64 + +// Create 2D square mesh +Mesh.ElementOrder = 1; +Point(1) = {0, 0, 0}; +Point(2) = {1, 0, 0}; +Point(3) = {1, 1, 0}; +Point(4) = {0, 1, 0}; + +// create lines +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; + +// create surface +Line Loop(1) = {1, 2, 3, 4}; +Plane Surface(1) = 1; + +// create new points of second surface +Point(5) = {0, 2, 0}; +Point(6) = {1, 2, 0}; + +// create new lines +Line(5) = {4, 5}; +Line(6) = {5, 6}; +Line(7) = {6, 3}; + +// create second surface +Line Loop(2) = {3, 5, 6, 7}; +Plane Surface(2) = 2; + +// Define the physical groups +Physical Surface("group_1") = 1; +Physical Surface("group_2") = 2; diff --git a/tests/test_default_materials.py b/tests/test_default_materials.py index 0182eabf6..eec2c5d06 100644 --- a/tests/test_default_materials.py +++ b/tests/test_default_materials.py @@ -40,5 +40,5 @@ def test_default_structural_materials(self): # compare json files using custom dictionary comparison TestUtils.assert_dictionary_almost_equal( - test_dict, expected_material_parameters_json + expected_material_parameters_json, test_dict ) diff --git a/tests/test_geometry.py b/tests/test_geometry.py index ff6906ab6..19bd6a8d0 100644 --- a/tests/test_geometry.py +++ b/tests/test_geometry.py @@ -1,5 +1,6 @@ import pytest from gmsh_utils.gmsh_IO import GmshIO +import numpy.testing as npt from stem.geometry import * @@ -12,7 +13,7 @@ def expected_geo_data_0D(self): Expected geometry data for a 0D geometry group. The group is a geometry of a point Returns: - - expected_geo_data (Dict[str, Any]): dictionary containing the geometry data as generated by the gmsh_io + - Dict[str, Any]: dictionary containing the geometry data as generated by the gmsh_io """ expected_points = {1: [0, 0, 0], 2: [0.5, 0, 0]} return {"points": expected_points} @@ -24,7 +25,7 @@ def expected_geo_data_1D(self): Expected geometry data for a 1D geometry group. The group is a geometry of a line Returns: - - expected_geo_data (Dict[str, Any]): dictionary containing the geometry data as generated by the gmsh_io + - Dict[str, Any]: dictionary containing the geometry data as generated by the gmsh_io """ expected_points = {4: [0, 1.0, 0], 11: [0, 2.0, 0], 12: [0.5, 2.0, 0]} expected_lines = {13: [4, 11], 14: [11, 12]} @@ -38,7 +39,7 @@ def expected_geo_data_2D(self): Expected geometry data for a 2D geometry group. The group is a geometry of a square. Returns: - - expected_geo_data (Dict[str, Any]): dictionary containing the geometry data as generated by the gmsh_io + - Dict[str, Any]: dictionary containing the geometry data as generated by the gmsh_io """ expected_points = {3: [0.5, 1, 0], 4: [0, 1, 0], 11: [0, 2, 0], 12: [0.5, 2.0, 0]} expected_lines = { 7: [3, 4], 13: [4, 11], 14: [11, 12], 15: [12, 3]} @@ -54,7 +55,7 @@ def expected_geo_data_3D(self): Expected geometry data for a 3D geometry group. The group is a geometry of a cubic block Returns: - - expected_geo_data (Dict[str, Any]): dictionary containing the geometry data as generated by the gmsh_io + - Dict[str, Any]: dictionary containing the geometry data as generated by the gmsh_io """ expected_points = {3: [0.5, 1., 0.], 4: [0., 1., 0.], 11: [0., 2., 0.], 12: [0.5, 2., 0.], 18: [0.5, 1., -0.5], @@ -70,12 +71,12 @@ def expected_geo_data_3D(self): "surfaces": expected_surfaces, "volumes": expected_volumes} - def test_create_0d_geometry_from_gmsh_group(self, expected_geo_data_0D): + def test_create_0d_geometry_from_gmsh_group(self, expected_geo_data_0D: Dict[str, Any]): """ Test the creation of a 0D geometry from a gmsh group. Args: - - expected_geo_data_0D (Dict[int, Any]): expected geometry data for a 0D geometry group. + - expected_geo_data_0D (Dict[str, Any]): expected geometry data for a 0D geometry group. """ # Read the gmsh geo file @@ -88,10 +89,11 @@ def test_create_0d_geometry_from_gmsh_group(self, expected_geo_data_0D): # Assert that the geometry is created correctly assert len(geometry.points) == len(expected_geo_data_0D["points"]) - for point in geometry.points: - assert pytest.approx(point.coordinates) == expected_geo_data_0D["points"][point.id] + for point_id, point in geometry.points.items(): + assert point_id == point.id + npt.assert_allclose(point.coordinates, expected_geo_data_0D["points"][point.id]) - def test_create_1d_geometry_from_gmsh_group(self, expected_geo_data_1D): + def test_create_1d_geometry_from_gmsh_group(self, expected_geo_data_1D: Dict[str, Any]): """ Test the creation of a 1D geometry from a gmsh group. @@ -110,14 +112,16 @@ def test_create_1d_geometry_from_gmsh_group(self, expected_geo_data_1D): # Assert that the geometry is created correctly assert len(geometry.points) == len(expected_geo_data_1D["points"]) - for point in geometry.points: - assert pytest.approx(point.coordinates) == expected_geo_data_1D["points"][point.id] + for point_id, point in geometry.points.items(): + assert point_id == point.id + npt.assert_allclose(point.coordinates, expected_geo_data_1D["points"][point.id]) assert len(geometry.lines) == len(expected_geo_data_1D["lines"]) - for line in geometry.lines: - assert line.point_ids == expected_geo_data_1D["lines"][line.id] + for line_id, line in geometry.lines.items(): + assert line_id == line.id + npt.assert_equal(line.point_ids, expected_geo_data_1D["lines"][line.id]) - def test_create_2d_geometry_from_gmsh_group(self, expected_geo_data_2D): + def test_create_2d_geometry_from_gmsh_group(self, expected_geo_data_2D: Dict[str, Any]): """ Test the creation of a 2D geometry from a gmsh group. @@ -136,24 +140,26 @@ def test_create_2d_geometry_from_gmsh_group(self, expected_geo_data_2D): # Assert that the geometry is created correctly assert len(geometry.points) == len(expected_geo_data_2D["points"]) - for point in geometry.points: - assert pytest.approx(point.coordinates) == expected_geo_data_2D["points"][point.id] + for point_id, point in geometry.points.items(): + assert point_id == point.id + npt.assert_allclose(point.coordinates, expected_geo_data_2D["points"][point.id]) assert len(geometry.lines) == len(expected_geo_data_2D["lines"]) - for line in geometry.lines: - assert line.point_ids == expected_geo_data_2D["lines"][line.id] + for line_id, line in geometry.lines.items(): + assert line_id == line.id + npt.assert_equal(line.point_ids, expected_geo_data_2D["lines"][line.id]) assert len(geometry.surfaces) == len(expected_geo_data_2D["surfaces"]) - for surface in geometry.surfaces: - assert surface.line_ids == expected_geo_data_2D["surfaces"][surface.id] - + for surface_id, surface in geometry.surfaces.items(): + assert surface_id == surface.id + npt.assert_equal(surface.line_ids, expected_geo_data_2D["surfaces"][surface.id]) def test_create_3d_geometry_from_gmsh_group(self, expected_geo_data_3D: Dict[str, Any]): """ Test the creation of a 3D geometry from a gmsh group. Args: - - expected_geo_data_3D (Dict[int, Any]): expected geometry data for a 3D geometry group. + - expected_geo_data_3D (Dict[str, Any]): expected geometry data for a 3D geometry group. """ @@ -167,19 +173,56 @@ def test_create_3d_geometry_from_gmsh_group(self, expected_geo_data_3D: Dict[str # Assert that the geometry is created correctly assert len(geometry.points) == len(expected_geo_data_3D["points"]) - for point in geometry.points: - assert pytest.approx(point.coordinates) == expected_geo_data_3D["points"][point.id] + for point_id, point in geometry.points.items(): + assert point_id == point.id + npt.assert_allclose(point.coordinates, expected_geo_data_3D["points"][point.id]) assert len(geometry.lines) == len(expected_geo_data_3D["lines"]) - for line in geometry.lines: - assert line.point_ids == expected_geo_data_3D["lines"][line.id] + for line_id, line in geometry.lines.items(): + assert line_id == line.id + npt.assert_equal(line.point_ids, expected_geo_data_3D["lines"][line.id]) assert len(geometry.surfaces) == len(expected_geo_data_3D["surfaces"]) - for surface in geometry.surfaces: - assert surface.line_ids == expected_geo_data_3D["surfaces"][surface.id] + for surface_id, surface in geometry.surfaces.items(): + assert surface_id == surface.id + npt.assert_equal(surface.line_ids, expected_geo_data_3D["surfaces"][surface.id]) assert len(geometry.volumes) == len(expected_geo_data_3D["volumes"]) - for volume in geometry.volumes: - assert volume.surface_ids == expected_geo_data_3D["volumes"][volume.id] + for volume_id, volume in geometry.volumes.items(): + assert volume_id == volume.id + npt.assert_equal(volume.surface_ids, expected_geo_data_3D["volumes"][volume.id]) + + def test_create_geometry_from_geo_data(self, expected_geo_data_3D: Dict[str, Any]): + """ + Test the creation of a 3D geometry from a geo_data dictionary. + + Args: + - expected_geo_data_3D (Dict[str, Any]): expected geometry data for a 3D geometry group. + """ + + geo_data = expected_geo_data_3D + + # Create the geometry from the gmsh group + geometry = Geometry().create_geometry_from_geo_data(geo_data) + # Assert that the geometry is created correctly + assert len(geometry.points) == len(expected_geo_data_3D["points"]) + for point_id, point in geometry.points.items(): + assert point_id == point.id + npt.assert_allclose(point.coordinates, expected_geo_data_3D["points"][point.id]) + + assert len(geometry.lines) == len(expected_geo_data_3D["lines"]) + for line_id, line in geometry.lines.items(): + assert line_id == line.id + npt.assert_equal(line.point_ids, expected_geo_data_3D["lines"][line.id]) + + assert len(geometry.surfaces) == len(expected_geo_data_3D["surfaces"]) + for surface_id, surface in geometry.surfaces.items(): + assert surface_id == surface.id + npt.assert_equal(surface.line_ids, expected_geo_data_3D["surfaces"][surface.id]) + + assert len(geometry.volumes) == len(expected_geo_data_3D["volumes"]) + for volume_id, volume in geometry.volumes.items(): + assert volume_id == volume.id + npt.assert_equal(volume.surface_ids, expected_geo_data_3D["volumes"][volume.id]) diff --git a/tests/test_kratos_additional_processes_io.py b/tests/test_kratos_additional_processes_io.py index 00da44a88..c532a112e 100644 --- a/tests/test_kratos_additional_processes_io.py +++ b/tests/test_kratos_additional_processes_io.py @@ -7,7 +7,7 @@ from tests.utils import TestUtils -class KratosAdditionalProcessesIO: +class TestKratosAdditionalProcessesIO: def test_create_additional_processes_dictionaries(self): """ @@ -17,7 +17,7 @@ def test_create_additional_processes_dictionaries(self): # define constraints # Absorbing boundaries - excavation_parameters = Excavation(deactivate_soil_part=True) + excavation_parameters = Excavation(deactivate_body_model_part=True) # collect the part names and parameters into a dictionary # TODO: change later when model part is implemented diff --git a/tests/test_kratos_boundaries_io.py b/tests/test_kratos_boundaries_io.py index 83b927ccd..95ccff749 100644 --- a/tests/test_kratos_boundaries_io.py +++ b/tests/test_kratos_boundaries_io.py @@ -68,5 +68,5 @@ def test_create_boundary_condition_dictionaries(self): # assert the objects to be equal TestUtils.assert_dictionary_almost_equal( - test_dictionary, expected_load_parameters_json + expected_load_parameters_json, test_dictionary ) \ No newline at end of file diff --git a/tests/test_kratos_loads_io.py b/tests/test_kratos_loads_io.py index f2d038891..1be82affb 100644 --- a/tests/test_kratos_loads_io.py +++ b/tests/test_kratos_loads_io.py @@ -38,7 +38,7 @@ def test_create_load_process_dict(self): # collect the part names and parameters into a dictionary # TODO: change later when model part is implemented - all_boundary_parameters = { + all_load_parameters = { "test_point_load": point_load_parameters, "test_line_load": line_load_parameters, "test_surface_load": surface_load_parameters, @@ -56,7 +56,7 @@ def test_create_load_process_dict(self): # TODO: when model part are implemented, generate file through kratos_io boundaries_io = KratosLoadsIO(domain="PorousDomain") - for part_name, part_parameters in all_boundary_parameters.items(): + for part_name, part_parameters in all_load_parameters.items(): _parameters = boundaries_io.create_load_dict( part_name=part_name, parameters=part_parameters ) @@ -69,5 +69,5 @@ def test_create_load_process_dict(self): # assert the objects to be equal TestUtils.assert_dictionary_almost_equal( - test_dictionary, expected_load_parameters_json + expected_load_parameters_json, test_dictionary ) diff --git a/tests/test_kratos_material_io.py b/tests/test_kratos_material_io.py index 4697c80ca..442d150ea 100644 --- a/tests/test_kratos_material_io.py +++ b/tests/test_kratos_material_io.py @@ -142,7 +142,7 @@ def test_write_soil_material_dict(self): # compare json files using custom dictionary comparison TestUtils.assert_dictionary_almost_equal( - test_dict, expected_material_parameters_json + expected_material_parameters_json, test_dict ) def test_write_structural_material_dict(self): @@ -216,5 +216,5 @@ def test_write_structural_material_dict(self): # compare json files using custom dictionary comparison TestUtils.assert_dictionary_almost_equal( - test_dict, expected_material_parameters_json + expected_material_parameters_json, test_dict ) diff --git a/tests/test_kratos_outputs_io.py b/tests/test_kratos_outputs_io.py index bf238651b..cd126a4af 100644 --- a/tests/test_kratos_outputs_io.py +++ b/tests/test_kratos_outputs_io.py @@ -131,5 +131,5 @@ def test_create_output_process_dictionary(self): # assert the objects to be equal TestUtils.assert_dictionary_almost_equal( - test_output, expected_load_parameters_json + expected_load_parameters_json, test_output ) diff --git a/tests/test_kratos_solver_io.py b/tests/test_kratos_solver_io.py index e09872473..7a2144638 100644 --- a/tests/test_kratos_solver_io.py +++ b/tests/test_kratos_solver_io.py @@ -48,11 +48,9 @@ def test_create_settings_dictionary(self): problem_data = Problem(problem_name="test", number_of_threads=2, settings=solver_settings) # create model parts - model_part1 = ModelPart() - model_part1.name = "ModelPart1" + model_part1 = ModelPart("ModelPart1") - body_model_part1 = BodyModelPart() - body_model_part1.name = "BodyModelPart1" + body_model_part1 = BodyModelPart("BodyModelPart1") model_parts = [model_part1, body_model_part1] diff --git a/tests/test_kratos_water_boundaries_io.py b/tests/test_kratos_water_boundaries_io.py new file mode 100644 index 000000000..5a6870e8e --- /dev/null +++ b/tests/test_kratos_water_boundaries_io.py @@ -0,0 +1,63 @@ +from tests.utils import TestUtils +import json + +from stem.IO.kratos_water_boundaries_io import KratosWaterBoundariesIO +from stem.water_boundaries import WaterBoundary, InterpolateLineBoundary, PhreaticMultiLineBoundary, PhreaticLine + + +class TestKratosWaterBoundariesIO: + + def test_create_water_boundary_process_dict(self): + """ + + Test the creation of the water boundary process dictionary for the + ProjectParameters.json file + + """ + multi_line_boundary = PhreaticMultiLineBoundary( + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + water_pressure=0, + x_coordinates=[-40.0, -11.4, 0.0, 9.0, 21.5, 95.0], + y_coordinates=[0.44, 0.44, 3.0, 3.0, -0.5, -0.5], + surfaces_assigment=["domain a", "domain b", "domain c"], + specific_weight=9.81, + ) + water_boundary = WaterBoundary(multi_line_boundary, name="water_soils_1") + # use the kratos io to create the dictionary + kratos_io = KratosWaterBoundariesIO(domain="PorousDomain") + # set the interpolation type + interpolation_type = InterpolateLineBoundary( + surfaces_assigment=["domain d"], + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + ) + water_boundary_interpolate = WaterBoundary(interpolation_type, name="water_soils_2") + # check phreatic line + phreatic_line = PhreaticLine( + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + value=0, + first_reference_coordinate=[0.0,1.0,0.0], + second_reference_coordinate=[1.0,0.5,0.0], + specific_weight=10000.0, + surfaces_assigment=["domain e"] + ) + water_boundary_phreatic_line = WaterBoundary(phreatic_line, name="water_soils_3") + + # check the dictionary + # read the expected dictionary from the json + with open("tests/test_data/expected_water_lines.json") as json_file: + expected_water_boundary_json = json.load(json_file) + # compare the dictionaries + TestUtils.assert_dictionary_almost_equal(expected_water_boundary_json['test'][0], + kratos_io.create_water_boundary_dict( + water_boundary + )) + TestUtils.assert_dictionary_almost_equal(expected_water_boundary_json['test'][1], + kratos_io.create_water_boundary_dict(water_boundary_interpolate)) + TestUtils.assert_dictionary_almost_equal(expected_water_boundary_json['test'][2], + kratos_io.create_water_boundary_dict(water_boundary_phreatic_line)) diff --git a/tests/test_mesh.py b/tests/test_mesh.py new file mode 100644 index 000000000..fd0f98fed --- /dev/null +++ b/tests/test_mesh.py @@ -0,0 +1,229 @@ +import pytest + +from stem.mesh import * + + +class TestMesh: + + def test_create_0d_mesh_from_gmsh_group(self): + """ + Test the creation of a 0D mesh from a gmsh group. + + """ + + # Set up the mesh data + mesh_data = {"ndim": 0, + "nodes": {1: [0, 0, 0], 2: [0.5, 0, 0]}, + "elements": {"POINT_1N": {1: [1], 2: [2]}}, + "physical_groups": {"points_group": {"ndim": 0, + 'element_ids': [1, 2], + "node_ids": [1, 2], + "element_type": "POINT_1N"}}} + + # Create the mesh from the gmsh group + generated_mesh = Mesh.create_mesh_from_gmsh_group(mesh_data, "points_group") + + # set expected mesh + expected_nodes = [Node(1, [0, 0, 0]), Node(2, [0.5, 0, 0])] + expected_elements = [Element(1, "POINT_1N", [1]), Element(2, "POINT_1N", [2])] + expected_mesh = Mesh(0) + expected_mesh.nodes = expected_nodes + expected_mesh.elements = expected_elements + + # Check the generated mesh + assert generated_mesh.ndim == expected_mesh.ndim + + # Check the nodes + for generated_node, expected_node in zip(generated_mesh.nodes, expected_mesh.nodes): + assert generated_node.id == expected_node.id + assert pytest.approx(generated_node.coordinates) == expected_node.coordinates + + # Check the elements + for generated_element, expected_element in zip(generated_mesh.elements, expected_mesh.elements): + assert generated_element.id == expected_element.id + assert generated_element.element_type == expected_element.element_type + assert generated_element.node_ids == expected_element.node_ids + + def test_create_1d_mesh_from_gmsh_group(self): + """ + Test the creation of a 1D mesh from a gmsh group. + + """ + + # Set up the mesh data + mesh_data = {"ndim": 1, + "nodes": {1: [0, 0, 0], 2: [0.5, 0, 0], 3: [1, 0, 0]}, + "elements": {"LINE_2N": {1: [1, 2], 2: [2, 3]}}, + "physical_groups": {"lines_group": {"ndim": 1, + 'element_ids': [1, 2], + "node_ids": [1, 2, 3], + "element_type": "LINE_2N"}}} + + # Create the mesh from the gmsh group + generated_mesh = Mesh.create_mesh_from_gmsh_group(mesh_data, "lines_group") + + # set expected mesh + expected_nodes = [Node(1, [0, 0, 0]), Node(2, [0.5, 0, 0]), Node(3, [1, 0, 0])] + expected_elements = [Element(1, "LINE_2N", [1, 2]), Element(2, "LINE_2N", [2, 3])] + expected_mesh = Mesh(1) + expected_mesh.nodes = expected_nodes + expected_mesh.elements = expected_elements + + # Check the generated mesh + assert generated_mesh.ndim == expected_mesh.ndim + + # Check the nodes + for generated_node, expected_node in zip(generated_mesh.nodes, expected_mesh.nodes): + assert generated_node.id == expected_node.id + assert pytest.approx(generated_node.coordinates) == expected_node.coordinates + + # Check the elements + for generated_element, expected_element in zip(generated_mesh.elements, expected_mesh.elements): + assert generated_element.id == expected_element.id + assert generated_element.element_type == expected_element.element_type + assert generated_element.node_ids == expected_element.node_ids + + def test_create_2d_mesh_from_gmsh_group(self): + """ + Test the creation of a 2D mesh from a gmsh group. + + """ + + # Set up the mesh data + mesh_data = {"ndim": 2, + "nodes": {1: [0, 0, 0], 2: [0.5, 0, 0], 3: [1, 0, 0], + 4: [0, 0.5, 0], 5: [0.5, 0.5, 0], 6: [1, 0.5, 0]}, + "elements": {"TRIANGLE_3N": {1: [1, 2, 4], 2: [2, 3, 5], 3: [3, 6, 5]}}, + "physical_groups": {"triangles_group": {"ndim": 2, + 'element_ids': [1, 2, 3], + "node_ids": [1, 2, 3, 4, 5, 6], + "element_type": "TRIANGLE_3N"}}} + + # Create the mesh from the gmsh group + generated_mesh = Mesh.create_mesh_from_gmsh_group(mesh_data, "triangles_group") + + # set expected mesh + expected_nodes = [Node(1, [0, 0, 0]), Node(2, [0.5, 0, 0]), Node(3, [1, 0, 0]), + Node(4, [0, 0.5, 0]), Node(5, [0.5, 0.5, 0]), Node(6, [1, 0.5, 0])] + expected_elements = [Element(1, "TRIANGLE_3N", [1, 2, 4]), + Element(2, "TRIANGLE_3N", [2, 3, 5]), + Element(3, "TRIANGLE_3N", [3, 6, 5])] + + expected_mesh = Mesh(2) + expected_mesh.nodes = expected_nodes + expected_mesh.elements = expected_elements + + # Check the generated mesh + assert generated_mesh.ndim == expected_mesh.ndim + + # Check the nodes + for generated_node, expected_node in zip(generated_mesh.nodes, expected_mesh.nodes): + assert generated_node.id == expected_node.id + assert pytest.approx(generated_node.coordinates) == expected_node.coordinates + + # Check the elements + for generated_element, expected_element in zip(generated_mesh.elements, expected_mesh.elements): + assert generated_element.id == expected_element.id + assert generated_element.element_type == expected_element.element_type + assert generated_element.node_ids == expected_element.node_ids + + def test_create_3d_mesh_from_gmsh_group(self): + """ + Test the creation of a 3D mesh from a gmsh group. + + """ + + # Set up the mesh data + mesh_data = {"ndim": 3, + "nodes": {1: [0, 0, 0], 2: [0.5, 0, 0], 3: [1, 0, 0], + 4: [0, 0.5, 0], 5: [0.5, 0.5, 0], 6: [1, 0.5, 0], + 7: [0, 0, 0.5], 8: [0.5, 0, 0.5], 9: [1, 0, 0.5], + 10: [0, 0.5, 0.5], 11: [0.5, 0.5, 0.5], 12: [1, 0.5, 0.5]}, + "elements": {"TETRAHEDRON_4N": {1: [1, 2, 4, 7], 2: [2, 3, 5, 8], 3: [3, 6, 5, 9], + 4: [4, 5, 7, 10], 5: [5, 6, 8, 11], 6: [6, 9, 11, 8]}}, + "physical_groups": {"tetrahedral_group": {"ndim": 3, + 'element_ids': [1, 2, 3, 4, 5, 6], + "node_ids": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], + "element_type": "TETRAHEDRON_4N"}}} + + # Create the mesh from the gmsh group + generated_mesh = Mesh.create_mesh_from_gmsh_group(mesh_data, "tetrahedral_group") + + # set expected mesh + expected_nodes = [Node(1, [0, 0, 0]), Node(2, [0.5, 0, 0]), Node(3, [1, 0, 0]), + Node(4, [0, 0.5, 0]), Node(5, [0.5, 0.5, 0]), Node(6, [1, 0.5, 0]), + Node(7, [0, 0, 0.5]), Node(8, [0.5, 0, 0.5]), Node(9, [1, 0, 0.5]), + Node(10, [0, 0.5, 0.5]), Node(11, [0.5, 0.5, 0.5]), Node(12, [1, 0.5, 0.5])] + + expected_elements = [Element(1, "TETRAHEDRON_4N", [1, 2, 4, 7]), + Element(2, "TETRAHEDRON_4N", [2, 3, 5, 8]), + Element(3, "TETRAHEDRON_4N", [3, 6, 5, 9]), + Element(4, "TETRAHEDRON_4N", [4, 5, 7, 10]), + Element(5, "TETRAHEDRON_4N", [5, 6, 8, 11]), + Element(6, "TETRAHEDRON_4N", [6, 9, 11, 8])] + + expected_mesh = Mesh(3) + expected_mesh.nodes = expected_nodes + expected_mesh.elements = expected_elements + + # Check the generated mesh + assert generated_mesh.ndim == expected_mesh.ndim + + # Check the nodes + for generated_node, expected_node in zip(generated_mesh.nodes, expected_mesh.nodes): + assert generated_node.id == expected_node.id + assert pytest.approx(generated_node.coordinates) == expected_node.coordinates + + # Check the elements + for generated_element, expected_element in zip(generated_mesh.elements, expected_mesh.elements): + assert generated_element.id == expected_element.id + assert generated_element.element_type == expected_element.element_type + assert generated_element.node_ids == expected_element.node_ids + + def test_create_mesh_from_non_existing_group(self): + """ + Test the creation of a mesh from a non-existing gmsh group. Expected to raise a ValueError. + + """ + + # Set up the mesh data + mesh_data = {"ndim": 0, + "nodes": {1: [0, 0, 0], 2: [0.5, 0, 0]}, + "elements": {"POINT_1N": {1: [1], 2: [2]}}, + "physical_groups": {"points_group": {"ndim": 0, + 'element_ids': [1, 2], + "node_ids": [1, 2], + "element_type": "POINT_1N"}}} + + # Create the mesh from the gmsh group + with pytest.raises(ValueError): + Mesh.create_mesh_from_gmsh_group(mesh_data, "non_existing_group") + + +class TestMeshSettings: + """ + Test the mesh settings class. + """ + + def test_validation_element_order_at_initialisation_expected_raise(self): + """ + Test the validation of the element order. Expected to raise a ValueError when the element order is not 1 or 2. + + """ + + # test if ValueError is raised when element_order is not 1 or 2 + with pytest.raises(ValueError): + + MeshSettings(element_order=3) + + def test_validation_element_order_after_initialisation_expected_raise(self): + """ + Test the validation of the element order. Expected to raise a ValueError when the element order is not 1 or 2. + + """ + + # test if ValueError is raised when element_order is not 1 or 2 + mesh_settings = MeshSettings() + + with pytest.raises(ValueError): + mesh_settings.element_order = 3 \ No newline at end of file diff --git a/tests/test_model.py b/tests/test_model.py new file mode 100644 index 000000000..dde45d386 --- /dev/null +++ b/tests/test_model.py @@ -0,0 +1,1646 @@ +from typing import Tuple +import pickle + +import pytest +from gmsh_utils import gmsh_IO +import numpy.testing as npt + +from stem.model import * +from stem.geometry import * +from tests.utils import TestUtils +from stem.solver import * +from stem.boundary import * + + +class TestModel: + + @pytest.fixture + def expected_geo_data_0D(self): + """ + Expected geometry data for a 0D geometry group. The group is a geometry of a point + + Returns: + - Dict[str, Any]: dictionary containing the geometry data as generated by the gmsh_io + """ + expected_points = {1: [0, 0, 0], 2: [0.5, 0, 0]} + return {"points": expected_points} + + @pytest.fixture + def expected_geometry_single_layer_2D(self): + """ + Sets expected geometry data for a 2D geometry group. The group is a geometry of a square. + + Returns: + - :class:`stem.geometry.Geometry`: geometry of a 2D square + """ + + geometry = Geometry() + + geometry.points = {1: Point.create([0, 0, 0], 1), + 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3), + 4: Point.create([0, 1, 0], 4)} + + geometry.lines = {1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3), + 4: Line.create([4, 1], 4)} + + geometry.surfaces = {1: Surface.create([1,2,3,4], 1)} + + geometry.volumes = {} + + return geometry + + + @pytest.fixture + def expected_geometry_single_layer_3D(self): + """ + Sets expected geometry data for a 3D geometry group. The group is a geometry of a cube. + + Returns: + - :class:`stem.geometry.Geometry`: geometry of a 3D cube + """ + + geometry = Geometry() + + geometry.points = {1, Point.create([0, 0, 0], 1), + 5, Point.create([0, 0, 1], 5), + 6, Point.create([1, 0, 1], 6), + 2, Point.create([1, 0, 0], 2), + 7, Point.create([1, 1, 1], 7), + 3, Point.create([1, 1, 0], 3), + 8, Point.create([0, 1, 1], 8), + 4, Point.create([0, 1, 0], 4)} + + geometry.lines = {5, Line.create([1, 5], 5), + 7, Line.create([5, 6], 7), + 6, Line.create([2, 6], 6), + 1, Line.create([1, 2], 1), + 9, Line.create([6, 7], 9), + 8, Line.create([3, 7], 8), + 2, Line.create([2, 3], 2), + 11, Line.create([7, 8], 11), + 10, Line.create([4, 8], 10), + 3, Line.create([3, 4], 3), + 12, Line.create([8, 5], 12), + 4, Line.create([4, 1], 4)} + + geometry.surfaces = {2, Surface.create([5, 7, -6, -1], 2), + 3, Surface.create([6, 9, -8, -2], 3), + 4, Surface.create([8,11, -10, -3], 4), + 5, Surface.create([10, 12, -5, -4], 5), + 1, Surface.create([1, 2, 3, 4], 1), + 6, Surface.create([7, 9, 11, 12], 6)} + + geometry.volumes = {1, Volume.create([-2, -3, -4, -5, -1, 6], 1)} + + return geometry + + @pytest.fixture + def expected_geometry_single_layer_3D(self): + """ + Sets expected geometry data for a 2D geometry group. The group is a geometry of a square. + + Returns: + - :class:`stem.geometry.Geometry`: geometry of a 2D square + """ + + geometry = Geometry() + + geometry.points = {1: Point.create([0, 0, 0], 1), + 5: Point.create([0, 0, 1], 5), + 6: Point.create([1, 0, 1], 6), + 2: Point.create([1, 0, 0], 2), + 7: Point.create([1, 1, 1], 7), + 3: Point.create([1, 1, 0], 3), + 8: Point.create([0, 1, 1], 8), + 4: Point.create([0, 1, 0], 4)} + + geometry.lines = {5: Line.create([1, 5], 5), + 7: Line.create([5, 6], 7), + 6: Line.create([2, 6], 6), + 1: Line.create([1, 2], 1), + 9: Line.create([6, 7], 9), + 8: Line.create([3, 7], 8), + 2: Line.create([2, 3], 2), + 11: Line.create([7, 8], 11), + 10: Line.create([4, 8], 10), + 3: Line.create([3, 4], 3), + 12: Line.create([8, 5], 12), + 4: Line.create([4, 1], 4)} + + geometry.surfaces = {2: Surface.create([5, 7, -6, -1], 2), + 3: Surface.create([6, 9, -8, -2], 3), + 4: Surface.create([8, 11, -10, -3], 4), + 5: Surface.create([10, 12, -5, -4], 5), + 1: Surface.create([1, 2, 3, 4], 1), + 6: Surface.create([7, 9, 11, 12], 6)} + + # The volumes list converted to a dictionary + geometry.volumes = {1: Volume.create([-2, -3, -4, -5, -1, 6], 1)} + + + return geometry + + @pytest.fixture + def expected_geometry_two_layers_2D(self): + """ + Sets expected geometries for 2 attached 2D squares. + + Returns: + - Tuple[:class:`stem.geometry.Geometry`,:class:`stem.geometry.Geometry`]: \ + geometries of 2 attached 2D squares + + """ + + # geometry_1 + geometry_1 = Geometry() + geometry_1.points = {1: Point.create([0, 0, 0], 1), + 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3), + 4: Point.create([0, 1, 0], 4)} + + geometry_1.lines = {1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3), + 4: Line.create([4, 1], 4)} + + geometry_1.surfaces = {1: Surface.create([1, 2, 3, 4], 1)} + + geometry_1.volumes = {} + + # geometry_2 + geometry_2 = Geometry() + + geometry_2.points = {5: Point.create([1, 2, 0], 5), + 6: Point.create([0, 2, 0], 6), + 4: Point.create([0, 1, 0], 4), + 3: Point.create([1, 1, 0], 3)} + + geometry_2.lines = {5: Line.create([5, 6],5), + 6: Line.create([6, 4], 6), + 3: Line.create([3, 4], 3), + 7: Line.create([3, 5], 7)} + + geometry_2.surfaces = {2: Surface.create([5, 6, -3, 7], 2)} + + geometry_2.volumes = {} + + return geometry_1, geometry_2 + + @pytest.fixture + def expected_geometry_two_layers_2D_after_sync(self): + """ + Sets expected geometry of two model parts and the whole model after synchronising the geometry. + + Returns: + - Tuple[:class:`stem.geometry.Geometry`,:class:`stem.geometry.Geometry`, \ + :class:`stem.geometry.Geometry`]: geometries of 2 attached 2D squares and the whole model + """ + + # create expected geometry layer 1 + geometry_1 = Geometry() + + geometry_1.points = { + 1: Point.create([0, 0, 0], 1), + 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3), + 4: Point.create([0.5, 1, 0], 4), + 5: Point.create([0, 1, 0], 5) + } + + geometry_1.lines = { + 1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3), + 4: Line.create([4, 5], 4), + 5: Line.create([5, 1], 5) + } + + geometry_1.surfaces = { + 1: Surface.create([1, 2, 3, 4, 5], 1) + } + + geometry_2 = Geometry() + geometry_2.points = { + 6: Point.create([1.0, 2.0, 0.0], 6), + 7: Point.create([0.5, 2.0, 0.0], 7), + 4: Point.create([0.5, 1, 0], 4), + 3: Point.create([1, 1, 0], 3) + } + + geometry_2.lines = { + 6: Line.create([6, 7], 6), + 7: Line.create([7, 4], 7), + 3: Line.create([3, 4], 3), + 8: Line.create([3, 6], 8) + } + + geometry_2.surfaces = { + 2: Surface.create([6, 7, -3, 8], 2) + } + + geometry_2.volumes = {} + + # create expected full geometry + full_geometry = Geometry() + full_geometry.points = { + 1: Point.create([0, 0, 0], 1), + 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3), + 4: Point.create([0.5, 1, 0], 4), + 5: Point.create([0, 1, 0], 5), + 6: Point.create([1, 2, 0], 6), + 7: Point.create([0.5, 2, 0], 7) + } + + full_geometry.lines = { + 1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3), + 4: Line.create([4, 5], 4), + 5: Line.create([5, 1], 5), + 6: Line.create([6, 7], 6), + 7: Line.create([7, 4], 7), + 8: Line.create([3, 6], 8) + } + + full_geometry.surfaces = { + 1: Surface.create([1, 2, 3, 4, 5], 1), + 2: Surface.create([6, 7, -3, 8], 2) + } + + full_geometry.volumes = {} + + return geometry_1, geometry_2, full_geometry + + @pytest.fixture + def expected_geometry_line_load(self): + """ + Sets expected geometry data for a 1D geometry group. The group is a geometry of a multi-line. + + Returns: + - :class:`stem.geometry.Geometry`: geometry of a 1D multi-line + """ + + geometry = Geometry() + + geometry.points = { + 1: Point.create([0, 0, 0], 1), + 2: Point.create([3, 0, 0], 2), + 3: Point.create([4, -1, 0], 3), + 4: Point.create([10, -1, 0], 4) + } + + geometry.lines = { + 1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3) + } + + geometry.surfaces = {} + + geometry.volumes = {} + + return geometry + + @pytest.fixture + def create_default_2d_soil_material(self): + """ + Create a default soil material for a 2D geometry. + + Returns: + - :class:`stem.soil_material.SoilMaterial`: default soil material + + """ + # define soil material + ndim = 2 + soil_formulation = OnePhaseSoil(ndim, IS_DRAINED=True, DENSITY_SOLID=2650, POROSITY=0.3) + constitutive_law = LinearElasticSoil(YOUNG_MODULUS=100e6, POISSON_RATIO=0.3) + soil_material = SoilMaterial(name="soil", soil_formulation=soil_formulation, constitutive_law=constitutive_law, + retention_parameters=SaturatedBelowPhreaticLevelLaw()) + return soil_material + + @pytest.fixture + def create_default_3d_soil_material(self): + """ + Create a default soil material for a 3D geometry. + + Returns: + - :class:`stem.soil_material.SoilMaterial`: default soil material + + """ + # define soil material + ndim = 3 + soil_formulation = OnePhaseSoil(ndim, IS_DRAINED=True, DENSITY_SOLID=2650, POROSITY=0.3) + constitutive_law = LinearElasticSoil(YOUNG_MODULUS=100e6, POISSON_RATIO=0.3) + soil_material = SoilMaterial(name="soil", soil_formulation=soil_formulation, constitutive_law=constitutive_law, + retention_parameters=SaturatedBelowPhreaticLevelLaw()) + return soil_material + + @pytest.fixture + def create_default_point_load_parameters(self): + """ + Create a default point load parameters. + + Returns: + - :class:`stem.load.PointLoad`: default point load + + """ + # define soil material + return PointLoad(active=[False, True, False], value=[0, -200, 0]) + + @pytest.fixture + def create_default_line_load_parameters(self): + """ + Create a default line load parameters. + + Returns: + - :class:`stem.load.PointLoad`: default point load + + """ + # define soil material + return LineLoad(active=[False, True, False], value=[0, -20, 0]) + + @pytest.fixture + def create_default_surface_load_parameters(self): + """ + Create a default surface load properties. + + Returns: + - :class:`stem.load.SurfaceLoad`: default surface load + + """ + # define soil material + return SurfaceLoad(active=[False, True, False], value=[0, -2, 0]) + + @pytest.fixture + def create_default_moving_load_parameters(self): + """ + Create a default surface load properties. + + Returns: + - :class:`stem.load.SurfaceLoad`: default surface load + + """ + # define soil material + return MovingLoad( + origin=[3.5, -0.5, 0.0], + load=[0.0, -10.0, 0.0], + velocity=5.0, + offset=3.0, + direction=[1, 1, 1] + ) + + @pytest.fixture + def expected_geometry_two_layers_3D_extruded(self): + """ + Expected geometry data for a 3D geometry create from 2D extrusion. The geometry is 2 stacked blocks, where the + top and bottom blocks are in different groups. + + Returns: + - Tuple[:class:`stem.geometry.Geometry`,:class:`stem.geometry.Geometry`]: expected geometry data + """ + + geometry_1 = Geometry() + geometry_1.points = { + 1: Point.create([0, 0, 0], 1), + 2: Point.create([0, 0, 1], 2), + 4: Point.create([1, 0, 1], 4), + 3: Point.create([1, 0, 0], 3), + 6: Point.create([1, 1, 1], 6), + 5: Point.create([1, 1, 0], 5), + 8: Point.create([0, 1, 1], 8), + 7: Point.create([0, 1, 0], 7) + } + + geometry_1.lines = { + 1: Line.create([1, 2], 1), + 4: Line.create([2, 4], 4), + 2: Line.create([3, 4], 2), + 3: Line.create([1, 3], 3), + 7: Line.create([4, 6], 7), + 5: Line.create([5, 6], 5), + 6: Line.create([3, 5], 6), + 10: Line.create([6, 8], 10), + 8: Line.create([7, 8], 8), + 9: Line.create([5, 7], 9), + 12: Line.create([8, 2], 12), + 11: Line.create([7, 1], 11) + } + + geometry_1.surfaces = { + 1: Surface.create([1, 4, -2, -3], 1), + 2: Surface.create([2, 7, -5, -6], 2), + 3: Surface.create([5, 10, -8, -9], 3), + 4: Surface.create([8, 12, -1, -11], 4), + 5: Surface.create([3, 6, 9, 11], 5), + 6: Surface.create([4, 7, 10, 12], 6) + } + + geometry_1.volumes = { + 1: Volume.create([-1, -2, -3, -4, -5, 6], 1) + } + + geometry_2 = Geometry() + + geometry_2.points = { + 9: Point.create([1.0, 2.0, 0.0], 9), + 10: Point.create([1., 2., 1.], 10), + 12: Point.create([0.0, 2., 1.], 12), + 11: Point.create([0, 2., 0.], 11), + 8: Point.create([0., 1., 1], 8), + 7: Point.create([0., 1., 0], 7), + 5: Point.create([1, 1., 0], 5), + 6: Point.create([1, 1., 1], 6) + } + + geometry_2.lines = { + 13: Line.create([9, 10], 13), + 16: Line.create([10, 12], 16), + 14: Line.create([11, 12], 14), + 15: Line.create([9, 11], 15), + 18: Line.create([12, 8], 18), + 8: Line.create([7, 8], 8), + 17: Line.create([11, 7], 17), + 5: Line.create([5, 6], 5), + 10: Line.create([6, 8], 10), + 9: Line.create([5, 7], 9), + 20: Line.create([6, 10], 20), + 19: Line.create([5, 9], 19) + } + + geometry_2.surfaces = { + 7: Surface.create([13, 16, -14, -15], 7), + 8: Surface.create([14, 18, -8, -17], 8), + 3: Surface.create([5, 10, -8, -9], 3), + 9: Surface.create([5, 20, -13, -19], 9), + 10: Surface.create([15, 17, -9, 19], 10), + 11: Surface.create([16, 18, -10, 20], 11) + } + + geometry_2.volumes = { + 2: Volume.create([-7, -8, 3, -9, -10, 11], 2) + } + + return geometry_1, geometry_2 + + @pytest.fixture + def expected_geometry_two_layers_3D_geo_file(self): + """ + Expected geometry data for a 3D geometry create in a geo file. The geometry is 2 stacked blocks, where the top + and bottom blocks are in different groups. + + Returns: + - Tuple[:class:`stem.geometry.Geometry`,:class:`stem.geometry.Geometry`]: expected geometry data + """ + + geometry_1 = Geometry() + geometry_1.volumes = { + 1: Volume.create([-10, 39, 26, 30, 34, 38], 1) + } + + geometry_1.surfaces = { + 10: Surface.create([5, 6, 7, 8], 10), + 39: Surface.create([19, 20, 21, 22], 39), + 26: Surface.create([5, 25, -19, -24], 26), + 30: Surface.create([6, 29, -20, -25], 30), + 34: Surface.create([7, 33, -21, -29], 34), + 38: Surface.create([8, 24, -22, -33], 38) + } + + geometry_1.lines = { + 5: Line.create([1, 2], 5), + 6: Line.create([2, 3], 6), + 7: Line.create([3, 4], 7), + 8: Line.create([4, 1], 8), + 19: Line.create([13, 14], 19), + 20: Line.create([14, 18], 20), + 21: Line.create([18, 22], 21), + 22: Line.create([22, 13], 22), + 25: Line.create([2, 14], 25), + 24: Line.create([1, 13], 24), + 29: Line.create([3, 18], 29), + 33: Line.create([4, 22], 33) + } + + geometry_1.points = { + 1: Point.create([0., 0., 0.], 1), + 2: Point.create([0.5, 0., 0.], 2), + 3: Point.create([0.5, 1., 0.], 3), + 4: Point.create([0., 1., 0.], 4), + 13: Point.create([0., 0., -0.5], 13), + 14: Point.create([0.5, 0., -0.5], 14), + 18: Point.create([0.5, 1., -0.5], 18), + 22: Point.create([0., 1., -0.5], 22) + } + + geometry_2 = Geometry() + geometry_2.volumes = { + 2: Volume.create([-17, 61, -48, -34, -56, -60], 2) + } + + geometry_2.surfaces = { + 17: Surface.create([-13, -7, -15, -14], 17), + 61: Surface.create([41, -21, 43, 44], 61), + 48: Surface.create([-13, 33, -41, -46], 48), + 34: Surface.create([7, 33, -21, -29], 34), + 56: Surface.create([-15, 55, -43, -29], 56), + 60: Surface.create([-14, 46, -44, -55], 60) + } + + geometry_2.lines = { + 13: Line.create([4, 11], 13), + 7: Line.create([3, 4], 7), + 15: Line.create([12, 3], 15), + 14: Line.create([11, 12], 14), + 41: Line.create([23, 22], 41), + 21: Line.create([18, 22], 21), + 43: Line.create([18, 32], 43), + 44: Line.create([32, 23], 44), + 33: Line.create([4, 22], 33), + 46: Line.create([11, 23], 46), + 29: Line.create([3, 18], 29), + 55: Line.create([12, 32], 55) + } + + geometry_2.points = { + 4: Point.create([0., 1., 0.], 4), + 11: Point.create([0., 2., 0.], 11), + 3: Point.create([0.5, 1., 0.], 3), + 12: Point.create([0.5, 2., 0.], 12), + 23: Point.create([0., 2., -0.5], 23), + 22: Point.create([0., 1., -0.5], 22), + 18: Point.create([0.5, 1., -0.5], 18), + 32: Point.create([0.5, 2., -0.5], 32) + } + + return geometry_1, geometry_2 + + @pytest.fixture(autouse=True) + def close_gmsh(self): + """ + Initializer to close gmsh if it was not closed before. In case a test fails, the destroyer method is not called + on the Model object and gmsh keeps on running. Therefore, nodes, lines, surfaces and volumes ids are not + reset to one. This causes also the next test after the failed one to fail as well, which has nothing to do + the test itself. + + Returns: + - None + + """ + gmsh_IO.GmshIO().finalize_gmsh() + + def test_add_single_soil_layer_2D(self, expected_geometry_single_layer_2D: Geometry, + create_default_2d_soil_material: SoilMaterial): + """ + Test if a single soil layer is added correctly to the model in a 2D space. A single soil layer is generated + and a single soil material is created and added to the model. + + Args: + - expected_geometry_single_layer_2D (:class:`stem.geometry.Geometry`): expected geometry of the model + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): default soil material + + """ + + ndim = 2 + + layer_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + + # define soil material + soil_material = create_default_2d_soil_material + + # create model + model = Model(ndim) + + # add soil layer + model.add_soil_layer_by_coordinates(layer_coordinates, soil_material, "soil1") + + # check if layer is added correctly + assert len(model.body_model_parts) == 1 + assert model.body_model_parts[0].name == "soil1" + assert model.body_model_parts[0].material == soil_material + + # check if geometry is added correctly + generated_geometry = model.body_model_parts[0].geometry + expected_geometry = expected_geometry_single_layer_2D + + # check if points are added correctly + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_single_soil_layer_3D(self, expected_geometry_single_layer_3D: Geometry, + create_default_3d_soil_material: SoilMaterial): + """ + Test if a single soil layer is added correctly to the model in a 3D space. A single soil layer is generated + and a single soil material is created and added to the model. + + Args: + - expected_geometry_single_layer_3D (:class:`stem.geometry.Geometry`): expected geometry of the model + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): default soil material + + """ + + ndim = 3 + + layer_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + + # define soil material + soil_material = create_default_3d_soil_material + + # create model + model = Model(ndim) + model.extrusion_length = [0, 0, 1] + + # add soil layer + model.add_soil_layer_by_coordinates(layer_coordinates, soil_material, "soil1") + + # check if layer is added correctly + assert len(model.body_model_parts) == 1 + assert model.body_model_parts[0].name == "soil1" + assert model.body_model_parts[0].material == soil_material + + # check if geometry is added correctly + generated_geometry = model.body_model_parts[0].geometry + expected_geometry = expected_geometry_single_layer_3D + + # check if points are added correctly + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_multiple_soil_layers_2D(self, expected_geometry_two_layers_2D: Tuple[Geometry, Geometry], + create_default_2d_soil_material: SoilMaterial): + """ + Test if multiple soil layers are added correctly to the model in a 2D space. Multiple soil layers are generated + and multiple soil materials are created and added to the model. + + Args: + - expected_geometry_two_layers_2D (Tuple[:class:`stem.geometry.Geometry`, :class:`stem.geometry.Geometry`]): \ + expected geometry of the model + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): default soil material + + """ + + ndim = 2 + + layer1_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + layer2_coordinates = [(1, 1, 0), (0, 1, 0), (0, 2, 0), (1, 2, 0)] + + # define soil materials + soil_material1 = create_default_2d_soil_material + soil_material1.name = "soil1" + + soil_material2 = create_default_2d_soil_material + soil_material2.name = "soil2" + + # create model + model = Model(ndim) + + # add soil layers + model.add_soil_layer_by_coordinates(layer1_coordinates, soil_material1, "layer1") + model.add_soil_layer_by_coordinates(layer2_coordinates, soil_material2, "layer2") + + # check if layers are added correctly + assert len(model.body_model_parts) == 2 + assert model.body_model_parts[0].name == "layer1" + assert model.body_model_parts[0].material == soil_material1 + assert model.body_model_parts[1].name == "layer2" + assert model.body_model_parts[1].material == soil_material2 + + # check if geometry is added correctly for each layer + for i in range(len(model.body_model_parts)): + generated_geometry = model.body_model_parts[i].geometry + expected_geometry = expected_geometry_two_layers_2D[i] + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_multiple_soil_layers_3D(self, expected_geometry_two_layers_3D_extruded: Tuple[Geometry, Geometry], + create_default_3d_soil_material: SoilMaterial): + """ + Test if multiple soil layers are added correctly to the model in a 3D space. Multiple soil layers are generated + and multiple soil materials are created and added to the model. + + Args: + - expected_geometry_two_layers_3D_extruded (Tuple[:class:`stem.geometry.Geometry`, \ + :class:`stem.geometry.Geometry`]): expected geometry of the model which is created by extruding \ + a 2D geometry + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): default soil material + + """ + + ndim = 3 + + layer1_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + layer2_coordinates = [(1, 1, 0), (0, 1, 0), (0, 2, 0), (1, 2, 0)] + + # define soil materials + soil_material1 = create_default_3d_soil_material + soil_material1.name = "soil1" + + soil_material2 = create_default_3d_soil_material + soil_material2.name = "soil2" + + # create model + model = Model(ndim) + model.extrusion_length = [0, 0, 1] + + # add soil layers + model.add_soil_layer_by_coordinates(layer1_coordinates, soil_material1, "layer1") + model.add_soil_layer_by_coordinates(layer2_coordinates, soil_material2, "layer2") + + model.synchronise_geometry() + + # check if layers are added correctly + assert len(model.body_model_parts) == 2 + assert model.body_model_parts[0].name == "layer1" + assert model.body_model_parts[0].material == soil_material1 + assert model.body_model_parts[1].name == "layer2" + assert model.body_model_parts[1].material == soil_material2 + + # check if geometry is added correctly for each layer + for i in range(len(model.body_model_parts)): + generated_geometry = model.body_model_parts[i].geometry + expected_geometry = expected_geometry_two_layers_3D_extruded[i] + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_all_layers_from_geo_file_2D(self, expected_geometry_two_layers_2D: Tuple[Geometry, Geometry]): + """ + Tests if all layers are added correctly to the model in a 2D space. A geo file is read and all layers are + added to the model. + + Args: + - expected_geometry_two_layers_2D (Tuple[:class:`stem.geometry.Geometry`, :class:`stem.geometry.Geometry`]): \ + expected geometry of the model + + """ + + geo_file_name = "tests/test_data/gmsh_utils_two_blocks_2D.geo" + + # create model + model = Model(ndim=2) + model.add_all_layers_from_geo_file(geo_file_name, ["group_1"]) + + # check if body model parts are added correctly + assert len(model.body_model_parts) == 1 + assert model.body_model_parts[0].name == "group_1" + + # check if process model part is added correctly + assert len(model.process_model_parts) == 1 + assert model.process_model_parts[0].name == "group_2" + + # check if geometry is added correctly for each layer + for i in range(len(model.body_model_parts)): + generated_geometry = model.body_model_parts[i].geometry + expected_geometry = expected_geometry_two_layers_2D[i] + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_all_layers_from_geo_file_3D(self, expected_geometry_two_layers_3D_geo_file: Tuple[Geometry, Geometry]): + """ + Tests if all layers are added correctly to the model in a 3D space. A geo file is read and all layers are + added to the model. + + Args: + - expected_geometry_two_layers_3D_geo_file (Tuple[:class:`stem.geometry.Geometry`, \ + :class:`stem.geometry.Geometry`]): expected geometry of the model + + """ + + geo_file_name = "tests/test_data/gmsh_utils_column_3D_tetra4.geo" + + # create model + model = Model(ndim=3) + model.add_all_layers_from_geo_file(geo_file_name, ["group_1"]) + + # check if body model parts are added correctly + assert len(model.body_model_parts) == 1 + assert model.body_model_parts[0].name == "group_1" + + # check if process model part is added correctly + assert len(model.process_model_parts) == 1 + assert model.process_model_parts[0].name == "group_2" + + # check if geometry is added correctly + all_model_parts = [] + all_model_parts.extend(model.body_model_parts) + all_model_parts.extend(model.process_model_parts) + + # check if geometry is added correctly for each layer + for i in range(len(all_model_parts)): + generated_geometry = all_model_parts[i].geometry + expected_geometry = expected_geometry_two_layers_3D_geo_file[i] + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_synchronise_geometry_2D(self, expected_geometry_two_layers_2D_after_sync: Tuple[Geometry, Geometry], + create_default_2d_soil_material: SoilMaterial): + """ + Test if the geometry is synchronised correctly in 2D after adding a new layer to the model. Where the new layer + overlaps with the existing layer, the existing layer is cut and the overlapping part is removed. + + Args: + - expected_geometry_two_layers_2D_after_sync (Tuple[:class:`stem.geometry.Geometry`, \ + :class:`stem.geometry.Geometry`, :class:`stem.geometry.Geometry`]): The expected geometry after \ + synchronising the geometry. + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # define layer coordinates + ndim = 2 + layer1_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + layer2_coordinates = [(1, 1, 0), (0.5, 1, 0), (0.5, 2, 0), (1, 2, 0)] + + # define soil materials + soil_material1 = create_default_2d_soil_material + soil_material1.name = "soil1" + + soil_material2 = create_default_2d_soil_material + soil_material2.name = "soil2" + + # create model + model = Model(ndim) + + # add soil layers + model.add_soil_layer_by_coordinates(layer1_coordinates, soil_material1, "layer1") + model.add_soil_layer_by_coordinates(layer2_coordinates, soil_material2, "layer2") + + # synchronise geometry and recalculates the ids + model.synchronise_geometry() + + # collect all generated geometries + generated_geometries = [model.body_model_parts[0].geometry, model.body_model_parts[1].geometry, model.geometry] + + # check if geometry is added correctly for each layer + for generated_geometry, expected_geometry in zip(generated_geometries, + expected_geometry_two_layers_2D_after_sync): + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_synchronise_geometry_3D(self, create_default_3d_soil_material: SoilMaterial): + """ + Test if the geometry is synchronised correctly in 3D after adding a new layer to the model. Where the new layer + overlaps with the existing layer, the existing layer is cut and the overlapping part is removed. + + Args: + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # define layer coordinates + ndim = 3 + layer1_coordinates = [(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)] + layer2_coordinates = [(1, 1, 0), (0.5, 1, 0), (0.5, 2, 0), (1, 2, 0)] + + # define soil materials + soil_material1 = create_default_3d_soil_material + soil_material1.name = "soil1" + + soil_material2 = create_default_3d_soil_material + soil_material2.name = "soil2" + + # create model + model = Model(ndim) + model.extrusion_length = [0, 0, 1] + + # add soil layers + model.add_soil_layer_by_coordinates(layer1_coordinates, soil_material1, "layer1") + model.add_soil_layer_by_coordinates(layer2_coordinates, soil_material2, "layer2") + + # synchronise geometry and recalculates the ids + model.synchronise_geometry() + + with open("tests/test_data/expected_geometry_after_sync_3D.pickle", "rb") as f: + expected_geometry_two_layers_3D_after_sync = pickle.load(f) + + # collect all generated geometries + generated_geometries = [model.body_model_parts[0].geometry, model.body_model_parts[1].geometry, model.geometry] + + # check if geometry is added correctly for each layer + for generated_geometry, expected_geometry in zip(generated_geometries, + expected_geometry_two_layers_3D_after_sync): + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_point_loads_to_2_points(self, create_default_point_load_parameters: PointLoad): + """ + Test if a single soil point load is added correctly to the model. Two points are generated + and a single load is created and added to the model. + + Args: + - create_default_point_load_properties (:class:`stem.load.PointLoad`): default point load parameters + + """ + + ndim = 3 + + point_coordinates = [(-0.5, 0, 0), (0.5, 0, 0)] + + # define soil material + load_parameters = create_default_point_load_parameters + + # create model + model = Model(ndim) + # add point load + model.add_load_by_coordinates(point_coordinates, load_parameters, "point_load_1") + + # check if layer is added correctly + assert len(model.process_model_parts) == 1 + assert model.process_model_parts[0].name == "point_load_1" + TestUtils.assert_dictionary_almost_equal( + model.process_model_parts[0].parameters.__dict__, + load_parameters.__dict__ + ) + + # check if geometry is added correctly + generated_geometry = model.process_model_parts[0].geometry + expected_geometry = Geometry( + points={1:Point.create([-0.5, 0, 0], 1), 2: Point.create([0.5, 0, 0], 2)} + ) + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_line_load_to_3_edges(self, expected_geometry_line_load: Geometry, + create_default_line_load_parameters: PointLoad): + """ + Test if a line load is added correctly to the model when applied on 3 edges. 4 points are generated + and a single soil material is created and added to the model. + + Args: + - expected_geometry_line_load (:class:`stem.geometry.Geometry`): expected geometry of the model + - create_default_line_load_parameters (:class:`stem.load.LineLoad`): default line load parameters + + """ + + ndim = 3 + + point_coordinates = [(0, 0, 0), (3, 0, 0), (4, -1, 0), (10, -1, 0)] + + # define soil material + load_parameters = create_default_line_load_parameters + + # create model + model = Model(ndim) + # add line load + model.add_load_by_coordinates(point_coordinates, load_parameters, "line_load_1") + + # check if layer is added correctly + assert len(model.process_model_parts) == 1 + assert model.process_model_parts[0].name == "line_load_1" + TestUtils.assert_dictionary_almost_equal( + model.process_model_parts[0].parameters.__dict__, + load_parameters.__dict__ + ) + # check if geometry is added correctly + generated_geometry = model.process_model_parts[0].geometry + expected_geometry = expected_geometry_line_load + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_add_moving_point_load(self, expected_geometry_line_load: Geometry, + create_default_moving_load_parameters: MovingLoad): + """ + Test if a single soil point load is added correctly to the model. Two points are generated + and a single load is created and added to the model. + + Args: + - expected_geometry_line_load (:class:`stem.geometry.Geometry`): expected geometry of the model + - create_default_moving_load_parameters (:class:`stem.load.MovingLoad`): default moving load parameters + + """ + + ndim = 3 + + point_coordinates = [(0, 0, 0), (3, 0, 0), (4, -1, 0), (10, -1, 0)] + # origin is in (3.5, -0.5, 0) thus in the trajectory + + # define soil material + load_parameters = create_default_moving_load_parameters + + # create model + model = Model(ndim) + # add moving load + model.add_load_by_coordinates(point_coordinates, load_parameters, "moving_load_1") + + # check if layer is added correctly + assert len(model.process_model_parts) == 1 + assert model.process_model_parts[0].name == "moving_load_1" + TestUtils.assert_dictionary_almost_equal( + model.process_model_parts[0].parameters.__dict__, + load_parameters.__dict__ + ) + + # check if geometry is added correctly + generated_geometry = model.process_model_parts[0].geometry + expected_geometry = expected_geometry_line_load + + TestUtils.assert_almost_equal_geometries(expected_geometry, generated_geometry) + + def test_validation_coordinates(self): + """ + Test that validation raises and error if the points are not correctly specified. + """ + + ndim = 3 + model = Model(ndim=ndim) + + # test inputs for numpy arrays: + # test for 2D-array, correct number of coordinates (shape 3,2) + model.validate_coordinates(np.zeros((2,3))) + + # test for incorrect number of coordinates in array (shape 3,2) + with pytest.raises(ValueError, match=f"Coordinates should be 3D but 2 coordinates were given."): + model.validate_coordinates(np.zeros((3,2))) + + # test for incorrect number of dimension in array (1-D array) + with pytest.raises(ValueError, match=f"Coordinates are not a sequence of a sequence or a 2D array."): + model.validate_coordinates(np.arange(3)) + + # test inputs for sequence of floats: + # test for incorrect number of coordinates + with pytest.raises(ValueError, match=f"Coordinates should be 3D but 4 coordinates were given."): + model.validate_coordinates([(0.0, 0.0, 0.0, 4.0)]) + + # test for incorrect type (Sequence of float instead of Sequence[Sequence[float]]) + with pytest.raises(ValueError, match="Coordinates are not a sequence of a sequence or a 2D array."): + model.validate_coordinates([0.0, 0.0, 0.0]) + + def test_validation_moving_load(self, create_default_moving_load_parameters:MovingLoad): + """ + Test validation of moving load when points is not collinear to the trajectory. + + Args: + - create_default_moving_load_parameters (:class:`stem.load.MovingLoad`): default moving load parameters + + """ + + ndim = 3 + + point_coordinates = [(0.0, 0, 0), (1, 0, 0), (2, 0, 0), (4, 0, 0)] + # origin is in (1.5, 0.5, 0) thus not in the trajectory + + # define soil material + load_parameters = create_default_moving_load_parameters + # create model + model = Model(ndim) + + with pytest.raises(ValueError, match="Origin is not in the trajectory of the moving load."): + model.add_load_by_coordinates( + point_coordinates, load_parameters, "moving_load_1" + ) + + def test_generate_mesh_with_only_a_body_model_part_2d(self, create_default_2d_soil_material: SoilMaterial): + """ + Test if the mesh is generated correctly in 2D if there is only one body model part. + + Args: + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + model = Model(2) + + # add soil material + soil_material = create_default_2d_soil_material + + # add soil layers + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)], soil_material, "layer1") + model.synchronise_geometry() + + # generate mesh + model.generate_mesh() + + mesh = model.body_model_parts[0].mesh + + assert mesh.ndim == 2 + + unique_element_ids = [] + # check if mesh is generated correctly, i.e. if the number of elements is correct and if the element type is + # correct and if the element ids are unique and if the number of nodes per element is correct + assert len(mesh.elements) == 162 + for element in mesh.elements: + assert element.element_type == "TRIANGLE_3N" + assert element.id not in unique_element_ids + assert len(element.node_ids) == 3 + unique_element_ids.append(element.id) + + # check if nodes are generated correctly, i.e. if there are nodes in the mesh and if the node ids are unique + # and if the number of coordinates per node is correct + unique_node_ids = [] + assert len(mesh.nodes) == 98 + for node in mesh.nodes: + assert node.id not in unique_node_ids + assert len(node.coordinates) == 3 + unique_node_ids.append(node.id) + + def test_generate_mesh_with_only_a_body_model_part_3d(self, create_default_3d_soil_material: SoilMaterial): + """ + Test if the mesh is generated correctly in 3D if there is only one body model part. + + Args: + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + model = Model(3) + model.extrusion_length = [0, 0, 1] + + # add soil material + soil_material = create_default_3d_soil_material + + # add soil layers + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)], soil_material, "layer1") + model.synchronise_geometry() + + # generate mesh + model.generate_mesh() + + mesh = model.body_model_parts[0].mesh + + assert mesh.ndim == 3 + + unique_element_ids = [] + # check if mesh is generated correctly, i.e. if the number of elements is correct and if the element type is + # correct and if the element ids are unique and if the number of nodes per element is correct + assert len(mesh.elements) == 1120 + for element in mesh.elements: + assert element.element_type == "TETRAHEDRON_4N" + assert element.id not in unique_element_ids + assert len(element.node_ids) == 4 + unique_element_ids.append(element.id) + + # check if nodes are generated correctly, i.e. if there are nodes in the mesh and if the node ids are unique + # and if the number of coordinates per node is correct + unique_node_ids = [] + assert len(mesh.nodes) == 340 + for node in mesh.nodes: + assert node.id not in unique_node_ids + assert len(node.coordinates) == 3 + unique_node_ids.append(node.id) + + def test_generate_mesh_with_body_and_process_model_part(self, create_default_2d_soil_material: SoilMaterial): + """ + Test if the mesh is generated correctly in the body model part and a process model part. + + Args: + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + """ + model = Model(2) + + # add soil material + soil_material = create_default_2d_soil_material + + # add soil layers + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)], soil_material, "layer1") + + # add process geometry + gmsh_process_input = {"process_0d": {"coordinates": [[0, 0.5, 0]], "ndim": 0}} + model.gmsh_io.generate_geometry(gmsh_process_input, "") + + # create process model part + process_model_part = ModelPart("process_0d") + + # set the geometry of the process model part + process_model_part.get_geometry_from_geo_data(model.gmsh_io.geo_data, "process_0d") + + # add process model part + model.process_model_parts.append(process_model_part) + + # synchronise geometry and generate mesh + model.synchronise_geometry() + model.generate_mesh() + + # check mesh of body model part + mesh_body = model.body_model_parts[0].mesh + + assert mesh_body.ndim == 2 + + unique_element_ids = [] + # check if mesh is generated correctly, i.e. if the number of elements is correct and if the element type is + # correct and if the element ids are unique and if the number of nodes per element is correct + assert len(mesh_body.elements) == 162 + for element in mesh_body.elements: + assert element.element_type == "TRIANGLE_3N" + assert element.id not in unique_element_ids + assert len(element.node_ids) == 3 + unique_element_ids.append(element.id) + + # check if nodes are generated correctly, i.e. if there are nodes in the mesh and if the node ids are unique + # and if the number of coordinates per node is correct + unique_body_node_ids = [] + assert len(mesh_body.nodes) == 98 + for node in mesh_body.nodes: + assert node.id not in unique_body_node_ids + assert len(node.coordinates) == 3 + unique_body_node_ids.append(node.id) + + # check process model part + mesh_process = model.process_model_parts[0].mesh + + assert mesh_process.ndim == 0 + + # check elements of process model part, i.e. if the number of elements is correct and if the element type is + # correct and if the element ids are unique and if the number of nodes per element is correct + assert len(mesh_process.elements) == 1 + for element in mesh_process.elements: + assert element.element_type == "POINT_1N" + assert element.id == 1 + assert element.id not in unique_element_ids + assert len(element.node_ids) == 1 + unique_element_ids.append(element.id) + + # check nodes of process model part, i.e. if there is 1 node in the mesh and if the node ids are present in the + # body mesh and if the number of coordinates per node is correct + assert len(mesh_process.nodes) == 1 + for node in mesh_process.nodes: + + # check if node is also available in the body mesh + assert node.id in unique_body_node_ids + assert len(node.coordinates) == 3 + + def test_validate_expected_success(self): + """ + Test if the model is validated correctly. A model is created with two process model parts which both have + a unique name. + + """ + + model = Model(2) + + model_part1 = ModelPart("test1") + model_part2 = ModelPart("test2") + + model.process_model_parts = [model_part1, model_part2] + + model.validate() + + def test_validate_expected_fail_non_unique_names(self): + """ + Test if the model is validated correctly. A model is created with two process model parts which both have + the same name. This should raise a ValueError. + + """ + + model = Model(2) + + model_part1 = ModelPart("test") + model_part2 = ModelPart("test") + + model.process_model_parts = [model_part1, model_part2] + + pytest.raises(ValueError, model.validate) + + def test_validate_expected_fail_no_name(self): + """ + Test if the model is validated correctly. A model is created with a process model part which does not contain + a name. This should raise a ValueError. + + """ + + model = Model(2) + + model_part1 = ModelPart(None) + model.process_model_parts = [model_part1] + + pytest.raises(ValueError, model.validate) + + def test_add_boundary_condition_by_geometry_ids(self,create_default_3d_soil_material: SoilMaterial): + """ + Test if a boundary condition is added correctly to the model. A boundary condition is added to the model by + specifying the geometry ids to which the boundary condition should be applied. + + Args: + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + + # create a 3D model + model = Model(3) + model.extrusion_length = [0, 0, 1] + + # create multiple boundary condition parameters + no_rotation_parameters = RotationConstraint(active=[True, True, True], is_fixed=[True, True, True], + value=[0, 0, 0]) + + absorbing_parameters = AbsorbingBoundary(absorbing_factors=[1,1], virtual_thickness=0) + + no_displacement_parameters = DisplacementConstraint(active=[True, True, True], is_fixed=[True, True, True], + value=[0, 0, 0]) + + # add body model part + soil_material = create_default_3d_soil_material + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0)], soil_material, "test_soil") + + + # add boundary conditions in 0d, 1d and 2d + model.add_boundary_condition_by_geometry_ids(0, [1, 2], no_rotation_parameters, "no_rotation") + model.add_boundary_condition_by_geometry_ids(1, [8], absorbing_parameters, "absorbing") + model.add_boundary_condition_by_geometry_ids(2, [1, 2], no_displacement_parameters, "no_displacement") + + model.synchronise_geometry() + + # set expected parameters of the boundary conditions + expected_0d_model_part_parameters = RotationConstraint(active=[True, True, True], is_fixed=[True, True, True], + value=[0, 0, 0]) + + expected_1d_model_part_parameters = AbsorbingBoundary(absorbing_factors=[1, 1], virtual_thickness=0) + + expected_2d_model_part_parameters = DisplacementConstraint(active=[True, True, True], + is_fixed=[True, True, True], value=[0, 0, 0]) + + # set expected geometry 0d boundary condition + expected_boundary_points = {1: Point.create([0, 0, 0], 1), 2: Point.create([1, 0, 0], 2)} + expected_boundary_lines = {1: Line.create([1, 2], 1)} + expected_boundary_surfaces = {} + expected_boundary_volumes = {} + + expected_boundary_geometry_0d = Geometry(expected_boundary_points, expected_boundary_lines, + expected_boundary_surfaces, expected_boundary_volumes) + + # set expected geometry 1d boundary condition + expected_boundary_points = {3: Point.create([1, 1, 0], 3), 7: Point.create([1, 1, 1], 7)} + expected_boundary_lines = {8: Line.create([3, 7], 8)} + expected_boundary_surfaces = {} + expected_boundary_volumes = {} + + expected_boundary_geometry_1d = Geometry(expected_boundary_points, expected_boundary_lines, + expected_boundary_surfaces, expected_boundary_volumes) + + # set expected geometry 2d boundary condition + + expected_boundary_geometry_2d = Geometry() + expected_boundary_geometry_2d.points = { + 1: Point.create([0, 0, 0], 1), + 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3), + 4: Point.create([0, 1, 0], 4), + 5: Point.create([0, 0, 1], 5), + 6: Point.create([1, 0, 1], 6) + } + + expected_boundary_geometry_2d.lines = { + 1: Line.create([1, 2], 1), + 2: Line.create([2, 3], 2), + 3: Line.create([3, 4], 3), + 4: Line.create([4, 1], 4), + 5: Line.create([1, 5], 5), + 7: Line.create([5, 6], 7), + 6: Line.create([2, 6], 6) + } + + expected_boundary_geometry_2d.surfaces = { + 1: Surface.create([1, 2, 3, 4], 1), + 2: Surface.create([5, 7, -6, -1], 2) + } + + expected_boundary_geometry_2d.volumes = {} + + # collect all expected geometries + all_expected_geometries = [expected_boundary_geometry_0d, expected_boundary_geometry_1d, + expected_boundary_geometry_2d] + + # check 0d parameters + npt.assert_allclose(model.process_model_parts[0].parameters.active, expected_0d_model_part_parameters.active) + npt.assert_allclose(model.process_model_parts[0].parameters.is_fixed, expected_0d_model_part_parameters.is_fixed) + npt.assert_allclose(model.process_model_parts[0].parameters.value, expected_0d_model_part_parameters.value) + + # check 1d parameters + npt.assert_allclose(model.process_model_parts[1].parameters.absorbing_factors, + expected_1d_model_part_parameters.absorbing_factors) + npt.assert_allclose(model.process_model_parts[1].parameters.virtual_thickness, + expected_1d_model_part_parameters.virtual_thickness) + + # check 2d parameters + npt.assert_allclose(model.process_model_parts[2].parameters.active, expected_2d_model_part_parameters.active) + npt.assert_allclose(model.process_model_parts[2].parameters.is_fixed, expected_2d_model_part_parameters.is_fixed) + npt.assert_allclose(model.process_model_parts[2].parameters.value, expected_2d_model_part_parameters.value) + + for expected_geometry, model_part in zip(all_expected_geometries, model.process_model_parts): + + TestUtils.assert_almost_equal_geometries(expected_geometry, model_part.geometry) + + + def test_add_gravity_load_1d_and_2d(self, create_default_2d_soil_material: SoilMaterial): + """ + Test if a gravity load is added correctly to the model in a 2d space containing 1d and 2d elements. A gravity + load is generated and added to the model. + + Args: + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # create model + model = Model(2) + + # add a 2d layer + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], create_default_2d_soil_material, "soil1") + + # add a 1d layer + layer_settings = {"beam": {"ndim": 1, + "element_size": -1, + "coordinates": [[0, 0, 0], [1, 0, 0]]}} + + model.gmsh_io.generate_geometry(layer_settings, "") + model.synchronise_geometry() + + # add 1d model part to model + body_model_part = BodyModelPart("beam") + body_model_part.material = EulerBeam(ndim=2, YOUNG_MODULUS=1e6, POISSON_RATIO=0.3, DENSITY=1, CROSS_AREA=1, + I33=1) + body_model_part.get_geometry_from_geo_data(model.gmsh_io.geo_data, "beam") + + model.body_model_parts.append(body_model_part) + + # add gravity load + model._Model__add_gravity_load() + + assert len(model.process_model_parts) == 2 + assert model.process_model_parts[0].name == "gravity_load_1d" + assert model.process_model_parts[1].name == "gravity_load_2d" + + # setup expected geometries for 1d and 2d + expected_geometry_points_1d = {1: Point.create([0, 0, 0],1), 2: Point.create([1, 0, 0], 2)} + expected_geometry_lines_1d = {1: Line.create([1, 2], 1)} + expected_geometry_gravity_1d = Geometry(expected_geometry_points_1d, expected_geometry_lines_1d, {}, {}) + + expected_geometry_points_2d = {1: Point.create([0, 0, 0], 1), 2: Point.create([1, 0, 0], 2), + 3: Point.create([1, 1, 0], 3)} + expected_geometry_lines_2d = {1: Line.create([1, 2], 1), 2: Line.create([2, 3], 2), 3: Line.create([3, 1], 3)} + expected_geometry_surfaces_2d = {1: Surface.create([1, 2, 3], 1)} + expected_geometry_gravity_2d = Geometry(expected_geometry_points_2d, expected_geometry_lines_2d, + expected_geometry_surfaces_2d, {}) + + expected_geometries = [expected_geometry_gravity_1d, expected_geometry_gravity_2d] + + # check if all process model parts are correct + for model_part in model.process_model_parts: + + # check if parameters are added correctly + npt.assert_allclose(model_part.parameters.value, [0, -9.81, 0]) + npt.assert_allclose(model_part.parameters.active, [True, True, True]) + + # check if geometry is added correctly + generated_model_part = model_part.geometry + + TestUtils.assert_almost_equal_geometries(expected_geometries[0], generated_model_part) + + def test_add_gravity_load_two_layers_same_dimension(self, create_default_2d_soil_material: SoilMaterial): + """ + Test if a gravity load is added correctly to the model in a 2d space containing 2 layers. A gravity load is + generated and added to the model. + + Args: + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # create model + model = Model(2) + + # add a 2d layer + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], create_default_2d_soil_material, "soil1") + model.add_soil_layer_by_coordinates([(1, 0, 0), (0, 0, 0), (1, -1, 0)], create_default_2d_soil_material, "soil2") + + model.synchronise_geometry() + + # add gravity load + model._Model__add_gravity_load(-12,0) + + assert len(model.process_model_parts) == 1 + + generated_geometry = model.process_model_parts[0].geometry + + # check if number of points, lines, surfaces are correct, i.e. if the number of points, lines, surfaces are the + # same as the number of points, lines, surfaces of the model geometry + assert len(generated_geometry.points) == len(model.geometry.points) == 4 + assert len(generated_geometry.lines) == len(model.geometry.lines) == 5 + assert len(generated_geometry.surfaces) == len(model.geometry.surfaces) == 2 + + assert model.process_model_parts[0].name == "gravity_load_2d" + npt.assert_allclose(model.process_model_parts[0].parameters.value, [-12, 0, 0]) + npt.assert_allclose(model.process_model_parts[0].parameters.active, [True, True, True]) + + def test_add_gravity_load_3d(self, create_default_3d_soil_material): + """ + Test if a gravity load is added correctly to the model in a 3d space. A gravity load is generated and added to + the model. + + Args: + - create_default_3d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # create model + model = Model(3) + model.extrusion_length = [0, 0, 1] + + # add a 2d layer + model.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], create_default_3d_soil_material, "soil1") + + model.synchronise_geometry() + + # add gravity load + model._Model__add_gravity_load(vertical_axis=2, gravity_value=-10) + + assert len(model.process_model_parts) == 1 + + generated_geometry = model.process_model_parts[0].geometry + + # check if number of points, lines, surfaces are correct, i.e. if the number of points, lines, surfaces and + # volumes are the same as the number of points, lines, surfaces and volumes of the model geometry + assert len(generated_geometry.points) == len(model.geometry.points) == 6 + assert len(generated_geometry.lines) == len(model.geometry.lines) == 9 + assert len(generated_geometry.surfaces) == len(model.geometry.surfaces) == 5 + assert len(generated_geometry.volumes) == len(model.geometry.volumes) == 1 + + assert model.process_model_parts[0].name == "gravity_load_3d" + npt.assert_allclose(model.process_model_parts[0].parameters.value, [0, 0, -10]) + npt.assert_allclose(model.process_model_parts[0].parameters.active, [True, True, True]) + + def test_setup_stress_initialisation(self, create_default_2d_soil_material: SoilMaterial): + """ + Test if the stress initialisation is set up correctly. A model is created with a soil layer. It is checked if + gravity is added in case the K0 procedure or gravity loading is used. + + Args: + - create_default_2d_soil_material (:class:`stem.soil_material.SoilMaterial`): A default soil material. + + """ + + # set up solver settings + analysis_type = AnalysisType.MECHANICAL_GROUNDWATER_FLOW + + solution_type = SolutionType.QUASI_STATIC + + time_integration = TimeIntegration(start_time=0.0, end_time=1.0, delta_time=0.1, reduction_factor=0.5, + increase_factor=2.0, max_delta_time_factor=500) + + convergence_criterion = DisplacementConvergenceCriteria() + + stress_initialisation_type = StressInitialisationType.NONE + + solver_settings = SolverSettings(analysis_type=analysis_type, solution_type=solution_type, + stress_initialisation_type=stress_initialisation_type, + time_integration=time_integration, + is_stiffness_matrix_constant=True, are_mass_and_damping_constant=True, + convergence_criteria=convergence_criterion) + + # set up problem data + problem_data = Problem(problem_name="test", number_of_threads=2, settings=solver_settings) + + model_no_gravity = Model(2) + model_no_gravity.project_parameters = problem_data + + # set up soil material + soil_material = create_default_2d_soil_material + model_no_gravity.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], soil_material, "soil1") + model_no_gravity.synchronise_geometry() + + # setup_stress_initialisation + model_no_gravity._Model__setup_stress_initialisation() + + model_k0 = Model(2) + model_k0.project_parameters = problem_data + + model_k0.project_parameters.settings.stress_initialisation_type = StressInitialisationType.K0_PROCEDURE + model_k0.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], soil_material, "soil1") + model_k0.synchronise_geometry() + + # setup_stress_initialisation + model_k0._Model__setup_stress_initialisation() + + model_gravity_loading = Model(2) + model_gravity_loading.project_parameters = problem_data + + model_gravity_loading.project_parameters.settings.stress_initialisation_type = \ + StressInitialisationType.GRAVITY_LOADING + model_gravity_loading.add_soil_layer_by_coordinates([(0, 0, 0), (1, 0, 0), (1, 1, 0)], soil_material, "soil1") + model_gravity_loading.synchronise_geometry() + + # setup_stress_initialisation + model_gravity_loading._Model__setup_stress_initialisation() + + assert len(model_no_gravity.process_model_parts) == 0 + assert len(model_k0.process_model_parts) == 1 + assert len(model_gravity_loading.process_model_parts) == 1 + + assert model_k0.process_model_parts[0].name == "gravity_load_2d" + assert model_gravity_loading.process_model_parts[0].name == "gravity_load_2d" + + def test_setup_stress_initialisation_without_project_parameters(self): + """ + A model is created without project parameters. It is + checked if a ValueError is raised while setting up the stress initialisation. + + """ + # create model + model = Model(2) + + # test if value error is raised + with pytest.raises(ValueError, + match=r"Project parameters must be set before setting up the stress initialisation"): + model._Model__setup_stress_initialisation() + + @pytest.mark.skip("Not implemented yet") + def test_post_setup(self): + pass diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 000000000..2f48dfd5b --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,83 @@ + +import numpy as np + +from stem.utils import Utils + + +class TestUtilsStem: + + def test_is_clockwise(self): + """ + Test the check which checks if coordinates are given in clockwise order + """ + + coordinates = [[0, 0], [2, 0], [2, 2], [0, 2]] + + assert not Utils.are_2d_coordinates_clockwise(coordinates=coordinates) + assert Utils.are_2d_coordinates_clockwise(coordinates=coordinates[::-1]) + + def test_is_clockwise_non_convex(self): + """ + Test the check which checks if coordinates are given in clockwise order for a non-convex polygon + + """ + + coordinates = [[0, 0], [2, 0], [2, 2], [1, -1], [0, 2]] + + assert not Utils.are_2d_coordinates_clockwise(coordinates=coordinates) + assert Utils.are_2d_coordinates_clockwise(coordinates=coordinates[::-1]) + + def test_collinearity_2d(self): + """ + Check collinearity between 3 points in 2D + """ + p1 = np.array([0, 0]) + p2 = np.array([-2, -1]) + + p_test_1 = np.array([2, 1]) + p_test_2 = np.array([-5, 1]) + + assert Utils.is_collinear(point=p_test_1, start_point=p1, end_point=p2) + assert not Utils.is_collinear(point=p_test_2, start_point=p1, end_point=p2) + + def test_collinearity_3d(self): + """ + Check collinearity between 3 points in 3D + """ + + p1 = np.array([0, 0, 0]) + p2 = np.array([-2, -2, 2]) + + p_test_1 = np.array([2, 2, -2]) + p_test_2 = np.array([2, -2, 2]) + + assert Utils.is_collinear(point=p_test_1, start_point=p1, end_point=p2) + assert not Utils.is_collinear(point=p_test_2, start_point=p1, end_point=p2) + + def test_is_in_between_2d(self): + """ + Check if point is in between other 2 points in 2D + """ + + p1 = np.array([0, 0]) + p2 = np.array([-2, -2]) + + p_test_1 = np.array([2, 2]) + p_test_2 = np.array([-1, -1]) + + assert not Utils.is_point_between_points(point=p_test_1, start_point=p1, end_point=p2) + assert Utils.is_point_between_points(point=p_test_2, start_point=p1, end_point=p2) + + def test_is_in_between_3d(self): + """ + Check if point is in between other 2 points in 3D + """ + + p1 = np.array([0, 0, 0]) + p2 = np.array([-2, -2, 2]) + + p_test_1 = np.array([2, 2, -2]) + p_test_2 = np.array([-1, -1, 1]) + + assert not Utils.is_point_between_points(point=p_test_1, start_point=p1, end_point=p2) + assert Utils.is_point_between_points(point=p_test_2, start_point=p1, end_point=p2) \ No newline at end of file diff --git a/tests/test_water_boundaries.py b/tests/test_water_boundaries.py new file mode 100644 index 000000000..eb38ed2a3 --- /dev/null +++ b/tests/test_water_boundaries.py @@ -0,0 +1,44 @@ +import pytest + +from stem.water_boundaries import * + + +class TestWaterBoundaries: + + def test_raise_errors_for_water_boundaries(self): + + pytest.raises(ValueError, + PhreaticMultiLineBoundary, + x_coordinates=[0, 1, 2], + y_coordinates=[0, 1, 2, 3], + surfaces_assigment=["surface_1", "surface_2", "surface_3"], + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + specific_weight=9.81, + water_pressure=1000) + + pytest.raises(ValueError, + PhreaticMultiLineBoundary, + x_coordinates=[0, 1, 2, 3, 4], + y_coordinates=[0, 1, 2, 3], + surfaces_assigment=["surface_1", "surface_2", "surface_3"], + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + specific_weight=9.81, + water_pressure=1000 + ) + + pytest.raises(ValueError, + PhreaticMultiLineBoundary, + x_coordinates=[0, 1, 2, 3], + y_coordinates=[0, 1, 2, 3], + z_coordinates=[0, 1, 2, 3, 4], + surfaces_assigment=["surface_1", "surface_2", "surface_3"], + is_fixed=True, + gravity_direction=1, + out_of_plane_direction=2, + specific_weight=9.81, + water_pressure=1000 + ) \ No newline at end of file diff --git a/tests/utils.py b/tests/utils.py index 99044283f..acfd4cb4a 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -2,6 +2,9 @@ from typing import Dict, Any import numpy.testing as npt +import pytest + +from stem.geometry import Geometry class TestUtils: @@ -12,8 +15,8 @@ def assert_dictionary_almost_equal(expected: Dict[Any, Any], actual: Dict[Any, A Checks whether two dictionaries are equal. Args: - expected: Expected dictionary. - actual: Actual dictionary. + - expected: Expected dictionary. + - actual: Actual dictionary. """ @@ -36,4 +39,44 @@ def assert_dictionary_almost_equal(expected: Dict[Any, Any], actual: Dict[Any, A npt.assert_allclose(v_i, actual_i) else: - npt.assert_allclose(v, actual[k]) \ No newline at end of file + npt.assert_allclose(v, actual[k]) + + @staticmethod + def assert_almost_equal_geometries(expected_geometry: Geometry, actual_geometry:Geometry): + """ + Checks whether two Geometries are (almost) equal. + + Args: + - expected_geometry (:class:`stem.geometry.Geometry`): expected geometry of the model + - actual_geometry (:class:`stem.geometry.Geometry`): actual geometry of the model + + Returns: + + """ + # check if points are added correctly + for (generated_point_id, generated_point), (expected_point_id, expected_point) in \ + zip(actual_geometry.points.items(), expected_geometry.points.items()): + assert generated_point_id == expected_point_id + assert generated_point.id == expected_point.id + npt.assert_allclose(generated_point.coordinates, expected_point.coordinates) + + # check if lines are added correctly + for (generated_lines_id, generated_line), (expected_line_id, expected_line) in \ + zip(actual_geometry.lines.items(), expected_geometry.lines.items()): + assert generated_lines_id == expected_line_id + assert generated_line.id == expected_line.id + npt.assert_equal(generated_line.point_ids, expected_line.point_ids) + + # check if surfaces are added correctly + for (generated_surface_id, generated_surface), (expected_surface_id, expected_surface) in \ + zip(actual_geometry.surfaces.items(), expected_geometry.surfaces.items()): + assert generated_surface_id == expected_surface_id + assert generated_surface.id == expected_surface.id + npt.assert_equal(generated_surface.line_ids, expected_surface.line_ids) + + # check if volumes are added correctly + for (generated_volume_id, generated_volume), (expected_volume_id, expected_volume) in \ + zip(actual_geometry.volumes.items(), expected_geometry.volumes.items()): + assert generated_volume_id == expected_volume_id + assert generated_volume.id == expected_volume.id + npt.assert_equal(generated_volume.surface_ids, expected_volume.surface_ids)