Merge pull request #379 from GazzolaLab/feature/order-operation

Explicit operation order determined by user definition

docs/guide/workflow.md

@@ -89,7 +89,7 @@ This can be repeated to create multiple rods. Supported geometries are listed in
 The number of element (`n_elements`) and `base_length` determines the spatial discretization `dx`. More detail discussion is included [here](discretization.md).
 :::
 
-<h2>3. Define Boundary Conditions, Forcings, Damping and Connections</h2>
+<h2>3.a Define Boundary Conditions, Forcings, and Connections</h2>
 
 Now that we have added all our rods to `SystemSimulator`, we
 need to apply relevant boundary conditions.
@@ -123,6 +123,34 @@ SystemSimulator.add_forcing_to(rod1).using(
 )
 ```
 
+One last condition we can define is the connections between rods. See [this page](../api/connections.rst) for in-depth explanations and documentation.
+
+```python
+from elastica.connections import FixedJoint
+
+# Connect rod 1 and rod 2. '_connect_idx' specifies the node number that
+# the connection should be applied to. You are specifying the index of a
+# list so you can use -1 to access the last node.
+SystemSimulator.connect(
+    first_rod  = rod1,
+    second_rod = rod2,
+    first_connect_idx  = -1, # Connect to the last node of the first rod.
+    second_connect_idx =  0  # Connect to first node of the second rod.
+    ).using(
+        FixedJoint,  # Type of connection between rods
+        k  = 1e5,    # Spring constant of force holding rods together (F = k*x)
+        nu = 0,      # Energy dissipation of joint
+        kt = 5e3     # Rotational stiffness of rod to avoid rods twisting
+        )
+```
+
+:::{note}
+Version 0.3.3: The order of the operation is defined by the order of the definition. For example, if you define the connection before the forcing condition, the connection will be applied first. This is less important for the boundary condition, forcing, and connection since they do not depend on each other. However, it is important for friction, contact, or any custom boundary conditions since they depend on other boundary conditions.
+For example, friction should be defined after contact, since contact will define the normal force applied to the surface, which friction depends on. Contact should be defined before any other boundary conditions, since aggregated normal force is used to calculate the repelling force.
+:::
+
+<h2>3.b Define Damping </h2>
+
 Next, if required, in order to numerically stabilize the simulation,
 we can apply damping to the rods.
 See [this page](../api/damping.rst) for in-depth explanations and documentation.
@@ -146,26 +174,6 @@ SystemSimulator.dampin(rod2).using(
 )
 ```
 
-One last condition we can define is the connections between rods. See [this page](../api/connections.rst) for in-depth explanations and documentation.
-
-```python
-from elastica.connections import FixedJoint
-
-# Connect rod 1 and rod 2. '_connect_idx' specifies the node number that
-# the connection should be applied to. You are specifying the index of a
-# list so you can use -1 to access the last node.
-SystemSimulator.connect(
-    first_rod  = rod1,
-    second_rod = rod2,
-    first_connect_idx  = -1, # Connect to the last node of the first rod.
-    second_connect_idx =  0  # Connect to first node of the second rod.
-    ).using(
-        FixedJoint,  # Type of connection between rods
-        k  = 1e5,    # Spring constant of force holding rods together (F = k*x)
-        nu = 0,      # Energy dissipation of joint
-        kt = 5e3     # Rotational stiffness of rod to avoid rods twisting
-        )
-```
 
 <h2>4. Add Callback Functions (optional)</h2>
 

elastica/modules/base_system.py

@@ -5,16 +5,19 @@
 Basic coordinating for multiple, smaller systems that have an independently integrable
 interface (i.e. works with symplectic or explicit routines `timestepper.py`.)
 """
-from typing import Iterable, Callable, AnyStr
+from typing import AnyStr, Iterable
+from elastica.typing import OperatorType, OperatorCallbackType, OperatorFinalizeType
 
 from collections.abc import MutableSequence
 
 from elastica.rod import RodBase
 from elastica.rigidbody import RigidBodyBase
 from elastica.surface import SurfaceBase
-from elastica.modules.memory_block import construct_memory_block_structures
 from elastica._synchronize_periodic_boundary import _ConstrainPeriodicBoundaries
 
+from .memory_block import construct_memory_block_structures
+from .operator_group import OperatorGroupFIFO
+
 
 class BaseSystemCollection(MutableSequence):
     """
@@ -29,11 +32,9 @@ class BaseSystemCollection(MutableSequence):
         _systems: list
             List of rod-like objects.
 
-    """
-
-    """
     Developer Note
     -----
+
     Note
     ----
     We can directly subclass a list for the
@@ -45,13 +46,11 @@ def __init__(self):
         # Collection of functions. Each group is executed as a collection at the different steps.
         # Each component (Forcing, Connection, etc.) registers the executable (callable) function
         # in the group that that needs to be executed. These should be initialized before mixin.
