RodMeshSurfaceContact class init

elastica/contact_forces.py

@@ -426,31 +426,51 @@ def _calculate_contact_forces_self_rod(
 
 
 @numba.njit(cache=True)
-def apply_normal_force_numba(
-    faces_normals,
-    faces_centers,
-    element_position,
-    position_idx_array,
+def _calculate_contact_forces_rod_mesh_surface(
+    faces_normals: np.ndarray,
+    faces_centers: np.ndarray,
+    element_position: np.ndarray,
+    position_idx_array: np.ndarray,
     face_idx_array,
-    surface_tol,
-    k,
-    nu,
-    radius,
-    mass,
-    position_collection,
-    velocity_collection,
-    internal_forces,
-    external_forces,
-    lengths,
-):
+    surface_tol: float,
+    k: float,
+    nu: float,
+    radius: np.array,
+    mass: np.array,
+    velocity_collection: np.ndarray,
+    external_forces: np.ndarray,
+) -> tuple:
     """
     This function computes the plane force response on the element, in the
     case of contact. Contact model given in Eqn 4.8 Gazzola et. al. RSoS 2018 paper
     is used.
 
     Parameters
     ----------
-    system_one
+    faces_normals: np.ndarray
+        mesh cell's normal vectors
+    faces_centers: np.ndarray
+        mesh cell's center points
+    element_position: np.ndarray
+        rod element's center points
+    position_idx_array: np.ndarray
+        rod element's index array
+    face_idx_array: np.ndarray
+        mesh cell's index array
+    surface_tol: float
+        Penetration tolerance between the surface and the rod-like object
+    k: float
+        Contact spring constant
+    nu: float
+        Contact damping constant
+    radius: np.array
+        rod element's radius
+    mass: np.array
+        rod element's mass
+    velocity_collection: np.ndarray
+        rod element's velocity
+    external_forces: np.ndarray
+        rod element's external forces
 
     Returns
     -------
@@ -464,22 +484,22 @@ def apply_normal_force_numba(
 
     if len(face_idx_array) > 0:
         element_position_contacts = element_position[:, position_idx_array]
-        contact_facet_centers = faces_centers[:, face_idx_array]
+        contact_face_centers = faces_centers[:, face_idx_array]
         normals_on_elements = faces_normals[:, face_idx_array]
         radius_contacts = radius[position_idx_array]
         element_velocity_contacts = element_velocity[:, position_idx_array]
 
     else:
         element_position_contacts = element_position
-        contact_facet_centers = np.zeros_like(element_position)
+        contact_face_centers = np.zeros_like(element_position)
         normals_on_elements = np.zeros_like(element_position)
         radius_contacts = radius
         element_velocity_contacts = element_velocity
 
     # Elastic force response due to penetration
 
     distance_from_plane = _batch_dot(
-        normals_on_elements, (element_position_contacts - contact_facet_centers)
+        normals_on_elements, (element_position_contacts - contact_face_centers)
     )
     plane_penetration = (
         -np.abs(np.minimum(distance_from_plane - radius_contacts, 0.0)) ** 1.5
@@ -825,7 +845,7 @@ def apply_contact(
         )
 
 
-class RodMeshSurfaceContact(NoContact):
+class RodMeshSurfaceContactWithGridMethod(NoContact):
     """
     This class is for applying contact forces between rod-mesh_surface.
     First system is always rod and second system is always mesh_surface.
@@ -835,17 +855,16 @@ class RodMeshSurfaceContact(NoContact):
     How to define contact between rod and mesh_surface.
 
     >>> simulator.detect_contact_between(rod, mesh_surface).using(
-    ...    RodMeshSurfaceContact,
+    ...    RodMeshSurfaceContactWithGridMethod,
     ...    k=1e4,
     ...    nu=10,
     ...    surface_tol=1e-2,
     ... )
-
-
-    .. [1] Preclik T., Popa Constantin., Rude U., Regularizing a Time-Stepping Method for Rigid Multibody Dynamics, Multibody Dynamics 2011, ECCOMAS. URL: https://www10.cs.fau.de/publications/papers/2011/Preclik_Multibody_Ext_Abstr.pdf
     """
 
-    def __init__(self, k: float, nu: float, surface_tol=1e-4):
+    def __init__(
+        self, k: float, nu: float, faces_grid: dict, grid_size: float, surface_tol=1e-4
+    ):
         """
 
         Parameters
@@ -854,24 +873,33 @@ def __init__(self, k: float, nu: float, surface_tol=1e-4):
             Stiffness coefficient between the plane and the rod-like object.
         nu: float
             Dissipation coefficient between the plane and the rod-like object.
+        faces_grid: dict
+            Dictionary containing the grid information of the mesh surface.
+        grid_size: float
+            Grid size of the mesh surface.
         surface_tol: float
             Penetration tolerance between the surface and the rod-like object.
 
         """
-        super(RodMeshSurfaceContact, self).__init__()
+        super(RodMeshSurfaceContactWithGridMethod, self).__init__()
         # n_faces = faces.shape[-1]
         self.k = k
         self.nu = nu
+        self.faces_grid = faces_grid
+        self.grid_size = grid_size
         self.surface_tol = surface_tol
 
     def _check_systems_validity(
         self,
         system_one: SystemType,
         system_two: AllowedContactType,
+        faces_grid: dict,
+        grid_size: float,
     ) -> None:
         """
         This checks the contact order and type of a SystemType object and an AllowedContactType object.
-        For the RodMeshSurfaceContact class first_system should be a rod and second_system should be a mesh_surface.
+        For the RodMeshSurfaceContact class first_system should be a rod and second_system should be a mesh_surface;
+        morever, the imported grid's attributes should match imported rod-mesh_surface(in contact) grid's attributes.
 
         Parameters
         ----------
@@ -891,7 +919,22 @@ def _check_systems_validity(
                 )
             )
 
