12. Element and Node Utilities - ElemTran, NodeGeom, Corot

The element and node utilities perform useful auxiliary functions such as transformation and condensation of element degrees of freedom and accumulation of forces and moments about a point, etc. The ElemTran module performs operations on element matrices and vectors, the NodeGeom module performs operations which are a function of node coordinates. The Corot module provides a methodology for augmenting stiffness and reaction formations for infinitesimally-based finite elements to allow for large rotations.

12.1. Element Transformations - ElemTran

The ElemTran module is designed to transform and/or condense the coefficients of an element vector or matrix. Transformations are performed by specifying direction cosine matrices selectively at each node to be transformed. This capability is useful for enabling the application of skewed boundary conditions. Condensations and expansions (uncondensations) are performed by setting a degree of freedom condensation “map”. The condensation map specifies which element degrees of freedom are to be condensed out. The condensation feature can be used to remove internal element degrees of freedom or to enable pin connections in elements with rotational degrees of freedom.

The methods associated with a ElemTran object are the following.

Instance a ElemTran object using vfe_ElemTranBegin(). Once a ElemTran is instanced, set the current element degree of freedom map using vfe_ElemTranSetMap().

Transform selected nodal freedoms of an element vector using vfe_ElemTranVector(). Transform selected nodal freedoms of an element matrix using vfe_ElemTranMatrix(). The transformation matrices are given by a set of direction cosine matrices indicating the coordinate system to which nodal freedoms are to be transformed.

Selected degrees of freedom may be condensed from element matrices and vectors. The degrees of freedom to be condensed are set using vfe_ElemTranSetCond(). By default, this is set to condense all internal element degrees of freedom. Use vfe_ElemTranCondVector() and vfe_ElemTranCondMatrix() to condense element vectors and matrices. Use the function vfe_ElemTranCondMap() to query for the degree of freedom map which reflects the degree of freedom type and ordering for the condensed element vectors and matrices.

Destroy an instance of a ElemTran object using vfe_ElemTranEnd().

12.2. Direction Cosine Matrix

Given that X’,Y’ and Z’ are three orthonormal vectors indicating the direction of the local coordinate axes in the global coordinate system (x,y,z), then the direction cosine matrix, tm for this local coordinate system is defined as:

tm[0][0] = X'x  tm[0][1] = X'y  tm[0][2] = X'z
tm[1][0] = Y'x  tm[1][1] = Y'y  tm[1][2] = Y'z
tm[2][0] = Z'x  tm[2][1] = Z'y  tm[2][2] = Z'z

12.3. Element Degree of Freedom Condensation and Elimination

Element degree of freedom condensation is a technique for statically condensing out a set of degrees of freedom within an element which do not couple with the degrees of freedom of adjacent elements. Typically condensation is used to remove internal element degrees which do not occur at nodes and therefore do not couple with adjacent elements. In this case element condensation of these internal degrees of freedom may be performed for economic reasons. Another use of element condensation is to eliminate nodal degrees of freedom from an element so that the forces associated with the element freedom do not couple or transmit to adjacent elements. For example the rotational degrees of freedom at the end of a beam element may be condensed resulting in a “pin” connection.

The process of element matrix condensation produces two associated recovery matrices which must be used as input to condense an associated element vector. In turn the process of element vector condensation produces an additional recovery vector which must be used along with the recovery matrices for uncondensing an element vector.

Assume that the uncondensed element matrix S is represented below in which the matrix has been partitioned into the retained freedoms followed by the condensed freedoms.

     --         --
     |  A     B  |
S =  |           |
     |  B(T)  C  |
     --         --

The condensed matrix, SCOND, is given by

SCOND = A - B(T) * C(-1) * B

The matrices B and C must be retained for performing any subsequent element vector condensations using the functions vfe_ElemTranCondVector() and vfe_ElemTranUncondVector() where they are referred to as the off diagonal recovery matrix bm and the block diagonal recovery matrix cm.

Assume that N is the number of degrees of freeom in the uncondensed matrix. This number is determined by the degree of freedom map specified using vfe_ElemTranSetMap(). Assume that M is the number of degrees of freeom in the condensed matrix. This number may be queried using vfe_ElemTranCondMap(). Let NC = N-M be the number of degrees of freedom which are condensed. The matrix B is stored column-wise and contains NC*M entries. The matrix C is stored in triangular format, already factorized:

C = L D L(T)

with the diagonal entries of L equal to 1. C is returned in cm and contains NC*(NC+1)/2 entries. The exact data structures used for B and C are not exposed.

