Merge pull request #120 from wprzadka/inherit-mesh

make Mesh inherit RawMesh

conmech/dynamics/dynamics.py

@@ -53,7 +53,7 @@ def reinitialize_matrices(self, elements_density: Optional[np.ndarray] = None):
             self._w_matrix,
             self._local_stifness_matrices,
         ) = get_basic_matrices(
-            elements=self.body.mesh.elements, nodes=self.body.mesh.initial_nodes
+            elements=self.body.mesh.elements, nodes=self.body.mesh.nodes
         )  # + self.displacement_old)
 
         if elements_density is not None:

conmech/mesh/mesh.py

@@ -7,6 +7,7 @@
 from conmech.mesh.boundaries_factory import BoundariesFactory
 from conmech.mesh.boundaries import Boundaries
 from conmech.properties.mesh_description import MeshDescription
+from conmech.mesh.zoo.raw_mesh import RawMesh
 
 
 @numba.njit
@@ -82,48 +83,46 @@ def get_base_seed_indices_numba(nodes):
     return base_seed_indices, int(np.argmin(errors))
 
 
-# TODO make it inherit from RawMesh
-class Mesh:
+class Mesh(RawMesh):
     def __init__(
         self,
         mesh_descr: MeshDescription,
         boundaries_description: BoundariesDescription,
     ):
-        self.initial_nodes: np.ndarray
-        self.elements: np.ndarray
         self.edges: np.ndarray
 
         self.boundaries: Boundaries
 
         input_nodes, input_elements = mesh_builders.build_mesh(mesh_descr=mesh_descr)
+        super().__init__(input_nodes, input_elements)
         unordered_nodes, unordered_elements = remove_unconnected_nodes_numba(
             input_nodes, input_elements
         )
         (
-            self.initial_nodes,
+            self.nodes,
             self.elements,
             self.boundaries,
         ) = BoundariesFactory.identify_boundaries_and_reorder_nodes(
             unordered_nodes, unordered_elements, boundaries_description
         )
-        edges_matrix = get_edges_matrix(nodes_count=len(self.initial_nodes), elements=self.elements)
+        edges_matrix = get_edges_matrix(nodes_count=len(self.nodes), elements=self.elements)
         self.edges = get_edges_list_numba(edges_matrix)
 
     @property
     def dimension(self):
-        return self.initial_nodes.shape[1]
+        return self.nodes.shape[1]
 
     @property
     def scale(self):
-        return np.max(self.initial_nodes, axis=0) - np.min(self.initial_nodes, axis=0)
+        return np.max(self.nodes, axis=0) - np.min(self.nodes, axis=0)
 
     def normalize_shift(self, vectors):
         _ = self
         return vectors - np.mean(vectors, axis=0)
 
     @property
     def normalized_initial_nodes(self):
-        return self.normalize_shift(self.initial_nodes)
+        return self.normalize_shift(self.nodes)
 
     @property
     def input_initial_nodes(self):
@@ -180,7 +179,7 @@ def boundary_indices(self):
 
     @property
     def initial_boundary_nodes(self):
-        return self.initial_nodes[self.boundary_indices]
+        return self.nodes[self.boundary_indices]
 
     @property
     def contact_indices(self):
@@ -207,7 +206,7 @@ def free_indices(self):
 
     @property
     def nodes_count(self):
-        return len(self.initial_nodes)
+        return len(self.nodes)
 
     @property
     def boundary_surfaces_count(self):

conmech/plotting/drawer.py

@@ -24,7 +24,7 @@ def __init__(self, state, config: Config):
         self.state = state
         self.config = config
         self.mesh = state.body.mesh
-        self.node_size = 2 + (300 / len(self.mesh.initial_nodes))
+        self.node_size = 2 + (300 / len(self.mesh.nodes))
         self.line_width = self.node_size / 2
         self.deformed_mesh_color = "k"
         self.original_mesh_color = "0.7"
