From 2725342bf69fffa4bb928bb147e3925212ac143d Mon Sep 17 00:00:00 2001
From: Rahul-JOON <rahuljoon16@gmail.com>
Date: Thu, 24 Aug 2023 21:20:10 +0530
Subject: [PATCH 1/5] RodMeshSurfaceContact class init

---
 elastica/contact_forces.py               | 227 ++++++++++++++++++++++-
 elastica/contact_utils.py                |  61 ++++++
 elastica/surface/__init__.py             |   2 +
 elastica/surface/mesh_surface.py         |  84 +++++++++
 examples/Binder/1_Timoshenko_Beam.ipynb  |   1 -
 examples/Binder/2_Slithering_Snake.ipynb |  10 +-
 tests/test_surface/test_mesh_surface.py  |  43 +++++
 7 files changed, 423 insertions(+), 5 deletions(-)
 create mode 100644 elastica/surface/mesh_surface.py
 create mode 100644 tests/test_surface/test_mesh_surface.py

diff --git a/elastica/contact_forces.py b/elastica/contact_forces.py
index 9d377373..3df65e83 100644
--- a/elastica/contact_forces.py
+++ b/elastica/contact_forces.py
@@ -3,15 +3,17 @@
 from elastica.typing import RodType, SystemType, AllowedContactType
 from elastica.rod import RodBase
 from elastica.rigidbody import Cylinder
+from elastica.surface import MeshSurface
 from elastica.contact_utils import (
     _dot_product,
     _norm,
     _find_min_dist,
     _prune_using_aabbs_rod_cylinder,
     _prune_using_aabbs_rod_rod,
+    find_contact_faces_idx,
 )
-
-from elastica._linalg import _batch_product_k_ik_to_ik
+from elastica.interaction import node_to_element_velocity, elements_to_nodes_inplace
+from elastica._linalg import _batch_product_k_ik_to_ik, _batch_dot, _batch_norm
 import numba
 import numpy as np
 
