File elasticRod.h

class elasticRod

Public Functions

elasticRod(int limb_idx, const Vector3d &start, const Vector3d &end, int num_nodes, double rho, double rod_radius, double youngs_modulus, double poisson_ratio, double mu)

Constructor for initializing a straight rod as defined by a start and end point.

Parameters:

limb_idx – The limb id of the rod. Should be a +1 increment of the previous rod. Other inputs will lead to undefined behavior.
start – The start point of the rod.
end – The end point of the rod.
num_nodes – The number of nodes the rod will be discretized to.
rho – Density of the rod [kg/m^3].
rod_radius – Radius of the rod [m].
youngs_modulus – Young’s modulus of the rod [N/m^2].
poisson_ratio – Poisson’s ratio of the rod.
mu – Friction coefficient of the rod.

elasticRod(int limb_idx, const vector<Vector3d> &nodes, double rho, double rod_radius, double youngs_modulus, double poisson_ratio, double mu)

Constructor for initializing an arbitrarily shaped rod as by a list of nodes.

Parameters:

limb_idx – The limb id of the rod. Should be a +1 increment of the previous rod. Other inputs will lead to undefined behavior.
nodes – A list of consecutive nodes of the rod.
rho – Density of the rod [kg/m^3].
rod_radius – Radius of the rod [m].
youngs_modulus – Young’s modulus of the rod [N/m^2].
poisson_ratio – Poisson’s ratio of the rod.
mu – Friction coefficient of the rod.

~elasticRod(): Default deconstructor.

void prepareForIteration()

Updates all necessary discrete values and frames for edges in preparation of next timestep.

Updates discrete tangents, reference frames, material frames, reference twists, edge lengths, and curvatures.

double updateNewtonX(double *dx, int offset, double alpha = 1.0)

Performs a Newton update iteration.

Parameters:

dx – The update vector provided by the time stepper and solver.
offset – A constant offset for the rod’s DOF index in dx with respect to the global system.
alpha – A Newton update coefficient.

Returns:

The max non-theta dx update.

void updateGuess(double weight, double dt)

Updates the current guess at the start of a Newton loop.

Updates all unconstrained DOFs using

x = x0 + weight + u * dt

Generally, performance seems best when the guess is close to x0, especially when friction is involved.

Parameters:

weight – The weight of the velocity update shown above.
dt – The timestep.

void enable2DSim() const

Enables 2D simulation along x-z plane.

Assumes all y and theta related DOFs to be constrained.

Vector3d getVertex(int k)

Get the a (3,) position vector of the kth node from current timestep.

Parameters:: k – Index of the node.
Returns:: The (3,) position vector.

Vector3d getPreVertex(int k)

Get the a (3,) position vector of the kth node from previous timestep.

Parameters:: k – Index of the node.
Returns:: The (3,) position vector.

Vector3d getVelocity(int k)

Get the a (3,) velocity vector of the kth node from current timestep.

Parameters:: k – Index of the node.
Returns:: The (3,) velocity vector.

int getIfConstrained(int k) const

Getter function for seeing if a certain DOF is constrained.

Parameters:: k – Index of the queried DOF.
Returns:: Value of 1 if constrained and 0 if unconstrained.

void updateMap()

Updates the constraint map according to current values of isConstrained and isDOFJoint.

Must be called when new boundary conditions are created or if boundary conditions are released. If the value of a preexisting boundary condition is updated, this function does not need to be called.

void addJoint(int node_num, bool attach_to_joint, int joint_node, int joint_limb)

Converts a node on the rod to a joint node.

This function should NOT be called by the user itself. The elasticJoint class will handle calling this function.

There are two types of applications to this function.

1) The node is being converted to a joint node that does not previously exist. To create a joint node for the first time, attach_to_joint must be set to False. In this setting, joint_node and joint_limb are ignored.

