Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Address #148 by adding ConfigurableFixedConstraint boundary condition class #143

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
48 commits
Select commit Hold shift + click to select a range
1836ba4
Start implementing ConfigurableFixedConstraint
mstoelzle Jul 15, 2022
aaca3f9
Replace deprecated np.bool with python bool
mstoelzle Jul 17, 2022
2ccef47
First working, preliminary implementation of `ConfigurableFixedConstr…
mstoelzle Jul 18, 2022
1c8551a
Add BC `configurable_fixed_constraint` example
mstoelzle Jul 18, 2022
c6c1018
Run make black
mstoelzle Jul 18, 2022
ba494fd
Add more code comments to `nb_constraint_rotational_values`
mstoelzle Jul 18, 2022
5b6a1b4
Fix bug in docstring example
mstoelzle Jul 18, 2022
448fcc1
Fix bug where allowed directors are not written to director collection
mstoelzle Jul 19, 2022
2e7bf34
Merge remote-tracking branch 'upstream/update-0.3.0' into feature/con…
mstoelzle Jul 19, 2022
14aaffe
Some reformatting using poetry black
mstoelzle Jul 19, 2022
bffd3cc
Rotate angular velocities to inertial frame before applying constraint
mstoelzle Jul 19, 2022
5bfdfff
Some minor suggestions to improve docstrings and code comments from c…
mstoelzle Jul 21, 2022
bd15e96
Merge branch 'update-0.3.0' into feature/configurable-custom-constraint
mstoelzle Jul 21, 2022
13c1995
Fix some merge bugs
mstoelzle Jul 21, 2022
4a9c71e
rename `ConfigurableFixedConstraint` to `ConfigurableConstraint`
mstoelzle Jul 21, 2022
2d13584
Some renaming in `nb_constrain_rotational_rates`
mstoelzle Jul 21, 2022
5d5e4b9
Merge instructions for new position values
mstoelzle Jul 21, 2022
970cf02
Rename `omega_collection_allowed` to `omega_collection_not_constrained`
mstoelzle Jul 21, 2022
0cc1705
Deactivate `nb_constraint_rotational_values` in `ConfigurableConstraint`
mstoelzle Jul 21, 2022
5451286
Change filename of position plot
mstoelzle Jul 21, 2022
0630bb8
Deactivate testing of `nb_constraint_rotational_values`
mstoelzle Jul 21, 2022
ab57775
Rename `configurable_fixed_constraint.py` to `configurable_constraint…
mstoelzle Jul 21, 2022
effff14
Rename `ConfigurableFixedConstraintSimulator` to `ConfigurableConstra…
mstoelzle Jul 21, 2022
4f0bcda
Add `test_configurable_constraint`
mstoelzle Jul 21, 2022
5a888fb
Fully remove `nb_constraint_rotational_values` from `ConfigurableCons…
mstoelzle Jul 21, 2022
1e5e0dc
Apply suggestions to improve docstrings
mstoelzle Jul 22, 2022
9486367
Merge branch 'update-0.3.0' into feature/configurable-custom-constraint
mstoelzle Jul 22, 2022
bb068fb
Rename constraint variables to appropiate names in `test_boundary_con…
mstoelzle Jul 22, 2022
33d42f1
Improve `test_configurable_constraint`
mstoelzle Jul 22, 2022
0eddbbf
Apply formatting
mstoelzle Jul 22, 2022
4e8c6e2
Add special case of 90 degree yaw to `test_configurable_constraint`
mstoelzle Jul 22, 2022
8d91580
Adapt setting of `ConfigurableConstraint` example
mstoelzle Jul 22, 2022
4f2e7a2
Some formatting
mstoelzle Jul 22, 2022
86d016c
Set simulation duration to 10s
mstoelzle Jul 22, 2022
ef7954f
Revert lock file to 0.3.0 branch
mstoelzle Jul 22, 2022
4d53264
Remove _inv_rotate import from `boundary_conditions.py`
mstoelzle Jul 23, 2022
44bb1b4
Use `batch_matrix_transpose` instead of `np.transpose`
mstoelzle Jul 23, 2022
51f3a91
Fix syntax errors of previous commits
mstoelzle Jul 23, 2022
4c06b0e
Rename `ConfigurableConstraint` to `GeneralConstraint`
mstoelzle Jul 23, 2022
7020480
Merge branch 'update-0.3.0' into feature/configurable-custom-constraint
mstoelzle Jul 23, 2022
b8eca47
Restore logic of legacy `FixedConstraint` implementation
mstoelzle Jul 23, 2022
bdec395
Add README entry for `BoundaryConditionsCases`
mstoelzle Jul 23, 2022
eb98cd7
Merge branch 'update-0.3.0' into feature/configurable-custom-constraint
mstoelzle Jul 24, 2022
a68f2e2
Implement some docstring and typing suggestions
mstoelzle Jul 24, 2022
5325532
Fix some bugs in previous commit
mstoelzle Jul 24, 2022
03fc150
Update initialisation of `translational_constraint_selector` and `rot…
mstoelzle Jul 24, 2022
2cdc2f9
Add `GeneralConstraint` to `docs/api/constraints.rst`
mstoelzle Jul 24, 2022
4682601
Add `GeneralConstraint` also to compatibility and description section
mstoelzle Jul 24, 2022
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions docs/api/constraints.rst
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ Constraints are equivalent to displacement boundary condition.
ConstraintBase
FreeBC
OneEndFixedBC
GeneralConstraint
FixedConstraint
HelicalBucklingBC
FreeRod
Expand All @@ -31,6 +32,7 @@ Constraint / Boundary Condition Rod Rigid Body
=============================== ==== ===========
FreeBC ✅ ✅
OneEndFixedBC ✅ ✅
GeneralConstraint ✅ ✅
FixedConstraint ✅ ✅
HelicalBucklingBC ✅ ❌
=============================== ==== ===========
Expand Down Expand Up @@ -62,6 +64,9 @@ Built-in Constraints
.. autoclass:: OneEndFixedBC
:special-members: __init__