In some situations degrees of freedom do not need to be statically condensed but simply struck from the element vector or matrix. This process is performed using vfe_ElemTranElimVector() and vfe_ElemTranElimMatrix(). For example, if the internal degrees of freedom of an element stiffness matrix are condensed, the internal freedoms of the associated element mass matrix only need to be eliminated.

12.4. Function Descriptions

The currently available ElemTran functions are described in detail in this section.

vfe_ElemTran *vfe_ElemTranBegin(void)

create an instance of an ElemTran object

Create an instance of a ElemTran object. Memory is allocated for the object private data and the pointer to the object is returned. The default condensation map is set to condense internal element degrees of freedom.

Destroy an instance of a ElemTran object using 

void vfe_ElemTranEnd (vfe_ElemTran *elemtran )

Return the current value of a ElemTran object error flag using 

Vint vfe_ElemTranError (vfe_ElemTran *elemtran )

Returns:

The function returns a pointer to the newly created ElemTran object. If the object creation fails, NULL is returned.

void vfe_ElemTranEnd(vfe_ElemTran *p)

destroy an instance of an ElemTran object

See vfe_ElemTranBegin()

Vint vfe_ElemTranError(vfe_ElemTran *p)

return the current value of an ElemTran object error flag

See vfe_ElemTranBegin()

void vfe_ElemTranCondVector(vfe_ElemTran *p, Vdouble v[], Vdouble bm[], Vdouble cm[], Vdouble vcond[], Vdouble bv[])

condense an element vector

Condense an element vector, v, and return the condensed vector in vcond. A recovery vector, bv, is also returned which is required if condensed element vectors are to be uncondensed using vfe_ElemTranUncondVector(). The element vector condensation procedure is preformed relative to an element matrix such as an element stiffness matrix. The recovery matrices bm and cm appearing as input are the recovery matrices returned by vfe_ElemTranCondMatrix() which must have been previously called to condense the element matrix associated with the element vector v.

Errors

SYS_ERROR_OPERATION is generated if the number of degrees of freedom in the condensed element is zero.

Parameters:

  • p – Pointer to ElemTran object.

  • v – Element vector to condense

  • bm – Off diagonal recovery matrix

  • cm – Block Diagonal recovery matrix

  • vcond – [out] Condensed element vector

  • bv – [out] Recovery vector

void vfe_ElemTranElimVector(vfe_ElemTran *p, Vdouble v[], Vdouble velim[])

eliminate freedoms from an element vector

Eliminate degrees of freedom from an element vector, v, and return the eliminated vector in velim. The rows associated with the freedoms specified by vfe_ElemTranSetCond() are struck from the vector.

Errors

SYS_ERROR_OPERATION is generated if the number of degrees of freedom in the condensed element is zero.

Parameters:

  • p – Pointer to ElemTran object.

  • v – Element vector to eliminate.

  • velim – [out] Eliminated element vector

void vfe_ElemTranUncondVector(vfe_ElemTran *p, Vdouble vcond[], Vdouble bm[], Vdouble cm[], Vdouble bv[], Vdouble v[])

expand (uncondense) an element vector

Uncondense an element vector, vcond, and return the expanded vector in v. The recovery matrices bm and cm appearing as input are the recovery matrices returned by vfe_ElemTranCondMatrix() which must have been previously called to condense the element matrix associated with the element vector v. The recovery vector, bv, is returned by a previous call to vfe_ElemTranCondVector().

Errors

SYS_ERROR_OPERATION is generated if the number of degrees of freedom in the condensed element is zero.

Parameters:

  • p – Pointer to ElemTran object.

  • vcond – Element vector to uncondense.

  • bm – Off diagonal recovery matrix

  • cm – Block Diagonal recovery matrix

  • bv – Recovery vector

  • v – [out] Uncondensed element vector

void vfe_ElemTranUnelimVector(vfe_ElemTran *p, Vdouble ul[], Vdouble vl[])

expand an element vector

In the un-elimination process the array of uneliminated degrees of freedom will match the input except for inactive degrees of freedom, whose value will be set to zero. Both input and output arrays have the same dimension and map to the same degrees of freedom.

Parameters:

  • p – Pointer to ElemTran object.

  • ul – Array of degrees of freedom

  • vl – [out] Uneliminated degrees of freedom

void vfe_ElemTranCondMatrix(vfe_ElemTran *p, Vdouble s[], Vdouble scond[], Vdouble bm[], Vdouble cm[])