-        self._feature_group_synchronize: Iterable[Callable[[float], None]] = []
-        self._feature_group_constrain_values: Iterable[Callable[[float], None]] = []
-        self._feature_group_constrain_rates: Iterable[Callable[[float], None]] = []
-        self._feature_group_callback: Iterable[Callable[[float, int, AnyStr], None]] = (
-            []
-        )
-        self._feature_group_finalize: Iterable[Callable] = []
+        self._feature_group_synchronize: Iterable[OperatorType] = OperatorGroupFIFO()
+        self._feature_group_constrain_values: Iterable[OperatorType] = []
+        self._feature_group_constrain_rates: Iterable[OperatorType] = []
+        self._feature_group_callback: Iterable[OperatorCallbackType] = []
+        self._feature_group_finalize: Iterable[OperatorFinalizeType] = []
         # We need to initialize our mixin classes
         super(BaseSystemCollection, self).__init__()
         # List of system types/bases that are allowed
@@ -169,34 +168,23 @@ def finalize(self):
 
         # Toggle the finalize_flag
         self._finalize_flag = True
-        # sort _feature_group_synchronize so that _call_contacts is at the end
-        _call_contacts_index = []
-        for idx, feature in enumerate(self._feature_group_synchronize):
-            if feature.__name__ == "_call_contacts":
-                _call_contacts_index.append(idx)
-
-        # Move to the _call_contacts to the end of the _feature_group_synchronize list.
-        for index in _call_contacts_index:
-            self._feature_group_synchronize.append(
-                self._feature_group_synchronize.pop(index)
-            )
 
     def synchronize(self, time: float):
         # Collection call _feature_group_synchronize
-        for feature in self._feature_group_synchronize:
-            feature(time)
+        for func in self._feature_group_synchronize:
+            func(time=time)
 
     def constrain_values(self, time: float):
         # Collection call _feature_group_constrain_values
-        for feature in self._feature_group_constrain_values:
-            feature(time)
+        for func in self._feature_group_constrain_values:
+            func(time=time)
 
     def constrain_rates(self, time: float):
         # Collection call _feature_group_constrain_rates
-        for feature in self._feature_group_constrain_rates:
-            feature(time)
+        for func in self._feature_group_constrain_rates:
+            func(time=time)
 
     def apply_callbacks(self, time: float, current_step: int):
         # Collection call _feature_group_callback
-        for feature in self._feature_group_callback:
-            feature(time, current_step)
+        for func in self._feature_group_callback:
+            func(time=time, current_step=current_step)

elastica/modules/connections.py

@@ -24,7 +24,6 @@ class Connections:
     def __init__(self):
         self._connections = []
         super(Connections, self).__init__()