2) The node is being attached to a preexisting joint node. As of right now, only the first and last nodes can be converted in this fashion. To attach this node to preexisting joint node, attach_to_joint must be set to True and joint_node and joint_limb must be set to the desired preexisting joint node.

Parameters:

node_num – The index for the node to be converted.
attach_to_joint – Whether or not this node is being attached to a preexisting joint.
joint_node – The joint node id of the preexisting joint.
joint_limb – The joint limb id of the preexisting joint.

void setVertexBoundaryCondition(Vector3d position, int k)

void setThetaBoundaryCondition(double desired_theta, int k)

void freeVertexBoundaryCondition(int k)

Vector3d getTangent(int k)

double getTheta(int k)

Public Members

int limb_idx: The limb id of the rod.

double rho: Density [kg/m^3].

double rod_radius

Cross-sectional radius of the rod [m].

Currently, assumes constant radius.

double cross_sectional_area: Cross-sectional area of the rod [m^2].

double mu: Friction coefficient.

double rod_length: Total relaxed length of the rod [m].

double youngM: Young’s modulus.

double shearM

Shear modulus.

See also

Computed in computeElasticStiffness().

double poisson_ratio: Poisson ratio.

double EA

Stretching stiffness.

See also

Computed in computeElasticStiffness().

double EI

Bending stiffness.

See also

Computed in computeElasticStiffness().

double GJ

Twisting stiffness.

See also

Computed in computeElasticStiffness().

int nv: Number of vertices (nodes).

int ne

Number of edges.

Number of edges is always one less than the number of nodes:

ne = nv-1

int ndof

Number of degrees of freedom (DOF).

Number of DOFs is always

ndof = 3*nv + ne

int ncons: Number of constrained DOF.

int uncons: Number of unconstrained DOF.

VectorXd x0

DOF vector from the previous timestep.

Vector of size (ndof, 1).

VectorXd x

DOF vector during and after the current timestep.

Vector of size (ndof, 1).

VectorXd x_ls

DOF vector used to save state during line search.

Vector of size (ndof, 1).

VectorXd u0

Velocity vector from the previous timestep.

Vector of size (ndof, 1).

VectorXd u

Velocity vector during and after the current timestep.

Vector of size (ndof, 1).

VectorXd edge_len

Vector of current edge lengths [m].

See also

Computed in computeEdgeLen().

VectorXd ref_len

Reference lengths of the discrete edges [m].

Computed only once during initialization. Starting configuration is assumed to be resting configuration.

See also

Computed in setReferenceLength().

VectorXd voronoi_len

Voronoi lengths of each node [m].

Computed only once during initialization. Computed from reference lengths.

See also

Computed in setReferenceLength().

MatrixXd kb

Curvature binormals of each discrete edge.

Computed during every iteration / timestep to compute the discrete curvatures necessary to compute bending energy.

See also

Computed in computeKappa().

MatrixXd d1

d1 axes of the reference frame during current iteration / timestep.

Matrix of size (ne, 3).

See also

Updated in computeTimeParallel().

See also

Used in getRefTwist().

MatrixXd d2

d2 axes of the reference frame during current iteration / timestep.

Matrix of size (ne, 3).

See also

Updated in computeTimeParallel().

See also

Used in getRefTwist().

MatrixXd d1_old

d1 axes of the reference frame from last timestep.

Matrix of size (ne, 3). Updated by time stepper classes when completing a timestep.

See also

Used in computeTimeParallel().

MatrixXd d2_old

d2 axes of the reference frame from last timestep.

Matrix of size (ne, 3). Updated by time stepper classes when completing a timestep.

MatrixXd m1

m1 axes of the material frame during current iteration / timestep.

Matrix of size (ne, 3).

See also

Updated in computeMaterialDirector().

See also

Used in computeKappa() to compute discrete curvatures.

MatrixXd m2

m2 axes of the material frame during current iteration / timestep.

Matrix of size (ne, 3).

See also

Updated in computeMaterialDirector().

See also