condense an element matrix

Condense an element matrix, s, and return the condensed matrix in scond. Two recovery matrices, bm and cm, are returned which are required to condense or uncondense any associated element vectors.

Errors

SYS_ERROR_OPERATION is generated if the number of degrees of freedom in the condensed element is zero.

Parameters:

  • p – Pointer to ElemTran object.

  • s – Element matrix to condense.

  • scond – [out] Condensed element matrix

  • bm – [out] Off diagonal recovery matrix

  • cm – [out] Block Diagonal recovery matrix

void vfe_ElemTranElimMatrix(vfe_ElemTran *p, Vdouble s[], Vdouble selim[])

eliminate freedoms from an element matrix

Eliminate degrees of freedom from an element matrix, s, and return the eliminated matrix in selim. The rows and columns associated with the freedoms specified by vfe_ElemTranSetCond() are struck from the matrix.

Errors

SYS_ERROR_OPERATION is generated if the number of degrees of freedom in the condensed element is zero.

Parameters:

  • p – Pointer to ElemTran object.

  • s – Element matrix to eliminate.

  • selim – [out] Eliminated element matrix

void vfe_ElemTranSetMap(vfe_ElemTran *p, Vint nedofs, Vint loc[], Vint tag[])

set current element degree of freedom map

Set the current number of element degrees of freedom and the degree of freedom map. The loc vector consists of nedofs one-based node indices into the element connectivity. A zero entry in loc denotes an internal elemental degree of freedom. The tag vector consists of nedofs defined constants which specify the degree of freedom type, eg. x-translation, y-rotation, etc. types.

Errors

  • SYS_ERROR_VALUE is generated if an improper nedofs or loc is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to ElemTran object.

  • nedofs – Number of element degrees of freedom

  • loc – Vector of degree of freedom locations

  • tag – Vector of degree of freedom types

void vfe_ElemTranSetCond(vfe_ElemTran *p, Vint nconds, Vint loc[], Vint tag[])

set degrees of freedom to be condensed

Set the current degree of freedom map of degrees of freedom to be condensed or eliminated from element vectors and matrices. The degree of freedom map representing the degree of freedom location and type of the condensed vectors and matrices may be queried using vfe_ElemTranCondMap(). By default element internal degrees of freedom are condensed. This is equivalent to nconds = 1, loc[0] = 0 and tag[0] = SYS_DOF_ELEM.

Errors

  • SYS_ERROR_VALUE is generated if an improper nconds or loc is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to ElemTran object.

  • nconds – Number of element degrees of freedom to condense

  • loc – Vector of degree of freedom locations

  • tag – Vector of degree of freedom types

void vfe_ElemTranCondMap(vfe_ElemTran *p, Vint *nedofs, Vint loc[], Vint tag[])

query condensed degree of freedom map

Query for the degree of freedom map of condensed element vectors and matrices. This degree of freedom map is the result of eliminating the condensed degrees of freedom specified by vfe_ElemTranSetCond() from the input element degree of freedom map specified by vfe_ElemTranSetMap().

Parameters:

  • p – Pointer to ElemTran object.

  • nedofs – [out] Number of element degrees of freedom in condensed element

  • loc – [out] Vector of degree of freedom locations

  • tag – [out] Vector of degree of freedom types

void vfe_ElemTranMatrix(vfe_ElemTran *p, Vint nin, Vint in[], Vdouble tm[][3][3], Vdouble s[])

transform an element matrix

Transform a selected set, in of freedoms at nin nodes of an element matrix, s. The element matrix is transformed in place. The transformation is specified by an array of nin direction cosine matrices.

Errors

SYS_ERROR_VALUE is generated if an improper nin is specified.

Parameters:

  • p – Pointer to ElemTran object.

  • nin – Number of nodes to transform

  • in – Vector of nin node indices to element connectivity

  • tm – Array of nin direction cosine matrices indicating the coordinate system to which nodal freedoms are to be transformed

  • s – Element matrix to transform.

void vfe_ElemTranVector(vfe_ElemTran *p, Vint nin, Vint in[], Vdouble tm[][3][3], Vdouble v[])

transform an element vector

Transform a selected set, in of freedoms at nin nodes of an element vector, v. The element vector is transformed in place. The transformation is specified by an array of nin direction cosine matrices.

Errors

SYS_ERROR_VALUE is generated if an improper nin is specified.