-    def apply_normal_force(
+        elif not faces_grid["grid_size"] == grid_size:
+            raise TypeError(
+                "Imported grid size does not match with the current rod-mesh_surface grid size. "
+            )
+
+        elif not faces_grid["model_path"] == system_two.model_path:
+            raise TypeError(
+                "Imported grid's model path does not match with the current mesh_surface model path. "
+            )
+
+        elif not faces_grid["surface_reorient"] == system_two.mesh_orientation:
+            raise TypeError(
+                "Imported grid's surface orientation does not match with the current mesh_surface rientation. "
+            )
+
+    def apply_contact(
         self, system_one: RodType, system_two: AllowedContactType
     ) -> tuple:
         """
@@ -924,14 +967,14 @@ def apply_normal_force(
             self.face_idx_array,
             self.element_position,
         ) = find_contact_faces_idx(
-            self.facets_grid,
+            self.faces_grid,
             self.mesh_surface_x_min,
             self.mesh_surface_y_min,
             self.grid_size,
             system_one.position_collection,
         )
 
-        return apply_normal_force_numba(
+        return _calculate_contact_forces_rod_mesh_surface(
             self.mesh_surface_face_normals,
             self.mesh_surface_face_centers,
             self.element_position,
@@ -942,9 +985,6 @@ def apply_normal_force(
             self.nu,
             system_one.radius,
             system_one.mass,
-            system_one.position_collection,
             system_one.velocity_collection,
-            system_one.internal_forces,
             system_one.external_forces,
-            system_one.lengths,
         )

elastica/contact_utils.py

@@ -194,56 +194,85 @@ def find_contact_faces_idx(faces_grid, x_min, y_min, grid_size, position_collect
     element_position = node_to_element_position(position_collection)
     n_element = element_position.shape[-1]
     position_idx_array = np.empty((0))
-    facet_idx_array = np.empty((0))
+    face_idx_array = np.empty((0))
     grid_position = np.round(
         (element_position[0:2, :] - np.array([x_min, y_min]).reshape((2, 1)))
         / grid_size
     )
 
-    # find facet neighborhood of each element position
+    # find face neighborhood of each element position
 
     for i in range(n_element):
         try:
-            facet_idx_1 = faces_grid[
+            face_idx_1 = faces_grid[
                 (int(grid_position[0, i]), int(grid_position[1, i]))
             ]  # first quadrant
         except Exception:
-            facet_idx_1 = np.empty((0))
+            face_idx_1 = np.empty((0))
         try:
-            facet_idx_2 = faces_grid[
+            face_idx_2 = faces_grid[
                 (int(grid_position[0, i] - 1), int(grid_position[1, i]))
             ]  # second quadrant
         except Exception:
-            facet_idx_2 = np.empty((0))
+            face_idx_2 = np.empty((0))
         try:
-            facet_idx_3 = faces_grid[
+            face_idx_3 = faces_grid[
                 (int(grid_position[0, i] - 1), int(grid_position[1, i] - 1))
             ]  # third quadrant
         except Exception:
-            facet_idx_3 = np.empty((0))
+            face_idx_3 = np.empty((0))
         try:
-            facet_idx_4 = faces_grid[
+            face_idx_4 = faces_grid[
                 (int(grid_position[0, i]), int(grid_position[1, i] - 1))
             ]  # fourth quadrant
         except Exception:
-            facet_idx_4 = np.empty((0))
-        facet_idx_element = np.concatenate(
-            (facet_idx_1, facet_idx_2, facet_idx_3, facet_idx_4)
+            face_idx_4 = np.empty((0))
+        face_idx_element = np.concatenate(
+            (face_idx_1, face_idx_2, face_idx_3, face_idx_4)
         )
-        facet_idx_element_no_duplicates = np.unique(facet_idx_element)
-        if facet_idx_element_no_duplicates.size == 0:
+        face_idx_element_no_duplicates = np.unique(face_idx_element)
+        if face_idx_element_no_duplicates.size == 0:
             raise RuntimeError(
-                "Snake outside surface boundary"
-            )  # a snake element is on four grids with no facets
+                "Rod object out of grid bounds"
+            )  # a snake element is on four grids with no faces
 
-        facet_idx_array = np.concatenate(
-            (facet_idx_array, facet_idx_element_no_duplicates)
+        face_idx_array = np.concatenate(
+            (face_idx_array, face_idx_element_no_duplicates)
         )
-        n_contacts = facet_idx_element_no_duplicates.shape[0]
+        n_contacts = face_idx_element_no_duplicates.shape[0]
         position_idx_array = np.concatenate(
             (position_idx_array, i * np.ones((n_contacts,)))
         )
 
     position_idx_array = position_idx_array.astype(int)
-    facet_idx_array = facet_idx_array.astype(int)
-    return position_idx_array, facet_idx_array, element_position
+    face_idx_array = face_idx_array.astype(int)
+    return position_idx_array, face_idx_array, element_position
+
+
+"""
+@numba.njit(cache=True)
+def surface_grid_numba(faces, grid_size, face_x_left, face_x_right, face_y_down, face_y_up):
+    x_min = np.min(faces[0, :, :])
+    y_min = np.min(faces[1, :, :])
+    n_x_positions = int(np.ceil((np.max(faces[0, :, :]) - x_min) / grid_size))
+    n_y_positions = int(np.ceil((np.max(faces[1, :, :]) - y_min) / grid_size))
+    faces_grid = dict()
+    for i in range(n_x_positions):
+        x_left = x_min + (i * grid_size)
+        x_right = x_min + ((i + 1) * grid_size)
+        for j in range(n_y_positions):
+            y_down = y_min + (j * grid_size)
+            y_up = y_min + ((j + 1) * grid_size)
+            if np.any(np.where(((face_y_down > y_up) + (face_y_up < y_down) + (face_x_right < x_left) + (face_x_left > x_right)) == 0)[0]):
+                faces_grid[(i, j)] = np.where(((face_y_down > y_up) + (face_y_up < y_down) + (face_x_right < x_left) + (face_x_left > x_right)) == 0)[0]
+    return faces_grid
+
+
+def surface_grid(faces, grid_size):
+    face_x_left = np.min(faces[0, :, :], axis=0)
+    face_y_down = np.min(faces[1, :, :], axis=0)
+    face_x_right = np.max(faces[0, :, :], axis=0)
+    face_y_up = np.max(faces[1, :, :], axis=0)
+
+    return dict(surface_grid_numba(faces, grid_size, face_x_left, face_x_right, face_y_down, face_y_up))
+"""

elastica/mesh/mesh_initializer.py

@@ -2,6 +2,7 @@
 
 import pyvista as pv
 import numpy as np
+import os
 
 
 class Mesh:
@@ -47,6 +48,9 @@ class Mesh:
         - Dimension: (3 spatial coordinates, 3 orthonormal basis vectors)
         - Initial value: [[1,0,0],[0,1,0],[0,0,1]]
 
+    mesh.model_path:
+            - Stores the path to the mesh file.
+
     Methods:
     --------
 
@@ -78,9 +82,16 @@ class Mesh:
     """
 
     def __init__(self, filepath: str) -> None:
-        self.mesh = pv.read(filepath)
-        self.orientation_cube = pv.Cube()
-        self.mesh_update()
+        if os.path.isfile(filepath):
+            self.mesh = pv.read(filepath)
+            self.model_path = filepath
+            self.orientation_cube = pv.Cube()
+            self.mesh_update()
+        else:
+            raise FileNotFoundError(
+                "Please check the filepath.\n"
+                " Please be sure to add .stl / .obj at the end of the filepath, if already present, ignore"
+            )
 
     def mesh_update(self) -> None:
         """

elastica/surface/__init__.py

@@ -3,4 +3,3 @@
 from elastica.surface.surface_base import SurfaceBase
 from elastica.surface.mesh_surface import MeshSurface
 from elastica.surface.plane import Plane
-