Skip to content

Commit

Permalink
Merge block feature for v0.2 (GazzolaLab#84)
Browse files Browse the repository at this point in the history
## Changes

* Rigid Body:
  * Rigid sphere and cylinder test cases
    * New test cases to test rotational and translational
  * Remove internal forces and torques in rigid body: not used
motion of rigid cylinder and sphere.
* Update gitignore for zip and csv file formats

## Enhancement
* Linalg numba faster functions for block structure.
* Factory function updated for new block structure.
* New wrapper to call memory block added
* Time-stepper interfaces are changed for block-structure implementation.
* Block-structure functions for Cosserat rod and reset ghosts.
* Memory block for the rigid body imp
* Memory block validity test cases added
* Governing equation tests for block structure

## Fix
* Rod intialization test fixed
* Stepper interface tests fixed for block structure
* Time-stepper analytic solution tests fixed
  * This commit fixes time-stepper tests using Harmonic oscillator test
cases for block structure implementation.
* Fix a bug in torque calculation. Transport term Jwxw was missing and with this commit we are adding this term.
* Fix: interaction routines for rigid body
* Fix: contact with rigid body
  * Due to the new implementation of rigid body block structure, radius and length are now are 1d arrays with 1 element, Before, they were scalar. This commit fixes the shape problems because of the array implementation.

## Note
This commit is the final version of rigid body block-structure implementation.
Issues and errors are fixed and code is tested.

## Order
Resolve GazzolaLab#9 

Co-authored-by: armantekinalp <[email protected]>
  • Loading branch information
skim0119 and armantekinalp authored Dec 21, 2021
1 parent 486cd5d commit 1305604
Show file tree
Hide file tree
Showing 62 changed files with 5,179 additions and 1,678 deletions.
6 changes: 6 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -229,3 +229,9 @@ outcmaes/*

# pickle files
*.pickle

# zip files
*.zip

# csv files
*.csv
6 changes: 6 additions & 0 deletions elastica/_calculus.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@
_average,
_clip_array,
_isnan_check,
_trapezoidal_for_block_structure,
_two_point_difference_for_block_structure,
)

position_difference_kernel = _difference
Expand All @@ -23,7 +25,11 @@
_clip_array,
_isnan_check,
_get_zero_array,
_trapezoidal_for_block_structure,
_two_point_difference_for_block_structure,
)

quadrature_kernel = _trapezoidal
difference_kernel = _two_point_difference
quadrature_kernel_for_block_structure = _trapezoidal_for_block_structure
difference_kernel_for_block_structure = _two_point_difference_for_block_structure
97 changes: 97 additions & 0 deletions elastica/_elastica_numba/_calculus.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@
from numpy import zeros, empty
import numba
from numba import njit
from elastica._elastica_numba._reset_functions_for_block_structure._reset_ghost_vector_or_scalar import (
_reset_vector_ghost,
)


@njit(cache=True)
Expand Down Expand Up @@ -48,6 +51,54 @@ def _trapezoidal(array_collection):
return temp_collection


@njit(cache=True)
def _trapezoidal_for_block_structure(array_collection, ghost_idx):
"""
Simple trapezoidal quadrature rule with zero at end-points, in a dimension agnostic way. This form
specifically for the block structure implementation and there is a reset function call, to reset
ghosts.
Parameters
----------
array_collection
ghost_idx
Returns
-------
Notes
-----
Micro benchmark results, for a block size of 100, using timeit
Python version: 8.21 µs per loop
This version: 1.03 µs per loop
User should use this function with extreme care, since this function is rewritten for
block structure.
"""

_reset_vector_ghost(array_collection, ghost_idx)

blocksize = array_collection.shape[1]
temp_collection = np.empty((3, blocksize + 1))

temp_collection[0, 0] = 0.5 * array_collection[0, 0]
temp_collection[1, 0] = 0.5 * array_collection[1, 0]
temp_collection[2, 0] = 0.5 * array_collection[2, 0]

temp_collection[0, blocksize] = 0.5 * array_collection[0, blocksize - 1]
temp_collection[1, blocksize] = 0.5 * array_collection[1, blocksize - 1]
temp_collection[2, blocksize] = 0.5 * array_collection[2, blocksize - 1]

for i in range(3):
for k in range(1, blocksize):
temp_collection[i, k] = 0.5 * (
array_collection[i, k] + array_collection[i, k - 1]
)

return temp_collection


@njit(cache=True)
def _two_point_difference(array_collection):
"""
Expand Down Expand Up @@ -83,6 +134,52 @@ def _two_point_difference(array_collection):
return temp_collection


@njit(cache=True)
def _two_point_difference_for_block_structure(array_collection, ghost_idx):
"""
This function does the differentiation, for Cosserat rod model equations. This form
specifically for the block structure implementation and there is a reset function call, to
reset ghosts.
Parameters
----------
array_collection
ghost_idx
Returns
-------
Note
----
Micro benchmark results showed that for a block size of 100, using timeit
Python version: 7.1 µs per loop
This version: 1.01 µs per loop
User should use this function with extreme care, since this function is rewritten for
block structure.
"""

_reset_vector_ghost(array_collection, ghost_idx)

blocksize = array_collection.shape[1]
temp_collection = np.empty((3, blocksize + 1))

temp_collection[0, 0] = array_collection[0, 0]
temp_collection[1, 0] = array_collection[1, 0]
temp_collection[2, 0] = array_collection[2, 0]

temp_collection[0, blocksize] = -array_collection[0, blocksize - 1]
temp_collection[1, blocksize] = -array_collection[1, blocksize - 1]
temp_collection[2, blocksize] = -array_collection[2, blocksize - 1]

for i in range(3):
for k in range(1, blocksize):
temp_collection[i, k] = array_collection[i, k] - array_collection[i, k - 1]

return temp_collection


@njit(cache=True)
def _difference(vector):
"""
Expand Down
Loading

0 comments on commit 1305604

Please sign in to comment.