Parameters:

  • p – Pointer to ElemTran object.

  • nin – Number of nodes to transform

  • in – Vector of nin node indices to element connectivity

  • tm – Array of nin direction cosine matrices indicating the coordinate system to which nodal freedoms are to be transformed

  • v – Element vector to transform.

vfe_NodeGeom *vfe_NodeGeomBegin(void)

create an instance of a NodeGeom object

Create an instance of a NodeGeom object. Memory is allocated for the object private data and the pointer to the object is returned. The default center is (0.,0.,0.) and all motions are zero.

Destroy an instance of a NodeGeom object using 

void vfe_NodeGeomEnd (vfe_NodeGeom *nodegeom )

Return the current value of a NodeGeom object error flag using 

Vint vfe_NodeGeomError (vfe_NodeGeom *nodegeom )

Returns:

The function returns a pointer to the newly created NodeGeom object. If the object creation fails, NULL is returned.

void vfe_NodeGeomEnd(vfe_NodeGeom *p)

destroy an instance of a NodeGeom object

See vfe_NodeGeomBegin()

Vint vfe_NodeGeomError(vfe_NodeGeom *p)

return the current value of a NodeGeom object error flag

See vfe_NodeGeomBegin()

void vfe_NodeGeomSetCenter(vfe_NodeGeom *p, Vdouble xc[3])

set center point

Set the current center point for rigid body calculation, force and moment summation and nodal acceleration vector calculation. By default the center point is (0.,0.,0.).

Parameters:

  • p – Pointer to NodeGeom object.

  • xc – Coordinates of center point

void vfe_NodeGeomSetMotion(vfe_NodeGeom *p, Vdouble a[3], Vdouble w[3], Vdouble wdot[3])

specify acceleration

Set acceleration state description. The translational acceleration vector, a, should contain the sum of all translational acceleration fields including gravity. The angular velocity vector, w and angular acceleration vector, wdot are assumed be about the current center as specified using vfe_NodeGeomSetCenter(). By default a, w and wdot are set to zero.

Parameters:

  • p – Pointer to NodeGeom object.

  • a – Translational acceleration vector

  • w – Angular velocity of rotating reference frame

  • wdot – Angular acceleration of rotating reference frame

void vfe_NodeGeomRigid(vfe_NodeGeom *p, Vint irigid, Vdouble x[3], Vint ntags, Vint tag[], Vdouble u[])

compute rigid body vector

Compute the rigid body motion of a point at x. The rigid body mode type is given by irigid. The number of degrees of freedom and the nodal degree of freedom types are input using ntags and tag. The vector u contains the output ntags rigid body motions. The point about which rotational rigid body modes are calculated is specified using vfe_NodeGeomSetCenter().

Errors

  • SYS_ERROR_VALUE is generated if an improper irigid or ntags is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to NodeGeom object.

  • irigid – Rigid body vector to compute 

    =1  Translation in X
=2  Translation in Y
=3  Translation in Z
=4  Rotation about X
=5  Rotation about Y
=6  Rotation about Z

  • x – Node coordinates

  • ntags – Number of degrees of freedom types used at node

  • tag – Vector of degree of freedom types at node

  • u – [out] Rigid body motion for degrees of freedom

void vfe_NodeGeomForceMoment(vfe_NodeGeom *p, Vdouble x[3], Vint ntags, Vint tag[], Vdouble f[], Vdouble force[3], Vdouble moment[3])

sum force and moment about center

Sum the input degree of freedom forces, f, about the current center to produce a resultant force and moment vector.

Errors

  • SYS_ERROR_VALUE is generated if an improper ntags is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to NodeGeom object.

  • x – Node coordinates

  • ntags – Number of degrees of freedom types used at nodes

  • tag – Vector of degree of freedom types at nodes

  • f – Force and moment at node degrees of freedom

  • force – [out] Force vector

  • moment – [out] Moment vector about center

void vfe_NodeGeomAcceleration(vfe_NodeGeom *p, Vint npts, Vdouble x[][3], Vdouble a[][3])

compute nodal acceleration

Compute the acceleration vectors, a, of a set of npts at locations x, The acceleration state parameters are set using vfe_NodeGeomSetMotion() and the center of rotation is set using vfe_NodeGeomSetCenter().

Errors

SYS_ERROR_VALUE is generated if an improper npts is specified.

Parameters:

  • p – Pointer to NodeGeom object.

  • npts – Number of node coordinates

  • x – Node coordinates

  • a – [out] Acceleration vectors at nodes

