Converting nodal velocities to element velocities

[armantekinalp]

In PyElastica we store velocity on nodes and for some computations (i.e. rod damping, plane damping or stokes flow) we have to convert velocity from nodes to elements. Currently we are averaging the neigboring node velocity to compute element velocity using the following function.

PyElastica/elastica/interaction.py

Lines 664 to 703 in fee87b9

	def node_to_element_pos_or_vel(vector_in_node_frame):
	"""
	This function computes the velocity of the elements.
	Here we define a separate function because benchmark results
	showed that using Numba, we get more than 3 times faster calculation.
	Parameters
	----------
	vector_in_node_frame: numpy.ndarray
	2D (dim, blocksize) array containing data with 'float' type.
	Returns
	-------
	vector_in_element_frame: numpy.ndarray
	2D (dim, blocksize) array containing data with 'float' type.
	"""
	"""
	Developer Notes
	-----
	Benchmark results, for a blocksize of 100,
	Python version: 3.5 µs ± 149 ns per loop
	This version: 729 ns ± 14.3 ns per loop
	"""
	n_elem = vector_in_node_frame.shape[1] - 1
	vector_in_element_frame = np.empty((3, n_elem))
	for k in range(n_elem):
	vector_in_element_frame[0, k] = 0.5 * (
	vector_in_node_frame[0, k + 1] + vector_in_node_frame[0, k]
	)
	vector_in_element_frame[1, k] = 0.5 * (
	vector_in_node_frame[1, k + 1] + vector_in_node_frame[1, k]
	)
	vector_in_element_frame[2, k] = 0.5 * (
	vector_in_node_frame[2, k + 1] + vector_in_node_frame[2, k]
	)
	return vector_in_element_frame

Although this function is correct for uniform rods (uniform cross-section, element length and density) for tapered rods, momentum is not conserved since neighboring node masses are not the same. Thus, we have to first compute the total momentum of element and divide by the element mass to compute element velocity.

How to fix :

1. Use the following function to compute element velocities

def node_to_element_velocity(mass, velocity_collection):
    """
    """
    n_elem = velocity_collection.shape[1] - 1
    element_velocity = np.empty((3, n_elem))
    for k in range(n_elem):
        element_velocity[0, k] = (
            mass[k+1] * velocity_collection[0, k + 1] + mass[k] * velocity_collection[0, k]
        )
        element_velocity[1, k] = (
            mass[k+1] * velocity_collection[1, k + 1] + mass[k] * velocity_collection[1, k]
        )
        element_velocity[2, k] =(
            mass[k+1] * velocity_collection[2, k + 1] + mass[k] * velocity_collection[2, k]
        )
        element_velocity[:,k] /= (mass[k+1] + mass[k])


    return element_velocity

2. Rename the function node_to_element_pos_or_vel as node_to_element_position, since this will be only used to compute the element positions.

[bhosale2]

Closed by #98.