Merge pull request GazzolaLab#328 from Ali-7800/113_dev_contact

Cylinder-plane contact class and tests

.coveragerc

@@ -8,6 +8,8 @@ exclude_lines =
     def __repr__
     from
     import
+    raise AssertionError
+    raise NotImplementedError
 show_missing = true
 
 [run]

elastica/_contact_functions.py

@@ -780,3 +780,69 @@ def _calculate_contact_forces_rod_plane_with_anisotropic_friction(
         director_collection,
         _batch_cross(torque_arm, static_friction_force_along_rolling_direction),
     )
+
+
+@numba.njit(cache=True)
+def _calculate_contact_forces_cylinder_plane(
+    plane_origin,
+    plane_normal,
+    surface_tol,
+    k,
+    nu,
+    length,
+    position_collection,
+    velocity_collection,
+    external_forces,
+):
+
+    # Compute plane response force
+    # total_forces = system.internal_forces + system.external_forces
+    total_forces = external_forces
+    force_component_along_normal_direction = _batch_product_i_ik_to_k(
+        plane_normal, total_forces
+    )
+    forces_along_normal_direction = _batch_product_i_k_to_ik(
+        plane_normal, force_component_along_normal_direction
+    )
+    # If the total force component along the plane normal direction is greater than zero that means,
+    # total force is pushing rod away from the plane not towards the plane. Thus, response force
+    # applied by the surface has to be zero.
+    forces_along_normal_direction[
+        ..., np.where(force_component_along_normal_direction > 0)[0]
+    ] = 0.0
+    # Compute response force on the element. Plane response force
+    # has to be away from the surface and towards the element. Thus
+    # multiply forces along normal direction with negative sign.
+    plane_response_force = -forces_along_normal_direction
+
+    # Elastic force response due to penetration
+    element_position = position_collection
+    distance_from_plane = _batch_product_i_ik_to_k(
+        plane_normal, (element_position - plane_origin)
+    )
+    plane_penetration = np.minimum(distance_from_plane - length / 2, 0.0)
+    elastic_force = -k * _batch_product_i_k_to_ik(plane_normal, plane_penetration)
+
+    # Damping force response due to velocity towards the plane
+    element_velocity = velocity_collection
+    normal_component_of_element_velocity = _batch_product_i_ik_to_k(
+        plane_normal, element_velocity
+    )
+    damping_force = -nu * _batch_product_i_k_to_ik(
+        plane_normal, normal_component_of_element_velocity
+    )
+
+    # Compute total plane response force
+    plane_response_force_total = plane_response_force + elastic_force + damping_force
+
+    # Check if the rigid body is in contact with plane.
+    no_contact_point_idx = np.where((distance_from_plane - length / 2) > 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[..., no_contact_point_idx] = 0.0
+    plane_response_force_total[..., no_contact_point_idx] = 0.0
+
+    # Update the external forces
+    external_forces += plane_response_force_total
+
+    return (_batch_norm(plane_response_force), no_contact_point_idx)

elastica/contact_forces.py

@@ -16,6 +16,7 @@
     _calculate_contact_forces_rod_sphere,
     _calculate_contact_forces_rod_plane,
     _calculate_contact_forces_rod_plane_with_anisotropic_friction,
+    _calculate_contact_forces_cylinder_plane,
 )
 import numpy as np
 
@@ -707,3 +708,89 @@ def apply_contact(self, system_one: RodType, system_two: SystemType) -> None:
             system_one.internal_torques,
             system_one.external_torques,
         )