void vfe_NodeGeomFreeMode(vfe_NodeGeom *p, Vdouble x[][3], Vint tmflag, Vdouble tm[][3][3], Vint ntags, Vint loc[], Vint tag[], Vint *ntra, Vdouble tra[3][3], Vint *nrot, Vdouble rot[3][6], Vdouble off[3][3])

compute unconstrained modes

Compute any possible translational and rotational rigid body modes existing under a set of constraints. The unique node coordinate locations, x, of the applied constraints are referenced by the ntags entries in the loc vector. The associated degree of freedom types of each constraint appear in the tag vector. The degree of freedom types are restricted to translations and rotations. The number of constraints, ntags, must be greater than 1.

If the degrees of freedom are not applied in the global Cartesian coordinate system then the tmflag must be enabled and a direction cosine matrix ,tm, input for each unique node. If the tmflag is disabled, then tm may be input as NULL.

The number of translational, ntra, and rotational, nrot, rigid body modes are computed and returned. If there are translational rigid body modes, the translations associated with each mode are returned in tra. If there are rotational rigid body modes, the translations and rotations associated with each mode are returned in rot. The center point about which rotational rigid body modes are calculated is specified using vfe_NodeGeomSetCenter(). The offset vector of a point relative to the center point through which the axis of rotation passes is returned in off. Note that the free modes are computed assuming the constraints are applied to a single continuous body.

Errors

  • SYS_ERROR_VALUE is generated if an improper ntags or loc is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to NodeGeom object.

  • x – Node coordinates

  • tmflag – Direction cosine matrices flag

  • tm – Node direction cosine matrices

  • ntags – Number of constrainted degrees of freedom

  • loc – Locations of constrained degrees of freedom

  • tag – Degree of freedom types of constraints

  • ntra – [out] Number of translational free modes

  • tra – [out] Vector of translational free mode translations

  • nrot – [out] Number of rotational free modes

  • rot – [out] Vector of rotational free mode translations and rotations

  • off – [out] Offset of axis of rotation from center

12.5. Node Geometry Utilities - NodeGeom

The NodeGeom module is designed to perform useful operations which are a function of the nodal point coordinates. These operations include computing rigid body modes, summing forces and moments, computing nodal acceleration vectors given information about a rotating frame of reference, etc.

The methods associated with a NodeGeom object are the following.

Instance a NodeGeom object using vfe_NodeGeomBegin(). Once a NodeGeom is instanced, set the center point about which rigid body modes are calculated, forces and moments are accumulated, etc. using vfe_NodeGeomSetCenter(). If nodal accelerations are to be computed using vfe_NodeGeomAcceleration(), then set the motion of the rotating reference frame and translational acceleration vector using vfe_NodeGeomSetMotion().

Rigid body modes may be calculated using vfe_NodeGeomRigid(). Forces and moments may be summed about the specified center point using vfe_NodeGeomForceMoment(). The function vfe_NodeGeomFreeMode() computes possible unconstrained motions or free modes given a set of translational and rotational constraints and their locations. This function is useful for determining if a set of user applied boundary conditions are sufficient to prevent rigid body motion and identify the rigid body motions if they exist.

12.6. Rigid Body Motion, Accelerations and Inertia Relief

A useful application of the NodeGeom module is in the use of the inertia relief technique to perform static analysis of unrestrained bodies under applied forces and moments. The basic notion of inertia relief is to compute a set of inertia forces, arising from the product of the mass matrix and nodal accelerations, which exactly balance the applied forces and moments. The nodal accelerations are assumed to be proportional to some linear combination of admissible rigid body motions, generally 3 translations and 3 rotations. The heart of the inertia relief calculation is to compute the 6 coefficients associated with the rigid body modes. The coefficients of each of these modes is chosen such as to exactly balance the applied loads. Certain functions in the NodeGeom have been designed to facilitate this calculation. The specific steps involved in the inertia relief calculation are outlined below:

  • Compute the 3 components of total force, Fext, and 3 components of moment, Mext of all the applied forces, Fapplied, about a reference point. Any point may serve as a suitable reference point, usually the origin (0.,0.,0.) of the global coordinate system is used. Use vfe_NodeGeomSetCenter() to set the center about which forces and moments are to be computed. Use vfe_NodeGeomForceMoment() to compute the force and moment about the specified center of a set of applied forces and moments at an arbitrary point.

  • Compute the rigid body displacements for each of 6 rigid body modes for each node in the finite element model. These rigid body displacements may be calculated using vfe_NodeGeomRigid() using the same center point used to accumulate the force and moment due to the applied forces above (a SysVector object in the VfsTools library is a convenient data structure to hold the rigid body displacements for all nodes for a particular rigid body mode).

  • Multiply each rigid body mode, by the nodal mass matrix (a SysMatrix object in the VfsTools library may be used to hold the mass matrix), to yield inertial forces associated with each mode, fx, fy, fz, frx, fry, frz, where fx are the forces due to translational acceleration in x, fry are the forces due to rotational acceleration about y, etc. The nodal forces of inertia, Finertia which are to exactly balance the applied force are a linear combination of these 6 vectors.

    Finertia = cx*fx + cy*fy + cz*fz + crx*frx + cry*fry + crz*frz

    Compute the total force and moment due to each rigid body mode about the center using vfe_NodeGeomForceMoment() and set this equal to the corresponding component of the total force and moment due to the applied load. Fxx is total force in x of fx, etc. This yields the following set of 6 equations and 6 unknowns to be solved for the coefficients cx, cy, etc.

    --     --     --                              --  --   --