.. autoclass:: GeneralConstraint
:special-members: __init__

.. autoclass:: FixedConstraint
:special-members: __init__

Expand Down
242 changes: 226 additions & 16 deletions elastica/boundary_conditions.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
"FreeRod", # Deprecated: remove v0.3.0
"OneEndFixedBC",
"OneEndFixedRod", # Deprecated: remove v0.3.0
"GeneralConstraint",
"FixedConstraint",
"HelicalBucklingBC",
]
Expand All @@ -19,9 +20,11 @@
import numba
from numba import njit

from elastica._linalg import _batch_matvec, _batch_matrix_transpose
from elastica._rotations import _get_rotation_matrix
from elastica.rod import RodBase
from elastica.rigidbody import RigidBodyBase
from elastica.typing import SystemType


class ConstraintBase(ABC):
Expand Down Expand Up @@ -258,32 +261,41 @@ class OneEndFixedRod(OneEndFixedBC):
)


class FixedConstraint(ConstraintBase):
class GeneralConstraint(ConstraintBase):
"""
This boundary condition class fixes the specified node or orientations.
This boundary condition class allows the specified node/link to have a configurable constraint.
Index can be passed to fix either or both the position or the director.
Constraining position is equivalent to setting 0 translational DOF.
Constraining director is equivalent to setting 0 rotational DOF.

Examples
--------
How to fix two ends of the rod:
How to fix all translational and rotational dof except allowing twisting around the z-axis in an inertial frame:

>>> simulator.constrain(rod).using(
... FixedConstraint,
... constrained_position_idx=(0,1,-2,-1),
... constrained_director_idx=(0,-1)
... GeneralConstraint,
... constrained_position_idx=(0,),
... constrained_director_idx=(0,),
... translational_constraint_selector=np.array([True, True, True]),
... rotational_constraint_selector=np.array([True, True, False]),
... )

How to pin the middle of the rod (10th node), without constraining the rotational DOF.
How to allow the end of the rod to move in the XY plane and allow all rotational dof:

>>> simulator.constrain(rod).using(
... FixedConstraint,
... constrained_position_idx=(10)
... GeneralConstraint,
... constrained_position_idx=(-1,),
... translational_constraint_selector=np.array([True, True, False]),
... )
"""