Used in computeKappa() to compute discrete curvatures.

MatrixXd tangent

Edge tangents during the current iteration / timestep.

Matrix of size (ne, 3).

See also

Updated in computeTangent().

See also

Used in computeTimeParallel().

MatrixXd tangent_old

Edge tangents from the previous timestep.

Matrix of size (ne, 3).

See also

Used in computeTimeParallel().

VectorXd ref_twist

Reference twists during the current iteration / timestep.

Vector of size (ne, 1).

See also

Updated in getRefTwist().

See also

Used in computeTwistBar().

VectorXd ref_twist_old

Reference twists from the previous timestep.

Vector of size (ne, 1).

VectorXd twist_bar

The twists in the resting configuration.

Vector of size (ne, 1).

See also

Updated in computeTwistBar().

VectorXd mass_array

Lumped mass array.

Vector of size (ndof, 1). Stores the lumped masses of each DOF.

See also

Updated in setMass().

MatrixXd kappa

Discrete curvatures of the inner nodes.

Matrix of size (nv, 2). Note that the first and outer values are not used since curvatures only apply to nodes with two edges attached.

See also

Updated in computeKappa().

MatrixXd kappa_bar

Discrete curvatures of the inner nodes in the resting configuration.

Matrix of size (nv, 2). Note that the first and outer values are not used since curvatures only apply to nodes with two edges attached. Is assumed to the curvatures computed from the initial starting configuration.

int *isConstrained

Array storing whether or not a certain DOF is constrained.

Size of ndof. Value of 1 if constrained and 0 if unconstrained.

int *unconstrainedMap

Array that maps unconstrained to full indices.

Size of uncons. Maps unconstrained DOF (by array location) to their respective index in x vector. The index of the mapping is stored.

See also

Gets updated in updateMap() and in setup().

See also

The inverse mapping is fullToUnconsMap.

int *fullToUnconsMap

Array that maps full to unconstrained indices.

Size of ndof. Maps all DOF from x vector (by array location) to their respective in unconstrained vector. The index of the mapping is stored. If a certain DOF is constrained, then the value at that array location is undefined.

See also

Gets updated in updateMap() and in setup().

See also

The inverse mapping is unconstrainedMap.

int *isDOFJoint

Array indicating whether or not a certain DOF is a joint node.

Size of ndof. Can be one of three possible values: 0: DOF is not a joint affiliated DOF. 1: DOF is a joint affiliated DOF and was attached to a preexisting joint. 2: DOF is a joint affiliated DOF and was used to create a joint node for the first time.

See also

Updated in addJoint().

int *isNodeJoint

Array indicating whether or not a certain node is a joint node.

Size of nv. Follows the same rules as isDOFJoint.

See also

Updated in addJoint().

int *isEdgeJoint

Array indicating whether or not a certain edge is adjacent to a joint node.

Size of ne. Follows the same rules as isDOFJoint.

See also

Updated in addJoint().

vector<pair<int, int>> joint_ids

Vector of joint node and limb ids.

Stores a list of pairs of joint node and joint limb ids. Used to reroute forces and Jacobians to their proper locations in the global force vector / Jacobian matrix for nodes that were attached to preexisting joint nodes. Currently of size nv.

See also

Updated in addJoint().

Private Functions

void setupMap()

void setup(const vector<Vector3d> &nodes)

void computeElasticStiffness()

void setMass()

void setReferenceLength()

void computeTimeParallel()

void computeTangent()

void computeSpaceParallel()

void computeMaterialDirector()

void computeKappa()

void computeTwistBar()

void computeEdgeLen()

void getRefTwist()

Private Static Functions

static void parallelTransport(const Vector3d &d1_1, const Vector3d &t1, const Vector3d &t2, Vector3d &d1_2)

static void rotateAxisAngle(Vector3d &v, const Vector3d &z, const double &theta)

static double signedAngle(const Vector3d &u, const Vector3d &v, const Vector3d &n)