| Fextx |     | Fxx  0    0    0    Fryx  Frzx |  | cx  |
| Fexty |     | 0    Fyy  0    Frxy 0     Frzy |  | cy  |
| Fextz |     | 0    0    Fzz  Frxz Fryz  0    |  | cz  |
| Mextx |     | 0    Myx  Mzx  Mxx  0     0    |  | crx |
| Mexty |     | Mxy  0    Mzy  0    Myy   0    |  | cry |
| Mextz |     | Mxz  Myz  0    0    0     Mzz  |  | crz |
--     --     --                              --  --   --

  • The inertial loads, Finertia, may now be computed and subtracted from the applied loads, Fapplied to yield the total loads to be applied to the body. Since the total load is now equilibrated, the body may now be constrained by any set of statically determinate constraints and no reaction forces will be felt by any of the constraints.

Another related application of NodeGeom is in the computation of nodal accelerations due to applied translational accelerations, rotational velocities and rotational accelerations. Set the center of rotation using vfe_NodeGeomSetCenter() and set the associated inertial state using vfe_NodeGeomSetMotion(). Then the nodal accelerations for any node location may be computed using vfe_NodeGeomAcceleration().

12.7. Function Descriptions

The currently available NodeGeom functions are described in detail in this section.

12.8. Corotational Utilities - Corot

The Corot module provides a convenient procedure which may be used to extend the capabilities of an existing linear finite element to allow for finite rotations. The procedure is an element independent corotational algorithm in which the element rigid body motion, consisting of translation and rotation, is separated from the total motion of the element leaving only the motion due to deformation. The Corot module may be applied to solid, shell, beam and truss finite elements with or without rotational degrees of freedom. It may also be applied to any other more abstract element for which a corotating frame may be defined, provided its topology has a dimensionality at least equal than a line.

The methods associated with a Corot object are the following.

Instance a Corot object using vfe_CorotBegin(). Once a Corot object is instanced, you must use vfe_CorotDef() to define the problem’s dimensionality. For two-dimensional problems, Corot will ensure that the corotational transformation always sets the third axis as the original z-axis (i.e., out of the plane, towards the reader). Set the topology of the element type for which a corotational frame is meant to be computed using vfe_CorotSetTopology(). Set the degree of freedom map for the element type using vfe_CorotSetMap(). At this point the state of the Corot object is ready to compute corotational frames and operate upon element stiffness and force vectors.

For each element there is, in general, a unique corotational frame which is attached to the element in its initial undeformed configuration and in its subsequent deformed configurations. At each iteration, use vfe_CorotPull() to pull-back the array of displacements and rotations degrees of freedom to an equivalent infinitesimal value applicable to the element’s initial reference configuration. vfe_CorotPull() will overwrite the incoming degrees of freedom with the new, infinitesimal values.

Once element matrices and vectors are computed using the original, undeformed element coordinates - an important point, since anisotropy in the element is defined in the reference configuration - they must be pushed forward to the current configuration using vfe_CorotPush() for assembly and solution, etc.

If the element formulation contains rotational degrees of freedom, there is an additional complication when large rotations are present. The update of the total displacement with the incremental displacement is not just a simple matter of addition when rotational degrees of freedom are present. This is due to the fact that the rotations are no longer, in general, infinitesimal. Use the function vfe_CorotUpdateRot() to add the increments of rotational degrees of freedom to the total rotation degree of freedom vector.