def __init__(self, *fixed_data, **kwargs):
def __init__(
self,
*fixed_data,
translational_constraint_selector: Optional[np.ndarray] = None,
rotational_constraint_selector: Optional[np.array] = None,
**kwargs,
):
"""

Initialization of the constraint. Any parameter passed to 'using' will be available in kwargs.
Expand All @@ -294,6 +306,12 @@ def __init__(self, *fixed_data, **kwargs):
Tuple of position-indices that will be constrained
constrained_director_idx : tuple
Tuple of director-indices that will be constrained
translational_constraint_selector: Optional[np.ndarray]
np.array of type bool indicating which translational degrees of freedom (dof) to constrain.
If entry is True, the corresponding dof will be constrained. If None, we constrain all dofs.
rotational_constraint_selector: Optional[np.ndarray]
np.array of type bool indicating which translational degrees of freedom (dof) to constrain.
If entry is True, the corresponding dof will be constrained.
"""
super().__init__(**kwargs)
pos, dir = [], []
Expand All @@ -317,6 +335,204 @@ def __init__(self, *fixed_data, **kwargs):
# transpose from (blocksize, dim, dim) to (dim, dim, blocksize)
self.fixed_directors = np.array(dir).transpose((1, 2, 0))

if translational_constraint_selector is None:
translational_constraint_selector = np.array([True, True, True])
if rotational_constraint_selector is None:
rotational_constraint_selector = np.array([True, True, True])
# properly validate the user-provided constraint selectors
assert (
type(translational_constraint_selector) == np.ndarray
and translational_constraint_selector.dtype == bool
and translational_constraint_selector.shape == (3,)
), "Translational constraint selector must be a 1D boolean array of length 3."
assert (
type(rotational_constraint_selector) == np.ndarray
and rotational_constraint_selector.dtype == bool
and rotational_constraint_selector.shape == (3,)
), "Rotational constraint selector must be a 1D boolean array of length 3."
# cast booleans to int
self.translational_constraint_selector = (
translational_constraint_selector.astype(int)
)
self.rotational_constraint_selector = rotational_constraint_selector.astype(int)

def constrain_values(self, system: SystemType, time: float) -> None:
if self.constrained_position_idx.size:
self.nb_constrain_translational_values(
system.position_collection,
self.fixed_positions,
self.constrained_position_idx,
self.translational_constraint_selector,
)

def constrain_rates(self, system: SystemType, time: float) -> None:
if self.constrained_position_idx.size:
self.nb_constrain_translational_rates(
system.velocity_collection,
self.constrained_position_idx,
self.translational_constraint_selector,
)
if self.constrained_director_idx.size:
self.nb_constrain_rotational_rates(
system.director_collection,
system.omega_collection,
self.constrained_director_idx,
self.rotational_constraint_selector,
)

@staticmethod
@njit(cache=True)
def nb_constrain_translational_values(
position_collection, fixed_position_collection, indices, constraint_selector
) -> None:
"""
Computes constrain values in numba njit decorator

Parameters
----------
position_collection : numpy.ndarray
2D (dim, blocksize) array containing data with `float` type.
fixed_position_collection : numpy.ndarray
2D (dim, blocksize) array containing data with `float` type.
indices : numpy.ndarray
1D array containing the index of constraining nodes
constraint_selector: numpy.ndarray
1D array of type int and size (3,) indicating which translational Degrees of Freedom (DoF) to constrain.
Entries are integers in {0, 1} (e.g. a binary values of either 0 or 1).
If entry is 1, the concerning DoF will be constrained, otherwise it will be free for translation.
Selector shall be specified in the inertial frame
"""
block_size = indices.size
for i in range(block_size):
k = indices[i]
# First term: add the old position values using the inverse constraint selector (e.g. DoF)
# Second term: add the fixed position values using the constraint selector (e.g. constraint dimensions)
position_collection[..., k] = (
1 - constraint_selector
) * position_collection[
..., k
] + constraint_selector * fixed_position_collection[
..., i
]

@staticmethod
@njit(cache=True)
def nb_constrain_translational_rates(
velocity_collection, indices, constraint_selector
) -> None:
"""
Compute constrain rates in numba njit decorator

Parameters
----------
velocity_collection : numpy.ndarray
2D (dim, blocksize) array containing data with `float` type.
indices : numpy.ndarray
1D array containing the index of constraining nodes
constraint_selector: numpy.ndarray
1D array of type int and size (3,) indicating which translational Degrees of Freedom (DoF) to constrain.
Entries are integers in {0, 1} (e.g. a binary values of either 0 or 1).
If entry is 1, the concerning DoF will be constrained, otherwise it will be free for translation.
Selector shall be specified in the inertial frame
"""

block_size = indices.size
for i in range(block_size):
k = indices[i]
# set the dofs to 0 where the constraint_selector mask is active
velocity_collection[..., k] = (
1 - constraint_selector
) * velocity_collection[..., k]
mstoelzle marked this conversation as resolved.
Show resolved Hide resolved