-        self._feature_group_synchronize.append(self._call_connections)
         self._feature_group_finalize.append(self._finalize_connections)
 
     def connect(
@@ -60,48 +59,57 @@ def connect(
         sys_dofs = [self._systems[idx].n_elems for idx in sys_idx]
 
         # Create _Connect object, cache it and return to user
-        _connector = _Connect(*sys_idx, *sys_dofs)
-        _connector.set_index(first_connect_idx, second_connect_idx)
-        self._connections.append(_connector)
+        _connect = _Connect(*sys_idx, *sys_dofs)
+        _connect.set_index(first_connect_idx, second_connect_idx)
+        self._connections.append(_connect)
+        self._feature_group_synchronize.append_id(_connect)
 
-        return _connector
+        return _connect
 
     def _finalize_connections(self):
         # From stored _Connect objects, instantiate the joints and store it
-
         # dev : the first indices stores the
         # (first rod index, second_rod_idx, connection_idx_on_first_rod, connection_idx_on_second_rod)
-        # to apply the connections to
-        # Technically we can use another array but it its one more book-keeping
-        # step. Being lazy, I put them both in the same array
-        self._connections[:] = [
-            (*connection.id(), connection()) for connection in self._connections
-        ]
+        # to apply the connections to.
+
+        for connection in self._connections:
+            first_sys_idx, second_sys_idx, first_connect_idx, second_connect_idx = (
+                connection.id()
+            )
+            connect_instance = connection.instantiate()
+
+            # FIXME: lambda t is included because OperatorType takes time as an argument
+            def apply_forces(time):
+                connect_instance.apply_forces(
+                    system_one=self._systems[first_sys_idx],
+                    index_one=first_connect_idx,
+                    system_two=self._systems[second_sys_idx],
+                    index_two=second_connect_idx,
+                )
+
+            def apply_torques(time):
+                connect_instance.apply_torques(
+                    system_one=self._systems[first_sys_idx],
+                    index_one=first_connect_idx,
+                    system_two=self._systems[second_sys_idx],
+                    index_two=second_connect_idx,
+                )
+
+            self._feature_group_synchronize.add_operators(
+                connection, [apply_forces, apply_torques]
+            )
+
+            self.warnings(connection)
+
+        self._connections = []
+        del self._connections
 
         # Need to finally solve CPP here, if we are doing things properly
         # This is to optimize the call tree for better memory accesses
         # https://brooksandrew.github.io/simpleblog/articles/intro-to-graph-optimization-solving-cpp/
 
-    def _call_connections(self, *args, **kwargs):
-        for (
-            first_sys_idx,
-            second_sys_idx,
-            first_connect_idx,
-            second_connect_idx,
-            connection,
-        ) in self._connections:
-            connection.apply_forces(
-                self._systems[first_sys_idx],
-                first_connect_idx,
-                self._systems[second_sys_idx],
-                second_connect_idx,
-            )
-            connection.apply_torques(
-                self._systems[first_sys_idx],
-                first_connect_idx,
-                self._systems[second_sys_idx],
-                second_connect_idx,
-            )
+    def warnings(self, connection):
+        pass
 
 
 class _Connect:
@@ -152,7 +160,7 @@ def __init__(
     def set_index(self, first_idx, second_idx):
         # TODO assert range
         # First check if the types of first rod idx and second rod idx variable are same.
-        assert type(first_idx) == type(
+        assert type(first_idx) is type(
             second_idx
         ), "Type of first_connect_idx :{}".format(
             type(first_idx)
@@ -265,7 +273,7 @@ def id(self):
             self.second_sys_connection_idx,
         )
 
-    def __call__(self, *args, **kwargs):
+    def instantiate(self):
         if not self._connect_cls:
             raise RuntimeError(
                 "No connections provided to link rod id {0}"

elastica/modules/contact.py

@@ -5,9 +5,11 @@
 Provides the contact interface to apply contact forces between objects
 (rods, rigid bodies, surfaces).
 """
-
+import logging
 from elastica.typing import SystemType, AllowedContactType
 
+logger = logging.getLogger(__name__)
+
 
 class Contact:
     """
@@ -23,7 +25,6 @@ class Contact:
     def __init__(self):
         self._contacts = []
         super(Contact, self).__init__()
-        self._feature_group_synchronize.append(self._call_contacts)
         self._feature_group_finalize.append(self._finalize_contact)
 
     def detect_contact_between(
@@ -51,6 +52,7 @@ def detect_contact_between(
         # Create _Contact object, cache it and return to user
         _contact = _Contact(*sys_idx)
         self._contacts.append(_contact)
+        self._feature_group_synchronize.append_id(_contact)
 
         return _contact
 
@@ -61,30 +63,36 @@ def _finalize_contact(self) -> None:
         # to apply the contacts to
         # Technically we can use another array but it its one more book-keeping
         # step. Being lazy, I put them both in the same array
-        self._contacts[:] = [(*contact.id(), contact()) for contact in self._contacts]
-
-        # check contact order
-        for (
-            first_sys_idx,
-            second_sys_idx,
-            contact,
-        ) in self._contacts:
-            contact._check_systems_validity(
-                self._systems[first_sys_idx],
-                self._systems[second_sys_idx],
-            )
+        for contact in self._contacts:
+            first_sys_idx, second_sys_idx = contact.id()
+            contact_instance = contact.instantiate()
 
-    def _call_contacts(self, time: float):
-        for (
-            first_sys_idx,
-            second_sys_idx,
-            contact,
-        ) in self._contacts:
-            contact.apply_contact(
+            contact_instance._check_systems_validity(
                 self._systems[first_sys_idx],
                 self._systems[second_sys_idx],
             )
 
+            def apply_contact(time):
+                contact_instance.apply_contact(
+                    system_one=self._systems[first_sys_idx],
+                    system_two=self._systems[second_sys_idx],
+                )
+
+            self._feature_group_synchronize.add_operators(contact, [apply_contact])
+
+            self.warnings(contact)
+
+        self._contacts = []
+        del self._contacts
+
+    def warnings(self, contact):
+        from elastica.contact_forces import NoContact
+
+        # Classes that should be used last
+        if not self._feature_group_synchronize.is_last(contact):
+            if isinstance(contact._contact_cls, NoContact):
+                logger.warning("Contact features should be instantiated lastly.")
+
 
 class _Contact:
     """
@@ -153,7 +161,7 @@ def id(self):
             self.second_sys_idx,
         )
 
-    def __call__(self, *args, **kwargs):
+    def instantiate(self, *args, **kwargs):
         if not self._contact_cls:
             raise RuntimeError(
                 "No contacts provided to to establish contact between rod-like object id {0}"