The function vfe_CorotSetParami() may be used to turn on and off the addition of the geometric “load” stiffness.

The user may want to toggle the calculation of the geometric load stiffness for nonlinear iteration convergence reasons. Disabling the geometric load stiffness term may provide a more stable but slower convergence rate.

12.9. Element Independent Corotation Frames

The obvious use of the Corot module is in the isolation of the effects of large rigid body motion, primarily rotation, from the underlying element formulation so that the analysis of large rotation phenomena may be modeled using finite elements which do not inherently support large rotations. These analyses are generally nonlinear and the corotational utilities must be applied to each element at each step in the nonlinear solution sequence. The following code fragment illustrates the basic procedure for a 4 node quadrilateral shell element.

                /* declare objects */
vfe_Corot *corot;
Vint nedofs, loc[24], tag[24];
Vdouble x[4][3], u[24];
Vdouble r[24], k[24*(24+1)/2];
Vdouble tm[3][3];
Vdouble ur[3], urd[3], urp[3];

                /* create corot object */
corot = vfe_CorotBegin ();

                /* define dimensionality */
vfe_CorotDef (corot,VFE_3D);

                /* specify element topology */
vfe_CorotSetTopology (corot,VIS_SHAPE_QUAD,2,0,0);

                /* set degree of freedom map */
vfe_CorotSetMap (corot,nedofs,loc,tag);

                /* at each load or time step */
...
                /* pull-back displacements to reference configuration */
vfe_CorotPull (corot,x,u,tm);

                /* form stiffness and force using modified u */
...

                /* push reactions and stiffness to current configuration */
vfe_CorotPush (corot,SYS_ON,r,SYS_ON,k);

                /* assemble forces, stiffnesses */
...
                /* solve for increment of displacement */
...
                /* for each rotation degree of freedom vector at a node */
                /* update rotation pseudovector */
vfe_CorotUpdateRot (corot,ur,urd,urp);

                /* go to next load or time step */
...
                /* delete corot object */
vfe_CorotEnd (corot);

12.10. Function Descriptions

The currently available Corot functions are described in detail in this section.

vfe_Corot *vfe_CorotBegin(void)

create an instance of a Corot object

Create an instance of a Corot object. Memory is allocated for the object private data and the pointer to the object is returned. The default topology is a linear quadrilateral.

Destroy an instance of a Corot object using 

void vfe_CorotEnd (vfe_Corot *corot )

Return the current value of a Corot object error flag using 

Vint vfe_CorotError (vfe_Corot *corot )

Returns:

The function returns a pointer to the newly created Corot object. If the object creation fails, NULL is returned.

void vfe_CorotEnd(vfe_Corot *p)

destroy an instance of a Corot object

See vfe_CorotBegin()

Vint vfe_CorotError(vfe_Corot *p)

return the current value of a Corot object error flag

See vfe_CorotBegin()

void vfe_CorotDef(vfe_Corot *p, Vint dim)

Define dimensionality.

Specifies the problem dimensionality. For two-dimensional problems, the corotational z-direction will always be equal to the original z-direction.

Errors

SYS_ERROR_ENUM is generated if an invalid value of dim is set.

Parameters:

  • p – Pointer to Corot object.

  • dim – Dimensionality 

    x=VFE_2D                 Two-dimensional problem
 =VFE_3D                 Three-dimensional problem

void vfe_CorotSetTopology(vfe_Corot *p, Vint shape, Vint maxi, Vint maxj, Vint maxk)

set element topology

Specify the element topology. The default topology is SYS_SHAPEQUAD with maxi = maxj = maxk = 0.

Errors

  • SYS_ERROR_VALUE is generated if an improper maxi, maxj or maxk is specified, or if a three-dimensional shape is specified for a two-dimensional problem.

  • SYS_ERROR_ENUM is generated if an improper shape is input.

Parameters:

  • p – Pointer to Corot object.

  • shape – The topological shape for the element 

    x=VIS_SHAPELINE          Line
 =VIS_SHAPETRI           Triangle
 =VIS_SHAPEQUAD          Quadrilateral
 =VIS_SHAPETET           Tetrahedron
 =VIS_SHAPEWED           Wedge
 =VIS_SHAPEHEX           Hexahedron

  • maxi – The number of points along the i direction. If maxi = 0 then the linear Serendipity element form of the specified shape is assumed.

  • maxj – The number of points along the j direction. If maxj = 0 then a Serendipity finite element is assumed. If 2 <= maxj <= 3 and 2 <= maxi <= 3, then a Lagrange finite element is assumed. If maxi != 0 then maxj = 0 or maxj = maxi.

  • maxk – The number of points along the k direction. Currently maxk must be equal to maxj.