@@ -90,22 +90,22 @@ def set_axes_limits(self, axes, foundation):
         # pylint: disable=nested-min-max
         if self.x_min is None:
             self.x_min = min(
-                min(self.state.body.mesh.initial_nodes[:, 0]), min(self.state.displaced_nodes[:, 0])
+                min(self.state.body.mesh.nodes[:, 0]), min(self.state.displaced_nodes[:, 0])
             )
         if self.x_max is None:
             self.x_max = max(
-                max(self.state.body.mesh.initial_nodes[:, 0]), max(self.state.displaced_nodes[:, 0])
+                max(self.state.body.mesh.nodes[:, 0]), max(self.state.displaced_nodes[:, 0])
             )
         dx = self.x_max - self.x_min
         x_margin = dx * 0.2
         xlim = (self.x_min - x_margin, self.x_max + x_margin)
         if self.y_min is None:
             self.y_min = min(
-                min(self.state.body.mesh.initial_nodes[:, 1]), min(self.state.displaced_nodes[:, 1])
+                min(self.state.body.mesh.nodes[:, 1]), min(self.state.displaced_nodes[:, 1])
             )
         if self.y_max is None:
             self.y_max = max(
-                max(self.state.body.mesh.initial_nodes[:, 1]), max(self.state.displaced_nodes[:, 1])
+                max(self.state.body.mesh.nodes[:, 1]), max(self.state.displaced_nodes[:, 1])
             )
         dy = self.y_max - self.y_min
         y_margin = dy * 0.2