@@ -423,6 +425,102 @@ def _calculate_contact_forces_self_rod(
                 external_forces_rod[..., j + 1] += net_contact_force
 
 
+@numba.njit(cache=True)
+def apply_normal_force_numba(
+    faces_normals,
+    faces_centers,
+    element_position,
+    position_idx_array,
+    face_idx_array,
+    surface_tol,
+    k,
+    nu,
+    radius,
+    mass,
+    position_collection,
+    velocity_collection,
+    internal_forces,
+    external_forces,
+    lengths,
+):
+    """
+    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
+
+    Returns
+    -------
+    magnitude of the plane response
+    """
+
+    # Damping force response due to velocity towards the plane
+    element_velocity = node_to_element_velocity(
+        mass=mass, node_velocity_collection=velocity_collection
+    )
+
+    if len(face_idx_array) > 0:
+        element_position_contacts = element_position[:, position_idx_array]
+        contact_facet_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)
+        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)
+    )
+    plane_penetration = (
+        -np.abs(np.minimum(distance_from_plane - radius_contacts, 0.0)) ** 1.5
+    )
+    elastic_force = -k * _batch_product_k_ik_to_ik(
+        plane_penetration, normals_on_elements
+    )
+
+    normal_component_of_element_velocity = _batch_dot(
+        normals_on_elements, element_velocity_contacts
+    )
+    damping_force = -nu * _batch_product_k_ik_to_ik(
+        normal_component_of_element_velocity, normals_on_elements
+    )
+
+    # Compute total plane response force
+    plane_response_force_contacts = elastic_force + damping_force
+
+    # Check if the rod elements are in contact with plane.
+    no_contact_point_idx = np.where(
+        (distance_from_plane - radius_contacts) > surface_tol
+    )[0]
+    # If rod element does not have any contact with plane, plane cannot apply response
+    # force on the element. Thus lets set plane response force to 0.0 for the no contact points.
+    plane_response_force_contacts[..., no_contact_point_idx] = 0.0
+
+    plane_response_forces = np.zeros_like(external_forces)
+    for i in range(len(position_idx_array)):
+        plane_response_forces[
+            :, position_idx_array[i]
+        ] += plane_response_force_contacts[:, i]
+
+    # Update the external forces
+    elements_to_nodes_inplace(plane_response_forces, external_forces)
+    return (
+        _batch_norm(plane_response_force_contacts),
+        no_contact_point_idx,
+        normals_on_elements,
+    )
+
+
 class RodRodContact(NoContact):
     """
     This class is for applying contact forces between rod-rod.
@@ -725,3 +823,128 @@ def apply_contact(
             self.k,
             self.nu,
         )
+
+
+class RodMeshSurfaceContact(NoContact):
+    """
+    This class is for applying contact forces between rod-mesh_surface.
+    First system is always rod and second system is always mesh_surface.
+
+    Examples
+    --------
+    How to define contact between rod and mesh_surface.
+
+    >>> simulator.detect_contact_between(rod, mesh_surface).using(
+    ...    RodMeshSurfaceContact,
+    ...    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):
+        """
+
+        Parameters
+        ----------
+        k: float
+            Stiffness coefficient between the plane and the rod-like object.
+        nu: float
+            Dissipation coefficient between the plane and the rod-like object.
+        surface_tol: float
+            Penetration tolerance between the surface and the rod-like object.
+
+        """
+        super(RodMeshSurfaceContact, self).__init__()
+        # n_faces = faces.shape[-1]
+        self.k = k
+        self.nu = nu
+        self.surface_tol = surface_tol
+
+    def _check_systems_validity(
+        self,
+        system_one: SystemType,
+        system_two: AllowedContactType,
+    ) -> 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.
+
+        Parameters
+        ----------
+        system_one
+            SystemType
+        system_two
+            AllowedContactType
+        """
+        if not issubclass(system_one.__class__, RodBase) or not issubclass(
+            system_two.__class__, MeshSurface
+        ):
+            raise TypeError(
+                "Systems provided to the contact class have incorrect order/type. \n"
+                " First system is {0} and second system is {1}. \n"
+                " First system should be a rod, second should be a mesh surface".format(
+                    system_one.__class__, system_two.__class__
+                )
+            )
+
+    def apply_normal_force(
+        self, system_one: RodType, system_two: AllowedContactType
+    ) -> tuple[np.array, np.array]:
+        """
+        In the case of contact with the plane, this function computes the plane reaction force on the element.
+
+        Parameters
+        ----------
+        system_one: object
+            Rod-like object.
+        system_two: Surface
+            Mesh surface.
+
+        Returns
+        -------
+        plane_response_force_mag : numpy.ndarray
+            1D (blocksize) array containing data with 'float' type.
+            Magnitude of plane response force acting on rod-like object.
+        no_contact_point_idx : numpy.ndarray
+            1D (blocksize) array containing data with 'int' type.
+            Index of rod-like object elements that are not in contact with the plane.
+        """
+
+        self.mesh_surface_faces = system_two.faces
+        self.mesh_surface_x_min = np.min(self.mesh_surface_faces[0, :, :])
+        self.mesh_surface_y_min = np.min(self.mesh_surface_faces[1, :, :])
+        self.mesh_surface_face_normals = system_two.face_normals
+        self.mesh_surface_face_centers = system_two.face_centers
+        (
+            self.position_idx_array,
+            self.face_idx_array,
+            self.element_position,
+        ) = find_contact_faces_idx(
+            self.facets_grid,
+            self.mesh_surface_x_min,
+            self.mesh_surface_y_min,
+            self.grid_size,
+            system_one.position_collection,
+        )
+
+        return apply_normal_force_numba(
+            self.mesh_surface_face_normals,
+            self.mesh_surface_face_centers,
+            self.element_position,
+            self.position_idx_array,
+            self.face_idx_array,
+            self.surface_tol,
+            self.k,
+            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,
+        )
diff --git a/elastica/contact_utils.py b/elastica/contact_utils.py
index 50159eab..a97f8f51 100644
--- a/elastica/contact_utils.py
+++ b/elastica/contact_utils.py
@@ -3,6 +3,7 @@
 from math import sqrt
 import numba
 import numpy as np
+from elastica.interaction import node_to_element_position
 
 
 @numba.njit(cache=True)
@@ -186,3 +187,63 @@ def _prune_using_aabbs_rod_rod(
         )
 
     return _aabbs_not_intersecting(aabb_rod_two, aabb_rod_one)
+
+
+def find_contact_faces_idx(faces_grid, x_min, y_min, grid_size, position_collection):
+
+    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))
+    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
+
+    for i in range(n_element):
+        try:
+            facet_idx_1 = faces_grid[
+                (int(grid_position[0, i]), int(grid_position[1, i]))
+            ]  # first quadrant
+        except Exception:
+            facet_idx_1 = np.empty((0))
+        try:
+            facet_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))
+        try:
+            facet_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))
+        try:
+            facet_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)
+        )
+        facet_idx_element_no_duplicates = np.unique(facet_idx_element)
+        if facet_idx_element_no_duplicates.size == 0:
+            raise RuntimeError(
+                "Snake outside surface boundary"
+            )  # a snake element is on four grids with no facets
+
+        facet_idx_array = np.concatenate(
+            (facet_idx_array, facet_idx_element_no_duplicates)
+        )
+        n_contacts = facet_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
diff --git a/elastica/surface/__init__.py b/elastica/surface/__init__.py
index ce755872..e33e254c 100644
--- a/elastica/surface/__init__.py
+++ b/elastica/surface/__init__.py
@@ -1,2 +1,4 @@
 __doc__ = """Surface classes"""
+
 from elastica.surface.surface_base import SurfaceBase
+from elastica.surface.mesh_surface import MeshSurface
diff --git a/elastica/surface/mesh_surface.py b/elastica/surface/mesh_surface.py
new file mode 100644
index 00000000..72c6dea7
--- /dev/null
+++ b/elastica/surface/mesh_surface.py
@@ -0,0 +1,84 @@
+__doc__ = """mesh surface class"""
+
+from elastica.surface.surface_base import SurfaceBase
+from elastica.mesh.mesh_initializer import Mesh
+import os
+
+
+class MeshSurface(SurfaceBase):
+    def __init__(self, filepath: str) -> None:
+        """
+        Mesh surface initializer.
+
+        Parameters
+        ----------
+        mesh file path (stl): str
+            - path to the mesh file
+
+        Attributes:
+        -----------
+
+        mesh_surface.faces:
+            - Stores the coordinates of the 3 vertices of each of the n faces of the imported mesh.
+            - Dimension: (3 spatial coordinates, 3 vertices, n faces)
+
+        mesh_surface.face_normals:
+            - Stores the coordinates of the unit normal vector of each of the n faces.
+            - Dimension: (3 spatial coordinates, n faces)
+
+        mesh_surface.face_centers:
+            - Stores the coordinates of the position vector of each of the n face centers.
+            - Dimension: (3 spatial coordinates, n faces)
+
+        mesh_surface.mesh_scale:
+            - Stores the 3 dimensions of the smallest box that could fit the mesh.
+            - Dimension: (3 spatial lengths)
+
+        mesh_surface.mesh_center:
+            - Stores the coordinates of the position vector of the center of the smallest box that could fit the mesh.
+            - Dimension: (3 spatial coordinates)
+
+        mesh_surface.mesh_orientation:
+            - store the 3 orthonormal basis vectors that define the mesh orientation.
+            - Dimension: (3 spatial coordinates, 3 orthonormal basis vectors)
+            - Initial value: [[1,0,0],[0,1,0],[0,0,1]]
+
+        Methods:
+        --------
+
+        mesh_surface.mesh_update():
+        Parameters: None
+            - This method updates/refreshes the mesh attributes in pyelastica geometry.
+            - By default this method is called at initialization and after every method that might change the mesh attributes.
+
+        mesh_surface.visualize():
+        Parameters: None
+            - This method visualizes the mesh using pyvista.
+
+        mesh_surface.translate():
+        Parameters: {numpy.ndarray-(3 spatial coordinates)}
+        ex : mesh.translate(np.array([1,1,1]))
+            - This method translates the mesh by a given vector.
+            - By default, the mesh's center is at the origin;
+            by calling this method, the mesh's center is translated to the given vector.
+
+        mesh_surface.scale():
+        Parameters: {numpy.ndarray-(3 spatial constants)}
+        ex : mesh.scale(np.array([1,1,1]))
+            - This method scales the mesh by a given factor in respective axes.
+
+        mesh_surface.rotate():
+        Parameters: {rotation_axis: unit vector[numpy.ndarray-(3 spatial coordinates)], angle: in degrees[float]}
+        ex : mesh.rotate(np.array([1,0,0]), 90)
+            - This method rotates the mesh by a given angle about a given axis.
+        """
+
+        if os.path.isfile(filepath):
+            mesh_surface = Mesh(filepath)
+            mesh_surface.mesh_update()  # using a mesh method so flake doesn't throw the error of unused variable
+
+        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"
+            )
diff --git a/examples/Binder/1_Timoshenko_Beam.ipynb b/examples/Binder/1_Timoshenko_Beam.ipynb
index 3176827b..158a2865 100644
--- a/examples/Binder/1_Timoshenko_Beam.ipynb
+++ b/examples/Binder/1_Timoshenko_Beam.ipynb
@@ -194,7 +194,6 @@
    },
    "outputs": [],
    "source": [
-    "\n",
     "dl = base_length / n_elem\n",
     "dt = 0.01 * dl\n",
     "timoshenko_sim.dampen(shearable_rod).using(\n",
diff --git a/examples/Binder/2_Slithering_Snake.ipynb b/examples/Binder/2_Slithering_Snake.ipynb
index fd41d9c7..23b0a94e 100644
--- a/examples/Binder/2_Slithering_Snake.ipynb
+++ b/examples/Binder/2_Slithering_Snake.ipynb
@@ -32,7 +32,13 @@
     "import numpy as np\n",
     "\n",
     "# import modules\n",
-    "from elastica.modules import BaseSystemCollection, Constraints, Forcing, CallBacks, Damping\n",
+    "from elastica.modules import (\n",
+    "    BaseSystemCollection,\n",
+    "    Constraints,\n",
+    "    Forcing,\n",
+    "    CallBacks,\n",
+    "    Damping,\n",
+    ")\n",
     "\n",
     "# import rod class, damping and forces to be applied\n",
     "from elastica.rod.cosserat_rod import CosseratRod\n",
@@ -105,7 +111,7 @@
     ")\n",
     "\n",
     "# Add rod to the snake system\n",
-    "snake_sim.append(shearable_rod)\n"
+    "snake_sim.append(shearable_rod)"
    ]
   },
   {
diff --git a/tests/test_surface/test_mesh_surface.py b/tests/test_surface/test_mesh_surface.py
new file mode 100644
index 00000000..6db08e62
--- /dev/null
+++ b/tests/test_surface/test_mesh_surface.py
@@ -0,0 +1,43 @@
+__doc__ = """Tests for mesh surface class"""
+
+import pytest
+from elastica.surface.mesh_surface import MeshSurface
+from sys import platform
+
+"""
+A dummy cube mesh stl file is used for testing at tests/cube.stl
+This dummy file was created using the open source code for creating meshes using 'numpy-stl'
+in numpy-stl documentation (https://numpy-stl.readthedocs.io/en/latest/usage.html#initial-usage)
+"""
+
+
+def test_mesh_surface_init():
+    """
+    Testing mesh_surface initialization by providing a valid path to a stl mesh file.
+    """
+    if platform == "win32":
+        path = r"tests\cube.stl"
+    else:
+        path = r"tests/cube.stl"
+
+    test_mesh_surface = MeshSurface(path)
+
+    assert isinstance(test_mesh_surface, MeshSurface), True
+
+
+def test_mesh_surface_init_error_message():
+    """
+    Testing mesh_surface error message by providing an invalid path to a stl mesh file.
+    """
+    if platform == "win32":
+        path = r"cube.stl"
+    else:
+        path = r"cube.stl"
+
+    with pytest.raises(FileNotFoundError) as excinfo:
+        test_mesh_surface = MeshSurface(path)
+
+    assert "Please check the filepath.\n"
+    " Please be sure to add .stl / .obj at the end of the filepath, if already present, ignore" in str(
+        excinfo.value
+    )

From 780dfcdc40b7708a159da0b64e289b407e073d20 Mon Sep 17 00:00:00 2001
From: Rahul-JOON <rahuljoon16@gmail.com>
Date: Thu, 24 Aug 2023 21:37:30 +0530
Subject: [PATCH 2/5] typehint error

---
 elastica/contact_forces.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/elastica/contact_forces.py b/elastica/contact_forces.py
index 3df65e83..8020dd3f 100644
--- a/elastica/contact_forces.py
+++ b/elastica/contact_forces.py
@@ -893,7 +893,7 @@ def _check_systems_validity(
 
     def apply_normal_force(
         self, system_one: RodType, system_two: AllowedContactType
-    ) -> tuple[np.array, np.array]:
+    ) -> tuple:
         """
         In the case of contact with the plane, this function computes the plane reaction force on the element.
 

From 24d919f5ead9359f220a825c2527cc8dc54f531f Mon Sep 17 00:00:00 2001
From: Rahul-JOON <rahuljoon16@gmail.com>
Date: Fri, 25 Aug 2023 23:00:58 +0530
Subject: [PATCH 3/5] Feedback1

---
 elastica/contact_forces.py              | 53 ++++++++++++------
 elastica/contact_utils.py               | 71 +++++++++++++++++--------
 elastica/mesh/mesh_initializer.py       | 17 ++++--
 elastica/surface/__init__.py            |  1 -
 elastica/surface/mesh_surface.py        | 19 +++----
 tests/test_mesh_init.py                 | 33 ++++++++++++
 tests/test_surface/test_mesh_surface.py | 20 +------
 7 files changed, 144 insertions(+), 70 deletions(-)

diff --git a/elastica/contact_forces.py b/elastica/contact_forces.py
index 8020dd3f..b0ae05d8 100644
--- a/elastica/contact_forces.py
+++ b/elastica/contact_forces.py
@@ -426,7 +426,7 @@ def _calculate_contact_forces_self_rod(
 
 
 @numba.njit(cache=True)
-def apply_normal_force_numba(
+def _calculate_contact_forces_rod_mesh_surface(
     faces_normals,
     faces_centers,
     element_position,
@@ -464,14 +464,14 @@ 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
@@ -479,7 +479,7 @@ def apply_normal_force_numba(
     # 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 +825,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 +835,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 +853,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 +899,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 +947,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,
diff --git a/elastica/contact_utils.py b/elastica/contact_utils.py
index a97f8f51..9477b3f6 100644
--- a/elastica/contact_utils.py
+++ b/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))
+"""
diff --git a/elastica/mesh/mesh_initializer.py b/elastica/mesh/mesh_initializer.py
index 2320c185..39925885 100644
--- a/elastica/mesh/mesh_initializer.py
+++ b/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:
         """
diff --git a/elastica/surface/__init__.py b/elastica/surface/__init__.py
index aee9393d..2fb544a7 100644
--- a/elastica/surface/__init__.py
+++ b/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
-
diff --git a/elastica/surface/mesh_surface.py b/elastica/surface/mesh_surface.py
index 72c6dea7..ade361da 100644
--- a/elastica/surface/mesh_surface.py
+++ b/elastica/surface/mesh_surface.py
@@ -2,10 +2,9 @@
 
 from elastica.surface.surface_base import SurfaceBase
 from elastica.mesh.mesh_initializer import Mesh
-import os
 
 
-class MeshSurface(SurfaceBase):
+class MeshSurface(SurfaceBase, Mesh):
     def __init__(self, filepath: str) -> None:
         """
         Mesh surface initializer.
@@ -43,6 +42,9 @@ def __init__(self, filepath: str) -> None:
             - Dimension: (3 spatial coordinates, 3 orthonormal basis vectors)
             - Initial value: [[1,0,0],[0,1,0],[0,0,1]]
 
+        mesh_surface.model_path:
+            - Stores the path to the mesh file.
+
         Methods:
         --------
 
@@ -72,13 +74,8 @@ def __init__(self, filepath: str) -> None:
         ex : mesh.rotate(np.array([1,0,0]), 90)
             - This method rotates the mesh by a given angle about a given axis.
         """
+        SurfaceBase.__init__(self)  # inhereting the parent class attributes, methods
+        Mesh.__init__(self, filepath)
 
-        if os.path.isfile(filepath):
-            mesh_surface = Mesh(filepath)
-            mesh_surface.mesh_update()  # using a mesh method so flake doesn't throw the error of unused variable
-
-        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"
-            )
+        mesh_surface = Mesh(filepath)
+        mesh_surface.mesh_update()  # using a mesh method so flake doesn't throw the error of unused variable
diff --git a/tests/test_mesh_init.py b/tests/test_mesh_init.py
index 038fb9d3..ba110d77 100644
--- a/tests/test_mesh_init.py
+++ b/tests/test_mesh_init.py
@@ -4,6 +4,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 from sys import platform
+import pytest
 
 
 """
@@ -13,6 +14,24 @@
 """
 
 
+def test_mesh_init_error_message():
+    """
+    Testing mesh_surface error message by providing an invalid path to a stl mesh file.
+    """
+    if platform == "win32":
+        path = r"cube.stl"
+    else:
+        path = r"cube.stl"
+
+    with pytest.raises(FileNotFoundError) as excinfo:
+        test_mesh_surface = Mesh(path)
+
+    assert "Please check the filepath.\n"
+    " Please be sure to add .stl / .obj at the end of the filepath, if already present, ignore" in str(
+        excinfo.value
+    )
+
+
 def cube_mesh_stl_init():
     """
     This function initializes a new cube mesh in stl file format.
@@ -242,3 +261,17 @@ def test_mesh_orientation():
     calculated_mesh_orientation = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
     assert_allclose(mockmeshstl.mesh_orientation, calculated_mesh_orientation)
     assert_allclose(mockmeshobj.mesh_orientation, calculated_mesh_orientation)
+
+
+def test_model_path():
+    """
+    This functions tests the model path of the cube mesh.
+    """
+    if platform == "win32":
+        path = r"tests\cube.stl"
+    else:
+        path = r"tests/cube.stl"
+
+    mockmeshstl = Mesh(path)
+
+    assert mockmeshstl.model_path == path
diff --git a/tests/test_surface/test_mesh_surface.py b/tests/test_surface/test_mesh_surface.py
index 6db08e62..ababd15a 100644
--- a/tests/test_surface/test_mesh_surface.py
+++ b/tests/test_surface/test_mesh_surface.py
@@ -1,6 +1,5 @@
 __doc__ = """Tests for mesh surface class"""
 
-import pytest
 from elastica.surface.mesh_surface import MeshSurface
 from sys import platform
 
@@ -23,21 +22,4 @@ def test_mesh_surface_init():
     test_mesh_surface = MeshSurface(path)
 
     assert isinstance(test_mesh_surface, MeshSurface), True
-
-
-def test_mesh_surface_init_error_message():
-    """
-    Testing mesh_surface error message by providing an invalid path to a stl mesh file.
-    """
-    if platform == "win32":
-        path = r"cube.stl"
-    else:
-        path = r"cube.stl"
-
-    with pytest.raises(FileNotFoundError) as excinfo:
-        test_mesh_surface = MeshSurface(path)
-
-    assert "Please check the filepath.\n"
-    " Please be sure to add .stl / .obj at the end of the filepath, if already present, ignore" in str(
-        excinfo.value
-    )
+    assert test_mesh_surface.model_path == path

From f01295152510d2c80ba969055aef159e73c840df Mon Sep 17 00:00:00 2001
From: Rahul-JOON <rahuljoon16@gmail.com>
Date: Fri, 8 Sep 2023 20:40:58 +0530
Subject: [PATCH 4/5] docs update

---
 elastica/contact_forces.py | 55 +++++++++++++++++++++++++-------------
 1 file changed, 36 insertions(+), 19 deletions(-)

diff --git a/elastica/contact_forces.py b/elastica/contact_forces.py
index b0ae05d8..7dd2168f 100644
--- a/elastica/contact_forces.py
+++ b/elastica/contact_forces.py
@@ -427,22 +427,19 @@ def _calculate_contact_forces_self_rod(
 
 @numba.njit(cache=True)
 def _calculate_contact_forces_rod_mesh_surface(
-    faces_normals,
-    faces_centers,
-    element_position,
-    position_idx_array,
+    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
@@ -450,7 +447,30 @@ def _calculate_contact_forces_rod_mesh_surface(
 
     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
     -------
@@ -965,9 +985,6 @@ def apply_contact(
             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,
         )

From 046a3f098d9ed354c69dbb1bf3514c8fda6c87a1 Mon Sep 17 00:00:00 2001
From: Rahul-JOON <rahuljoon16@gmail.com>
Date: Wed, 13 Sep 2023 20:33:24 +0530
Subject: [PATCH 5/5] Typo

---
 elastica/contact_utils.py | 46 ++++++++++++++++++++++++++++++++-------
 1 file changed, 38 insertions(+), 8 deletions(-)

diff --git a/elastica/contact_utils.py b/elastica/contact_utils.py
index 9477b3f6..c1b56db6 100644
--- a/elastica/contact_utils.py
+++ b/elastica/contact_utils.py
@@ -190,7 +190,6 @@ def _prune_using_aabbs_rod_rod(
 
 
 def find_contact_faces_idx(faces_grid, x_min, y_min, grid_size, position_collection):
-
     element_position = node_to_element_position(position_collection)
     n_element = element_position.shape[-1]
     position_idx_array = np.empty((0))
@@ -234,7 +233,7 @@ def find_contact_faces_idx(faces_grid, x_min, y_min, grid_size, position_collect
         if face_idx_element_no_duplicates.size == 0:
             raise RuntimeError(
                 "Rod object out of grid bounds"
-            )  # a snake element is on four grids with no faces
+            )  # a rod element is on four grids with no faces
 
         face_idx_array = np.concatenate(
             (face_idx_array, face_idx_element_no_duplicates)
@@ -249,9 +248,14 @@ def find_contact_faces_idx(faces_grid, x_min, y_min, grid_size, position_collect
     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):
+def surface_grid_numba(
+    faces, grid_size, face_x_left, face_x_right, face_y_down, face_y_up
+):
+    """
+    Computes the faces_grid dictionary for rod-meshsurface contact
+    Consider calling surface_grid for face_grid generation
+    """
     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))
@@ -263,16 +267,42 @@ def surface_grid_numba(faces, grid_size, face_x_left, face_x_right, face_y_down,
         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]
+            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):
+    """
+    Returns the faces_grid dictionary for rod-meshsurface contact
+    This function only creates the faces_grid dict;
+    the user must store the data in a binary file using pickle.dump
+    """
     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))
-"""
+    return dict(
+        surface_grid_numba(
+            faces, grid_size, face_x_left, face_x_right, face_y_down, face_y_up
+        )
+    )