void vfe_CorotSetMap(vfe_Corot *p, Vint nedofs, Vint loc[], Vint tag[])

set element degree of freedom map

Set the current number of element degrees of freedom and the degree of freedom map. The loc vector consists of nedofs one-based node indices into the element connectivity. A zero entry in loc denotes an internal elemental degree of freedom. The tag vector consists of nedofs defined constants which specify the degree of freedom type, eg. x-translation, y-rotation, etc. types. Degrees of freedom other than SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TZ, SYS_DOF_RX, SYS_DOF_RY, and SYS_DOF_RZ are ignored by the corotational transformations.

Errors

  • SYS_ERROR_VALUE is generated if an improper nedofs or loc is specified.

  • SYS_ERROR_ENUM is generated if an improper tag is specified.

Parameters:

  • p – Pointer to Corot object.

  • nedofs – Number of element degrees of freedom

  • loc – Vector of degree of freedom locations

  • tag – Vector of degree of freedom types

void vfe_CorotSetParami(vfe_Corot *p, Vint type, Vint iparam)

set corotational parameters

Set corotational parameters.

The addition of a geometric stiffness term to the element stiffness may be toggled using the VFE_GEOMSTIFF parameter. The geometric stiffness is added to the element stiffness in the function vfe_CorotPull() and may be interpreted as a “load” stiffness. The addition of the geometric stiffness due to load may be disabled by setting VFE_GEOMSTIFF to SYS_OFF. By default VFE_GEOMSTIFF is set to SYS_OFF.

For two-dimensional problems, you must set VFE_DIMENSION to 2. This will enforce a local corotational system whereby the local z-direction is always aligned with the global z-direction. By default VFE_DIMENSION is set to 3.

Errors

SYS_ERROR_ENUM is generated if an improper type or iparam is specified.

Parameters:

  • p – Pointer to Corot object.

  • type – Type of parameter to set 

    x=VFE_DIMENSION          Model dimension
 =VFE_GEOMSTIFF          Enable geometric stiffness

  • iparam – [out] Integer parameter value. 

    x=SYS_OFF                Disable
 =SYS_ON                 Enable
 =2 or 3                 Model dimension value

void vfe_CorotPull(vfe_Corot *p, Vdouble x[][3], Vdouble u[], Vdouble tm[3][3])

pull-back to reference coordinates

Transform global degree of freedoms into their equivalent infinitesimal values. As a side benefit, this function also returns the transformation matrix that relates the reference geometry to the corotational frame. This matrix can be used to transform stresses and strains from the reference to the current configurations. This function must to be called prior to vfe_CorotPush().

Errors

SYS_ERROR_COMPUTE is generated if the geometry defined by the coordinate array x is invalid - that is, nodes are coincident in space or three nodes do not form a proper triad.

Parameters:

  • p – Pointer to Corot object.

  • x – Node coordinates in reference frame.

  • u – Current degree of freedom vector

  • tm – Direction cosine matrix from reference to corotational frame

void vfe_CorotUpdateRot(vfe_Corot *p, Vdouble ur0[3], Vdouble dur[3], Vdouble ur[3])

update rotational degree of freedoms

Update a rotation vector ur0 represented as a pseudovector by an incremental rotation vector, dur, to produce a final updated rotation pseudovector ur.

Parameters:

  • p – Pointer to Corot object.

  • ur0 – Rotations at the end of the previous step

  • dur – Increment of rotations

  • ur – [out] Updated rotations

void vfe_CorotPush(vfe_Corot *p, Vint rflag, Vdouble r[], Vint kflag, Vdouble k[])

push-forward to current coordinates

Push the reaction force vector r and the element stiffness matrix k into the current configuration.

Errors

SYS_ERROR_OPERATION is generated if vfe_CorotPull() has not been called.

Parameters:

  • p – Pointer to Corot object.

  • rflag – Flag indicating whether reactions are to be pushed

  • r – [inout] in: Degree of freedom vector of reaction forces, out: Reaction forces are replaced by their pushed values.

  • kflag – Flag indicating whether stiffness matrix is to be pushed

  • k – [inout] in: Element stiffness matrix out: Element stiffness matrix is replaced by its pushed value.