@staticmethod
@njit(cache=True)
def nb_constrain_rotational_rates(
director_collection, omega_collection, indices, constraint_selector
) -> None:
"""
Compute constrain rates in numba njit decorator

Parameters
----------
director_collection : numpy.ndarray
2D (dim, blocksize) array containing data with `float` type.
omega_collection : numpy.ndarray
2D (dim, blocksize) array containing data with `float` type.
indices : numpy.ndarray
1D array containing the index of constraining nodes
constraint_selector: numpy.ndarray
1D array of type int and size (3,) indicating which rotational Degrees of Freedom (DoF) to constrain.
Entries are integers in {0, 1} (e.g. a binary values of either 0 or 1).
If an entry is 1, the rotation around the respective axis will be constrained,
otherwise the system can freely rotate around the axis.
The selector shall be specified in the lab frame
"""
directors = director_collection[..., indices]

# rotate angular velocities to lab frame
omega_collection_lab_frame = _batch_matvec(
_batch_matrix_transpose(directors), omega_collection[..., indices]
)

# apply constraint selector to angular velocities in lab frame
omega_collection_not_constrained = (
1 - np.expand_dims(constraint_selector, 1)
) * omega_collection_lab_frame

# rotate angular velocities vector back to local frame and apply to omega_collection
omega_collection[..., indices] = _batch_matvec(
directors, omega_collection_not_constrained
)


class FixedConstraint(GeneralConstraint):
"""
This boundary condition class fixes the specified node or orientations.
Index can be passed to fix either or both the position or the director.
Constraining position is equivalent to setting 0 translational DOF.
Constraining director is equivalent to setting 0 rotational DOF.

Examples
--------
How to fix two ends of the rod:

>>> simulator.constrain(rod).using(
... FixedConstraint,
... constrained_position_idx=(0,-1),
... constrained_director_idx=(0,-1)
... )

How to pin the middle of the rod (10th node), without constraining the rotational DOF.

>>> simulator.constrain(rod).using(
... FixedConstraint,
... constrained_position_idx=(10)
... )

See Also
---------
GeneralConstraint: Generalized constraint with configurable DOF.
"""
mstoelzle marked this conversation as resolved.
Show resolved Hide resolved

def __init__(self, *args, **kwargs):
"""

Initialization of the constraint. Any parameter passed to 'using' will be available in kwargs.

Parameters
----------
constrained_position_idx : tuple
Tuple of position-indices that will be constrained
constrained_director_idx : tuple
Tuple of director-indices that will be constrained
"""
super().__init__(
*args,
translational_constraint_selector=np.array([True, True, True]),
rotational_constraint_selector=np.array([True, True, True]),
**kwargs,
)

def constrain_values(
self, rod: Union[Type[RodBase], Type[RigidBodyBase]], time: float
) -> None:
Expand Down Expand Up @@ -354,7 +570,6 @@ def nb_constraint_rotational_values(
) -> None:
"""
Computes constrain values in numba njit decorator

Parameters
----------
director_collection : numpy.ndarray
Expand All @@ -363,7 +578,6 @@ def nb_constraint_rotational_values(
3D (dim, dim, blocksize) array containing data with `float` type.
indices : numpy.ndarray
1D array containing the index of constraining nodes

"""
block_size = indices.size
for i in range(block_size):
Expand All @@ -377,7 +591,6 @@ def nb_constrain_translational_values(
) -> None:
"""
Computes constrain values in numba njit decorator

Parameters
----------
position_collection : numpy.ndarray
Expand All @@ -386,7 +599,6 @@ def nb_constrain_translational_values(
2D (dim, blocksize) array containing data with `float` type.
indices : numpy.ndarray
1D array containing the index of constraining nodes

"""
block_size = indices.size
for i in range(block_size):
Expand All @@ -398,7 +610,6 @@ def nb_constrain_translational_values(
def nb_constrain_translational_rates(velocity_collection, indices) -> None:
"""
Compute constrain rates in numba njit decorator

Parameters
----------
velocity_collection : numpy.ndarray
Expand All @@ -419,7 +630,6 @@ def nb_constrain_translational_rates(velocity_collection, indices) -> None:
def nb_constrain_rotational_rates(omega_collection, indices) -> None:
"""
Compute constrain rates in numba njit decorator

Parameters
----------
omega_collection : numpy.ndarray
Expand Down
Loading