Allow importing meshes from file

conmech/mesh/mesh.py

@@ -115,6 +115,8 @@ def reinitialize_data(
             mesh_prop=mesh_prop,
             create_in_subprocess=create_in_subprocess,
         )
+        if mesh_prop.mesh_type == "msh_file":
+            self.fill_mesh_prop_data(input_nodes)
         unordered_nodes, unordered_elements = remove_unconnected_nodes_numba(
             input_nodes, input_elements
         )
@@ -131,6 +133,10 @@ def reinitialize_data(
         edges_matrix = get_edges_matrix(nodes_count=len(self.initial_nodes), elements=self.elements)
         self.edges = get_edges_list_numba(edges_matrix)
 
+    def fill_mesh_prop_data(self, nodes: np.ndarray):
+        self.mesh_prop.dimension = nodes.shape[1]
+        self.mesh_prop.scale = [max(nodes[:, j]) - min(nodes[:, j]) for j in range(nodes.shape[1])]
+
     def normalize_shift(self, vectors):
         _ = self
         return vectors - np.mean(vectors, axis=0)

conmech/mesh/mesh_builders.py

@@ -1,6 +1,7 @@
 from typing import Tuple
 
 import numpy as np
+import meshio
 
 from conmech.helpers import mph, nph
 from conmech.mesh.zoo import MeshZOO
@@ -40,6 +41,10 @@ def build_initial_mesh(
     mesh_prop: MeshProperties,
     create_in_subprocess=False,
 ) -> RawMesh:
+    if mesh_prop.mesh_type == "msh_file":
+        mesh = meshio.read(mesh_prop.path)
+        return RawMesh(nodes=mesh.points, elements=mesh.cells_dict["triangle"])
+
     if not create_in_subprocess:
         return MeshZOO.get_by_name(mesh_prop.mesh_type)(mesh_prop)
     return mph.run_process(lambda: MeshZOO.get_by_name(mesh_prop.mesh_type)(mesh_prop))

conmech/properties/mesh_properties.py

@@ -7,15 +7,37 @@
 @dataclass
 class MeshProperties:
     mesh_type: str
-    mesh_density: List[float]
-    scale: List[float]
-    dimension: int
+
+    mesh_density: Optional[List[float]] = None
+    grid_height: Optional[float] = None
+    scale: Optional[List[float]] = None
+    dimension: Optional[int] = None
+
+    path: Optional[str] = None
+
     initial_base: Optional[np.ndarray] = None
     initial_position: Optional[np.ndarray] = None
     mean_at_origin: bool = False
     corners_vector: Optional[np.ndarray] = None
     corner_mesh_data: Optional[np.ndarray] = None
 
+    def __post_init__(self):
+        self.do_sanity_checks()
+        if self.mesh_type != "msh_file":
+            self.scale = self.scale or [
+                (self.grid_height / self.mesh_density[1]) * elems_num
+                for elems_num in self.mesh_density
+            ]
+            self.dimension = self.dimension or 2
+
+    def do_sanity_checks(self):
+        if not (
+            self.mesh_type == "msh_file" and self.path is not None and self.mesh_density is None
+        ) and not (
+            self.mesh_type != "msh_file" and self.path is None and self.mesh_density is not None
+        ):
+            raise ValueError("Improper combination of mesh parameters")
+
     @staticmethod
     def _get_modulo(array, index):
         return array[index % len(array)]

conmech/scenarios/problems.py

@@ -3,10 +3,11 @@
 """
 from abc import ABC
 from dataclasses import dataclass
-from typing import Tuple, Union, Optional, Callable
+from typing import Optional, Callable
 
 import numpy as np
 
+from conmech.mesh.mesh import MeshProperties
 from conmech.mesh.boundaries_description import BoundariesDescription
 
 
@@ -35,13 +36,9 @@ def regularized_subderivative_tangential_direction(
 @dataclass
 class Problem(ABC):
     # pylint: disable=unused-argument
-    dimension = 2  # TODO #74 : Not used?
-    mesh_type: str
-    grid_height: float
+    mesh_prop: MeshProperties
     boundaries: BoundariesDescription
 
-    elements_number: Union[Tuple[int, int], Tuple[int, int, int]]  # number of triangles per aside
-
     @staticmethod
     def inner_forces(x: np.ndarray, t: Optional[float] = None) -> np.ndarray:
         return np.zeros_like(x)
@@ -85,8 +82,6 @@ def external_temperature(x: np.ndarray, t: Optional[float] = None) -> np.ndarray
 
 @dataclass
 class DisplacementProblem(Problem, ABC):
-    elements_number: Union[Tuple[int, int], Tuple[int, int, int]]  # number of triangles per aside
-
     mu_coef: float
     la_coef: float
     contact_law: ContactLaw

conmech/simulations/problem_solver.py

@@ -26,7 +26,6 @@
     ViscoelasticPiezoelectricProperties,
     ViscoelasticTemperatureProperties,
 )
-from conmech.properties.mesh_properties import MeshProperties
 from conmech.properties.schedule import Schedule
 from conmech.scenarios.problems import (
     TimeDependentProblem,
@@ -65,17 +64,8 @@ def __init__(self, problem: Problem, body_properties: BodyProperties):
         else:
             self.time_step = 0
 
-        grid_width: float = (
-            problem.grid_height / problem.elements_number[0]
-        ) * problem.elements_number[1]
-
         self.body: BodyForces = BodyForces(
-            mesh_prop=MeshProperties(
-                dimension=2,
-                mesh_type=problem.mesh_type,
-                mesh_density=[problem.elements_number[1], problem.elements_number[0]],
-                scale=[float(grid_width), float(problem.grid_height)],
-            ),
+            mesh_prop=problem.mesh_prop,
             body_prop=body_properties,
             schedule=Schedule(time_step=self.time_step, final_time=0.0),
             boundaries_description=problem.boundaries,

examples/Jureczka_Ochal_Bartman_2023.py

@@ -13,6 +13,7 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import ContactLaw, QuasistaticDisplacementProblem
 from conmech.simulations.problem_solver import TimeDependentSolver as QuasistaticProblemSolver
+from conmech.mesh.mesh import MeshProperties
 from examples.utils import viscoelastic_constitutive_law
 
 
@@ -56,8 +57,6 @@ def regularized_subderivative_tangential_direction(
 def make_setup(mesh_type_, boundaries_, contact_law_, elements_number_, friction_bound_):
     @dataclass()
     class QuasistaticSetup(QuasistaticDisplacementProblem):
-        grid_height: ... = 1.0
-        elements_number: ... = elements_number_
         mu_coef: ... = 40
         la_coef: ... = 100
         th_coef: ... = 40
@@ -89,7 +88,8 @@ def friction_bound(u_nu: float) -> float:
 
         boundaries: ... = boundaries_
 
-    return QuasistaticSetup(mesh_type=mesh_type_)
+    mesh_prop = MeshProperties(mesh_type=mesh_type_, mesh_density=elements_number_, grid_height=1.0)
+    return QuasistaticSetup(mesh_prop)
 
 
 def main(config: Config):

examples/Jureczka_and_Ochal_2019.py

@@ -9,6 +9,7 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import ContactLaw, StaticDisplacementProblem
 from conmech.simulations.problem_solver import StaticSolver as StaticProblemSolver
+from conmech.mesh.mesh import MeshProperties
 
 
 class JureczkaOchal2019(ContactLaw):
@@ -43,8 +44,6 @@ def regularized_subderivative_tangential_direction(
 
 @dataclass
 class StaticSetup(StaticDisplacementProblem):
-    grid_height: ... = 1.0
-    elements_number: ... = (8, 16)
     mu_coef: ... = 4
     la_coef: ... = 4
     contact_law: ... = JureczkaOchal2019
@@ -76,9 +75,11 @@ def main(config: Config):
 
     To see result of simulation you need to call from python `main(Config().init())`.
     """
-    setup = StaticSetup(mesh_type="cross")
+    mesh_prop = MeshProperties(mesh_type="cross", mesh_density=[16, 8], grid_height=1)
     if config.test:
-        setup.elements_number = (2, 4)
+        mesh_prop.mesh_density = [4, 2]
+
+    setup = StaticSetup(mesh_prop)
     runner = StaticProblemSolver(setup, "schur")
 
     state = runner.solve(

examples/Sofonea_Ochal_Bartman_2023.py

@@ -29,6 +29,7 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import RelaxationQuasistaticProblem
 from conmech.simulations.problem_solver import QuasistaticRelaxation
+from conmech.mesh.mesh import MeshProperties
 
 from examples.p_slope_contact_law import make_const_contact_law
 from examples.utils import elastic_relaxation_constitutive_law
@@ -46,8 +47,6 @@
 
 @dataclass
 class QuasistaticSetup(RelaxationQuasistaticProblem):
-    grid_height: ... = 1.0
-    elements_number: ... = (20, 20)
     mu_coef: ... = E / (1 + kappa)
     la_coef: ... = E * kappa / ((1 + kappa) * (1 - 2 * kappa))
     time_step: ... = 1 / 128
@@ -86,11 +85,12 @@ def main(config: Config):
 
     To see result of simulation you need to call from python `main(Config().init())`.
     """
-    setup = QuasistaticSetup(mesh_type="tunnel")
+    mesh_prop = MeshProperties(mesh_type="tunnel", grid_height=1.0, mesh_density=[20, 20])
+    setup = QuasistaticSetup(mesh_prop)
     if config.test:
-        setup.elements_number = (8, 8)
+        mesh_prop.mesh_density = [8, 8]
         setup.time_step *= 8
-    h = setup.elements_number[0]
+    h = mesh_prop.mesh_density[1]
 
     def sin_outer_forces(x, t=None):
         if x[1] <= oy:

examples/example_dynamic.py

@@ -10,13 +10,12 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import DynamicDisplacementProblem
 from conmech.simulations.problem_solver import TimeDependentSolver
+from conmech.mesh.mesh import MeshProperties
 from examples.p_slope_contact_law import make_slope_contact_law
 
 
 @dataclass()
 class DynamicSetup(DynamicDisplacementProblem):
-    grid_height: ... = 1.0
-    elements_number: ... = (2, 5)
     boundaries: ... = BoundariesDescription(
         contact=lambda x: x[1] == 0, dirichlet=lambda x: x[0] == 0
     )
@@ -47,7 +46,8 @@ def main(config: Config):
     To see result of simulation you need to call from python `main(Config().init())`.
     """
 
-    setup = DynamicSetup(mesh_type="cross")
+    mesh_prop = MeshProperties(mesh_type="cross", mesh_density=[5, 2], grid_height=1.0)
+    setup = DynamicSetup(mesh_prop)
     runner = TimeDependentSolver(setup, solving_method="schur")
     n_steps = 32 if not config.test else 10
 

examples/example_imported_mesh.py

@@ -0,0 +1,60 @@
+"""
+Created at 05.09.2023
+"""
+from dataclasses import dataclass
+
+import numpy as np
+from conmech.helpers.config import Config
+from conmech.mesh.boundaries_description import BoundariesDescription
+from conmech.plotting.drawer import Drawer
+from conmech.scenarios.problems import StaticDisplacementProblem
+from conmech.simulations.problem_solver import StaticSolver
+from conmech.mesh.mesh import MeshProperties
+
+
+from examples.p_slope_contact_law import make_slope_contact_law
+
+
+E = 10000
+kappa = 0.4
+
+
+@dataclass
+class StaticSetup(StaticDisplacementProblem):
+    mu_coef: ... = E / (1 + kappa)
+    la_coef: ... = E * kappa / ((1 + kappa) * (1 - 2 * kappa))
+    contact_law: ... = make_slope_contact_law(slope=1)
+
+    @staticmethod
+    def inner_forces(x, t=None):
+        return np.array([0, 0])
+
+    @staticmethod
+    def outer_forces(x, t=None):
+        return np.array([0, -1]) if x[0] > 1.9 and x[1] < 0.1 else np.zeros(2)
+
+    @staticmethod
+    def friction_bound(u_nu: float) -> float:
+        return 0
+
+    boundaries: ... = BoundariesDescription(
+        contact=lambda x: x[1] < 0.1 and x[0] < 0.5, dirichlet=lambda x: x[0] == 0
+    )
+
+
+def main(config: Config):
+    """
+    Entrypoint to example.
+
+    To see result of simulation you need to call from python `main(Config().init())`.
+    """
+    mesh_prop = MeshProperties(mesh_type="msh_file", path="examples/example_mesh.msh")
+    setup = StaticSetup(mesh_prop=mesh_prop)
+    runner = StaticSolver(setup, "schur")
+
+    state = runner.solve(verbose=True, initial_displacement=setup.initial_displacement)
+    Drawer(state=state, config=config).draw(show=config.show, save=config.save)
+
+
+if __name__ == "__main__":
+    main(Config().init())

examples/example_piezo_quasistatic.py

@@ -10,6 +10,7 @@
 from conmech.simulations.problem_solver import PiezoelectricTimeDependentSolver
 from conmech.scenarios.problems import PiezoelectricQuasistaticProblem
 from conmech.plotting.drawer import Drawer
+from conmech.mesh.mesh import MeshProperties
 
 from examples.p_slope_contact_law import make_slope_contact_law
 
@@ -30,8 +31,6 @@ def h_temp(u_tau):  # potential  # TODO # 48
 
 @dataclass()
 class PQuasistaticSetup(PiezoelectricQuasistaticProblem):
-    grid_height: ... = 1.0
-    elements_number: ... = (2, 2)
     mu_coef: ... = 45
     la_coef: ... = 105
     th_coef: ... = 4.5
@@ -86,7 +85,8 @@ def main(config: Config):
 
     To see result of simulation you need to call from python `main(Config().init())`.
     """
-    setup = PQuasistaticSetup(mesh_type="Barboteu2008")
+    mesh_prop = MeshProperties(mesh_type="Barboteu2008", mesh_density=[2, 2], grid_height=1.0)
+    setup = PQuasistaticSetup(mesh_prop)
     runner = PiezoelectricTimeDependentSolver(setup, solving_method="global")
     steps = 100 if not config.test else 10
     output = steps // 5

examples/example_piezoelectric_dynamic.py

@@ -11,6 +11,7 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import PiezoelectricDynamicProblem
 from conmech.simulations.problem_solver import PiezoelectricTimeDependentSolver
+from conmech.mesh.mesh import MeshProperties
 from examples.p_slope_contact_law import make_slope_contact_law
 
 
@@ -46,8 +47,6 @@ def h_temp(u_tau):  # potential
 
 @dataclass()
 class PDynamicSetup(PiezoelectricDynamicProblem):
-    grid_height: ... = 1.0
-    elements_number: ... = (4, 3)
     mu_coef: ... = 45
     la_coef: ... = 105
     th_coef: ... = 4.5
@@ -102,7 +101,8 @@ def main(config: Config):
 
     To see result of simulation you need to call from python `main(Config().init())`.
     """
-    setup = PDynamicSetup(mesh_type="Barboteu2008")
+    mesh_prop = MeshProperties(mesh_type="Barboteu2008", mesh_density=[3, 4], grid_height=1.0)
+    setup = PDynamicSetup(mesh_prop)
     runner = PiezoelectricTimeDependentSolver(setup, solving_method="global")
 
     steps = 100 if not config.test else 10

examples/example_poisson.py

@@ -7,13 +7,11 @@
 from conmech.plotting.drawer import Drawer
 from conmech.scenarios.problems import PoissonProblem
 from conmech.simulations.problem_solver import PoissonSolver
+from conmech.mesh.mesh import MeshProperties
 
 
 @dataclass()
 class StaticPoissonSetup(PoissonProblem):
-    grid_height: ... = 1
-    elements_number: ... = (8, 8)
-
     @staticmethod
     def external_temperature(x: np.ndarray, t: Optional[float] = None) -> np.ndarray:
         if x[1] == 0:
@@ -46,7 +44,8 @@ def main(config: Config):
 
     To see result of simulation you need to call from python `main(Config().init())`.
     """
-    setup = StaticPoissonSetup(mesh_type="cross")
+    mesh_prop = MeshProperties(mesh_type="cross", mesh_density=[8, 8], grid_height=1.0)
+    setup = StaticPoissonSetup(mesh_prop)
     runner = PoissonSolver(setup, "direct")
 
     state = runner.solve(verbose=True)