+
+
+class CylinderPlaneContact(NoContact):
+    """
+    This class is for applying contact forces between cylinder-plane.
+    First system is always cylinder and second system is always plane.
+    For more details regarding the contact module refer to
+    Eqn 4.8 of Gazzola et al. RSoS (2018).
+
+        How to define contact between cylinder and plane.
+    >>> simulator.detect_contact_between(cylinder, plane).using(
+    ...    CylinderPlaneContact,
+    ...    k=1e4,
+    ...    nu=10,
+    ... )
+    """
+
+    def __init__(
+        self,
+        k: float,
+        nu: float,
+    ):
+        """
+        Parameters
+        ----------
+        k : float
+            Contact spring constant.
+        nu : float
+            Contact damping constant.
+        """
+        super(CylinderPlaneContact, self).__init__()
+        self.k = k
+        self.nu = nu
+        self.surface_tol = 1e-4
+
+    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 RodPlaneContact class first_system should be a cylinder and second_system should be a plane.
+        Parameters
+        ----------
+        system_one
+            SystemType
+        system_two
+            AllowedContactType
+        """
+        if not issubclass(system_one.__class__, Cylinder) or not issubclass(
+            system_two.__class__, Plane
+        ):
+            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 cylinder, second should be a plane".format(
+                    system_one.__class__, system_two.__class__
+                )
+            )
+
+    def apply_contact(self, system_one: Cylinder, system_two: SystemType):
+        """
+        This function computes the plane force response on the rigid body, in the
+        case of contact. Contact model given in Eqn 4.8 Gazzola et. al. RSoS 2018 paper
+        is used.
+
+        Parameters
+        ----------
+        system_one: object
+            cylinder object.
+        system_two: object
+            Plane object.
+
+        """
+        return _calculate_contact_forces_cylinder_plane(
+            system_two.origin,
+            system_two.normal,
+            self.surface_tol,
+            self.k,
+            self.nu,
+            system_one.length,
+            system_one.position_collection,
+            system_one.velocity_collection,
+            system_one.external_forces,
+        )

elastica/interaction.py

@@ -4,18 +4,14 @@
 import numpy as np
 from elastica.external_forces import NoForces
 from numba import njit
-from elastica._linalg import (
-    _batch_norm,
-    _batch_product_i_k_to_ik,
-    _batch_product_i_ik_to_k,
-)
 from elastica.contact_utils import (
     _elements_to_nodes_inplace,
     _node_to_element_velocity,
 )
 from elastica._contact_functions import (
     _calculate_contact_forces_rod_plane,
     _calculate_contact_forces_rod_plane_with_anisotropic_friction,
+    _calculate_contact_forces_cylinder_plane,
 )
 
 
@@ -548,7 +544,7 @@ def apply_normal_force(self, system):
         -------
         magnitude of the plane response
         """
-        return apply_normal_force_numba_rigid_body(
+        return _calculate_contact_forces_cylinder_plane(
             self.plane_origin,
             self.plane_normal,
             self.surface_tol,
@@ -574,54 +570,8 @@ def apply_normal_force_numba_rigid_body(
     external_forces,
 ):
 
-    # Compute plane response force
-    # total_forces = system.internal_forces + system.external_forces
-    total_forces = external_forces
-    force_component_along_normal_direction = _batch_product_i_ik_to_k(
-        plane_normal, total_forces
-    )
-    forces_along_normal_direction = _batch_product_i_k_to_ik(
-        plane_normal, force_component_along_normal_direction
-    )
-    # If the total force component along the plane normal direction is greater than zero that means,
-    # total force is pushing rod away from the plane not towards the plane. Thus, response force
-    # applied by the surface has to be zero.
-    forces_along_normal_direction[
-        ..., np.where(force_component_along_normal_direction > 0)[0]
-    ] = 0.0
-    # Compute response force on the element. Plane response force
-    # has to be away from the surface and towards the element. Thus
-    # multiply forces along normal direction with negative sign.
-    plane_response_force = -forces_along_normal_direction
-
-    # Elastic force response due to penetration
-    element_position = position_collection
-    distance_from_plane = _batch_product_i_ik_to_k(
-        plane_normal, (element_position - plane_origin)
-    )
-    plane_penetration = np.minimum(distance_from_plane - length / 2, 0.0)
-    elastic_force = -k * _batch_product_i_k_to_ik(plane_normal, plane_penetration)
-
-    # Damping force response due to velocity towards the plane
-    element_velocity = velocity_collection
-    normal_component_of_element_velocity = _batch_product_i_ik_to_k(
-        plane_normal, element_velocity
-    )
-    damping_force = -nu * _batch_product_i_k_to_ik(
-        plane_normal, normal_component_of_element_velocity
+    raise NotImplementedError(
+        "This function is removed in v0.3.2. For cylinder plane contact please use: \n"
+        "elastica._contact_functions._calculate_contact_forces_cylinder_plane() \n"
+        "For detail, refer to issue #113."
     )
-
-    # Compute total plane response force
-    plane_response_force_total = plane_response_force + elastic_force + damping_force
-
-    # Check if the rigid body is in contact with plane.
-    no_contact_point_idx = np.where((distance_from_plane - length / 2) > 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[..., no_contact_point_idx] = 0.0
-    plane_response_force_total[..., no_contact_point_idx] = 0.0
-
-    # Update the external forces
-    external_forces += plane_response_force_total
-
-    return (_batch_norm(plane_response_force), no_contact_point_idx)