@@ -123,7 +123,7 @@ def set_axes_limits(self, axes, foundation):
     def draw_meshes(self, axes):
         if self.original_mesh_color is not None:
             self.draw_mesh(
-                self.mesh.initial_nodes,
+                self.mesh.nodes,
                 axes,
                 label="Original",
                 node_color=self.original_mesh_color,
@@ -200,7 +200,7 @@ def draw_forces(self, axes):
             neumann_nodes = list(set(neumann_nodes.flatten()))
             x = self.state.displaced_nodes[neumann_nodes]
             v = self.state.body.dynamics.force.outer.node_source(
-                self.state.body.mesh.initial_nodes, self.state.time
+                self.state.body.mesh.nodes, self.state.time
             )[neumann_nodes]
 
             scale = self.outer_forces_scale

conmech/scene/body_forces.py

@@ -66,17 +66,17 @@ def clear(self):
         self.outer.cache = None
 
     def get_integrated_inner_forces(self, time: float):
-        inner_forces = self.inner.node_source(self.body.mesh.initial_nodes, time)
+        inner_forces = self.inner.node_source(self.body.mesh.nodes, time)
         # TODO: should be only on boundary!
         return self.body.dynamics.volume_at_nodes @ inner_forces
 
     def get_integrated_outer_forces(self, time: float):
         neumann_surfaces = get_surface_per_boundary_node_numba(
             boundary_surfaces=self.body.mesh.neumann_boundary,
             considered_nodes_count=self.body.mesh.nodes_count,
-            moved_nodes=self.body.mesh.initial_nodes,
+            moved_nodes=self.body.mesh.nodes,
         )
-        outer_forces = self.outer.node_source(self.body.mesh.initial_nodes, time)
+        outer_forces = self.outer.node_source(self.body.mesh.nodes, time)
         return neumann_surfaces * outer_forces
 
     def get_integrated_field_sources_column(self, time: float):

conmech/simulations/problem_solver.py

@@ -375,7 +375,7 @@ def solve(self, *, initial_displacement: Callable, verbose: bool = False, **kwar
         """
         state = State(self.body)
         state.displacement = initial_displacement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
 
         self.step_solver.u_vector[:] = state.displacement.T.ravel().copy()
@@ -436,11 +436,9 @@ def solve(
         output_step = (0, *output_step) if output_step else (0, n_steps)  # 0 for diff
 
         state = State(self.body)
-        state.absement[:] = initial_absement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
-        )
+        state.absement[:] = initial_absement(self.body.mesh.nodes[: self.body.mesh.nodes_count])
         state.displacement[:] = initial_displacement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
 
         self.step_solver.b_vector[:] = state.absement.T.ravel().copy()
@@ -505,11 +503,9 @@ def solve(
 
         state: State = State(self.body)
         state.displacement[:] = initial_displacement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
-        )
-        state.velocity[:] = initial_velocity(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
+        state.velocity[:] = initial_velocity(self.body.mesh.nodes[: self.body.mesh.nodes_count])
 
         self.step_solver.u_vector[:] = state.displacement.T.ravel().copy()
         self.step_solver.v_vector[:] = state.velocity.T.ravel().copy()
@@ -576,13 +572,11 @@ def solve(
 
         state = TemperatureState(self.body)
         state.displacement[:] = initial_displacement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
-        )
-        state.velocity[:] = initial_velocity(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
+        state.velocity[:] = initial_velocity(self.body.mesh.nodes[: self.body.mesh.nodes_count])
         state.temperature[:] = initial_temperature(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
 
         solution = state.velocity.reshape(2, -1)
@@ -671,13 +665,11 @@ def solve(
 
         state = PiezoelectricState(self.body)
         state.displacement[:] = initial_displacement(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
-        )
-        state.velocity[:] = initial_velocity(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
+        state.velocity[:] = initial_velocity(self.body.mesh.nodes[: self.body.mesh.nodes_count])
         state.electric_potential[:] = initial_electric_potential(
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count]
         )
 
         solution = state.velocity.reshape(2, -1)

conmech/solvers/direct.py

@@ -57,7 +57,7 @@ def _solve_impl(self, initial_guess: np.ndarray, **kwargs) -> np.ndarray:
                 self.equation,
                 initial_guess,
                 args=(
-                    self.body.mesh.initial_nodes,
+                    self.body.mesh.nodes,
                     self.body.mesh.contact_boundary,
                     self.body.mesh.boundaries.contact_normals,
                     self.node_relations,

conmech/solvers/optimization/optimization.py

@@ -106,7 +106,7 @@ def _solve_impl(
         maxiter = kwargs.get("maxiter", int(len(initial_guess) * 1e9))
         tol = kwargs.get("tol", 1e-12)
         args = (
-            self.body.mesh.initial_nodes,
+            self.body.mesh.nodes,
             self.body.mesh.contact_boundary,
             self.body.mesh.boundaries.contact_normals,
             self.lhs,

conmech/solvers/validator.py

@@ -27,7 +27,7 @@ def validate(self, state: State, solution: np.ndarray) -> float:
         quality_inv = np.linalg.norm(
             self.rhs(
                 solution,
-                state.body.mesh.initial_nodes,
+                state.body.mesh.nodes,
                 state.body.mesh.contact_boundary,
                 state.body.mesh.boundaries.contact_normals,
                 self.elasticity,

conmech/state/state.py

@@ -14,7 +14,7 @@ def __init__(self, body):
         self.displacement: np.ndarray = np.zeros(
             (self.body.mesh.nodes_count, self.body.mesh.dimension)
         )
-        self.displaced_nodes: np.ndarray = np.copy(self.body.mesh.initial_nodes)
+        self.displaced_nodes: np.ndarray = np.copy(self.body.mesh.nodes)
         self.velocity: np.ndarray = np.zeros((self.body.mesh.nodes_count, self.body.mesh.dimension))
         self.setup = None
         self.__stress: np.ndarray = None
@@ -26,7 +26,7 @@ def set_displacement(
     ):
         self.displacement = displacement_vector.reshape((self.body.mesh.dimension, -1)).T
         self.displaced_nodes[: self.body.mesh.nodes_count, :] = (
-            self.body.mesh.initial_nodes[: self.body.mesh.nodes_count, :] + self.displacement[:, :]
+            self.body.mesh.nodes[: self.body.mesh.nodes_count, :] + self.displacement[:, :]
         )
         if update_absement:
             dt = time - self.time
@@ -39,8 +39,7 @@ def set_velocity(self, velocity_vector: np.ndarray, time: float, *, update_displ
             dt = time - self.time
             self.displacement += dt * self.velocity
             self.displaced_nodes[: self.body.mesh.nodes_count, :] = (
-                self.body.mesh.initial_nodes[: self.body.mesh.nodes_count, :]
-                + self.displacement[:, :]
+                self.body.mesh.nodes[: self.body.mesh.nodes_count, :] + self.displacement[:, :]
             )
         self.time = time
 
@@ -54,7 +53,7 @@ def stress(self):
                 absement=self.absement,
                 setup=self.setup,
                 elements=self.body.mesh.elements,
-                nodes=self.body.mesh.initial_nodes,
+                nodes=self.body.mesh.nodes,
                 time=self.time,
             )
         return self.__stress

examples/error_estimates.py

@@ -15,18 +15,18 @@
 def compare(ref: TemperatureState, sol: TemperatureState):
     ut = 0
     tt = 0
-    x = sol.body.mesh.initial_nodes[:, 0]
-    y = sol.body.mesh.initial_nodes[:, 1]
+    x = sol.body.mesh.nodes[:, 0]
+    y = sol.body.mesh.nodes[:, 1]
 
     soltri = tri.Triangulation(x, y, triangles=sol.body.mesh.elements)
     u1hi = tri.LinearTriInterpolator(soltri, sol.displacement[:, 0])
     u2hi = tri.LinearTriInterpolator(soltri, sol.displacement[:, 1])
     thi = tri.LinearTriInterpolator(soltri, sol.temperature)
 
     for element in ref.body.mesh.elements:
-        x0 = ref.body.mesh.initial_nodes[element[0]]
-        x1 = ref.body.mesh.initial_nodes[element[1]]
-        x2 = ref.body.mesh.initial_nodes[element[2]]
+        x0 = ref.body.mesh.nodes[element[0]]
+        x1 = ref.body.mesh.nodes[element[1]]
+        x2 = ref.body.mesh.nodes[element[2]]
         u1_0 = ref.displacement[element[0], 0]
         u1_1 = ref.displacement[element[1], 0]
         u1_2 = ref.displacement[element[2], 0]

examples/example_nonhomogenous_density.py

@@ -51,7 +51,7 @@ def main(config: Config):
 
     elem_centers = np.empty(shape=(len(runner.body.mesh.elements), 2))
     for idx, elem in enumerate(runner.body.mesh.elements):
-        verts = runner.body.mesh.initial_nodes[elem]
+        verts = runner.body.mesh.nodes[elem]
         elem_centers[idx] = np.sum(verts, axis=0) / len(elem)
     density = np.asarray([1 if x[0] < 1 else 0.2 for x in elem_centers])
 

examples/example_static.py

@@ -67,7 +67,7 @@ def main(config: Config):
         fig = plt.figure()
         axs = fig.add_subplot(111, projection="3d")
         # Draw nodes
-        nodes = state.body.mesh.initial_nodes
+        nodes = state.body.mesh.nodes
         axs.scatter(nodes[:, 0], nodes[:, 1], nodes[:, 2], c="b", marker="o")
 
         # Draw elements

examples/example_temperature_dynamic_2.py

@@ -24,8 +24,8 @@
 
 
 def compute_error(ref: TemperatureState, sol: TemperatureState):
-    x = sol.body.mesh.initial_nodes[:, 0]
-    y = sol.body.mesh.initial_nodes[:, 1]
+    x = sol.body.mesh.nodes[:, 0]
+    y = sol.body.mesh.nodes[:, 1]
 
     soltri = tri.Triangulation(x, y, triangles=sol.body.mesh.elements)
     u1hi = tri.LinearTriInterpolator(soltri, sol.velocity[:, 0])
@@ -35,9 +35,9 @@ def compute_error(ref: TemperatureState, sol: TemperatureState):
     total_abs_error = np.full_like(ref.temperature, fill_value=np.nan)
 
     for element in ref.body.mesh.elements:
-        x0 = ref.body.mesh.initial_nodes[element[0]]
-        x1 = ref.body.mesh.initial_nodes[element[1]]
-        x2 = ref.body.mesh.initial_nodes[element[2]]
+        x0 = ref.body.mesh.nodes[element[0]]
+        x1 = ref.body.mesh.nodes[element[1]]
+        x2 = ref.body.mesh.nodes[element[2]]
         if total_abs_error[element[0]] != np.nan:
             total_abs_error[element[0]] = (  # abs(ref.velocity[element[0], 0] - u1hi(*x0))
                 # + abs(ref.velocity[element[0], 1] - u2hi(*x0))

tests/test_conmech/regression/test_Jureczka_and_Ochal_2019.py

@@ -72,8 +72,8 @@ def test(solving_method):
         fixed_point_abs_tol=0.001, initial_displacement=setup.initial_displacement
     )
 
-    displacement = result.body.mesh.initial_nodes[:] - result.displaced_nodes[:]
-    std_ids = standard_boundary_nodes(runner.body.mesh.initial_nodes, runner.body.mesh.elements)
+    displacement = result.body.mesh.nodes[:] - result.displaced_nodes[:]
+    std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements)
 
     # print result
     np.set_printoptions(precision=8, suppress=True)

tests/test_conmech/regression/test_density_update.py

@@ -22,7 +22,7 @@ def solving_method(request):
 def get_elem_centers(runner: NonHomogenousSolver):
     elem_centers = np.empty(shape=(len(runner.body.mesh.elements), 2))
     for idx, elem in enumerate(runner.body.mesh.elements):
-        verts = runner.body.mesh.initial_nodes[elem]
+        verts = runner.body.mesh.nodes[elem]
         elem_centers[idx] = np.sum(verts, axis=0) / len(elem)
     return elem_centers
 
@@ -252,8 +252,8 @@ def test_nonhomogenous_solver(solving_method, setup, density_func, expected_disp
 
     result = runner.solve(initial_displacement=setup.initial_displacement)
 
-    displacement = result.displaced_nodes[:] - result.body.mesh.initial_nodes[:]
-    std_ids = standard_boundary_nodes(runner.body.mesh.initial_nodes, runner.body.mesh.elements)
+    displacement = result.displaced_nodes[:] - result.body.mesh.nodes[:]
+    std_ids = standard_boundary_nodes(runner.body.mesh.nodes, runner.body.mesh.elements)
 
     # print result
     np.set_printoptions(precision=8, suppress=True)