Rigid bodies

A rigid body couples the motion of many surface (visualization) nodes to a single reference point (RP). The surface nodes are not independent degrees of freedom: they carry no field variables and are never written to a node field. Instead, their current configuration is a pure function of the RP’s kinematics, so the only unknowns a rigid body contributes to the global system are the RP displacement and (in 3D) rotation DOFs. Rigid bodies are consumed by contact constraints such as surfacetodiscreterigidbodypenalty and are rendered by the output managers via their surface nodes.

Kinematics

The reference point owns a total displacement \(\mathbf{u}_{RP}\) and, in 3D, a rotation pseudo-vector \(\boldsymbol{\theta}_{RP}\) (an axis-angle vector: direction = rotation axis, magnitude = rotation angle in radians). The orientation is recovered through the SO(3) exponential map (Rodrigues’ formula)

\[\mathbf{R}(\boldsymbol{\theta}) = \exp(\mathrm{skew}(\boldsymbol{\theta})) = \mathbf{I} + \frac{\sin\varphi}{\varphi}\,\mathrm{skew}(\boldsymbol{\theta}) + \frac{1-\cos\varphi}{\varphi^2}\,\mathrm{skew}(\boldsymbol{\theta})^2 , \qquad \varphi = \lVert\boldsymbol{\theta}\rVert ,\]

where \(\mathrm{skew}(\mathbf{v})\) is the skew-symmetric matrix with \(\mathrm{skew}(\mathbf{v})\,\mathbf{x} = \mathbf{v}\times\mathbf{x}\). The stored DOF is the raw total pseudo-vector; because \(\exp(\boldsymbol{\theta}+\delta\boldsymbol{\theta}) \neq \exp(\boldsymbol{\theta})\exp(\delta\boldsymbol{\theta})\) in general, linearizing any quantity with respect to \(\boldsymbol{\theta}_{RP}\) requires the right Jacobian \(\mathbf{J}_r(\boldsymbol{\theta}_{RP})\) of the exponential map (this is what makes the discrete-rigid-body contact tangent geometrically exact — see Constraints).

A surface node with initial position offset \(\boldsymbol{\rho}_i = \mathbf{X}_i - \mathbf{X}_{RP}\) relative to the RP moves rigidly to the current configuration

\[\mathbf{x}_i = (\mathbf{X}_{RP} + \mathbf{u}_{RP}) + \mathbf{R}\,\boldsymbol{\rho}_i .\]

updateKinematics() writes these current coordinates onto the surface nodes after every converged increment, which is what the output managers read to draw the moving body. Because the surface nodes are not DOFs, this is a display / broadphase convenience, not part of the solve: every nonlinear solver invokes it through updateRigidBodies(), so all solver variants (serial, parallel, arc-length, explicit) keep the visualization geometry consistent. Contact constraints do not rely on it — they query the surface with the RP pose taken directly from the current Newton iterate (see below). That pose is provided by poseFromDofs(). The integrated contact constraint searches the body-frame triangles of referenceTriangles(), which requires a closed surface without zero-area triangles, so that every triangle has a defined outward normal.

Mass properties

For dynamics (and for reporting), the exact volume, mass, center of mass and inertia tensor of a homogeneous solid bounded by the (closed, manifold) surface mesh are computed by the polyhedronmassproperties module. Each surface facet is fan-triangulated and, together with the origin, forms a signed tetrahedron; the divergence theorem then turns the volume integrals of \(1\), \(\mathbf{x}\), and \(\mathbf{x}\otimes\mathbf{x}\) into a sum over those tetrahedra, e.g. for the volume

\[V = \int_\Omega \mathrm{d}V = \frac{1}{6}\sum_{\text{facets}} \mathbf{a}\cdot(\mathbf{b}\times\mathbf{c}) ,\]

with \(\mathbf{a},\mathbf{b},\mathbf{c}\) the facet’s triangle vertices. Consistent outward face winding makes the signed tetrahedra cancel correctly for non-convex bodies, so the result is exact (no volumetric mesh required). The inertia tensor is returned about the center of mass, in the global (non-rotated) axes of the input vertices.

Surface query

Contact needs, for an arbitrary array of world-frame query points, the signed distance to the rigid surface (negative inside) and the outward normal of the closest facet, evaluated for the body’s current pose. This is provided by DiscreteSurfaceQuery, which wraps VTK’s vtkImplicitPolyDataDistance and a static cell locator built once on the reference mesh. Rather than re-transform the whole mesh each query, the query points are mapped into the body’s local (reference) frame,

\[\mathbf{p}_{\text{local}} = \mathbf{X}_{RP} + \mathbf{R}^{\top}\big(\mathbf{p} - (\mathbf{X}_{RP}+\mathbf{u}_{RP})\big) ,\]

the distance/closest-facet lookup is done there, and the returned normals are rotated back to the world frame by \(\mathbf{R}\).

Note

The implementation stores points as row vectors (shape (N, 3)), for which v.dot(M) computes \(\mathbf{M}^{\top}\mathbf{v}\) per row. Consequently the inverse map above is coded as localCoords.dot(rotation_matrix) (i.e. \(\mathbf{R}^{\top}\)) and the forward map of normals as normals.dot(rotation_matrix.T) (i.e. \(\mathbf{R}\)) — the opposite of what a naive read of the transpose suggests. This convention is documented in the source and was verified numerically against a known rotated mesh; do not “fix” the transpose without re-deriving it.

Outward normals are established from the mesh topology (compute_normals(..., auto_orient_normals=True)) rather than trusting the source file’s face winding, so the sign of the signed distance and the direction of the returned normals are reliable for any correctly closed input mesh.

Creating a discrete rigid body

Discrete rigid bodies are normally created from a surface mesh file (Exodus / STL / OBJ) with the discreterigidbodygenerator model generator, which loads the mesh, creates the RP and surface nodes, computes the mass properties, and registers the body in model.rigidBodies. See Generators for the generator options and an input-file example.

API

Rigid body base class

Module edelweissfe.rigidbodies.rigidbody

class edelweissfe.rigidbodies.rigidbody.RigidBody(name, model)

Abstract Base Class for all rigid bodies in the explicit and implicit solver framework.

Methods

getCurrentKinematics()	Retrieve the current kinematic state of the rigid body.
getVisualizationElements()	Returns the geometric elements for visualization (e.g., facets).
getVisualizationField(fieldName)	Returns the mapped results for visualization (e.g., displacement) on the nodes.
getVisualizationNodes()	Returns the nodes that should be visualized for this rigid body.
referenceCoordinatesOfMovedNodes()	The reference coordinates of the nodes whose coordinates this body moves along with it.
updateKinematics([timeStep])	Update the kinematics of the rigid body according to its prescribed or computed motion.

Parameters:

name (str)

abstractmethod getCurrentKinematics()

Retrieve the current kinematic state of the rigid body.

abstractmethod getVisualizationElements()

Returns the geometric elements for visualization (e.g., facets). Each element is a dict with ‘type’ and ‘nodes’ keys.

Return type:

List[Dict[str, Any]]

abstractmethod getVisualizationField(fieldName)

Returns the mapped results for visualization (e.g., displacement) on the nodes.

Parameters:

fieldName (str)

Return type:

ndarray

abstractmethod getVisualizationNodes()

Returns the nodes that should be visualized for this rigid body.

Return type:

List

referenceCoordinatesOfMovedNodes()

The reference coordinates of the nodes whose coordinates this body moves along with it.

A node’s coordinates are normally its reference coordinates, and the topology fingerprint (fingerprint()) hashes them. A rigid body that writes its current configuration into its nodes reports their reference coordinates here, so that the fingerprint describes the mesh and not the body’s current position.

Returns:

The reference coordinates, keyed by node label. Empty if this body moves no nodes.

Return type:

Dict[int, numpy.ndarray]

abstractmethod updateKinematics(timeStep=None)

Update the kinematics of the rigid body according to its prescribed or computed motion.

Discrete rigid body

Module edelweissfe.rigidbodies.discreterigidbody

class edelweissfe.rigidbodies.discreterigidbody.DiscreteRigidBody(name, model, surfaceMesh, *args, **kwargs)

Discrete Rigid Body entity.

It represents a rigid body defined by a discrete surface mesh. Instances are normally created via generateDiscreteRigidBodyFromMeshFile(), which handles loading the mesh and creating the surface/reference-point nodes – this class itself only deals with the rigid body’s kinematics.

Methods

getAABB([kinematics])	Axis-aligned bounding box of the surface nodes' current positions.
getCurrentKinematics()	Return the current rigid body motion.
getVisualizationElements()	Returns the geometric elements for visualization (e.g., facets).
getVisualizationField(fieldName)	Compute a field's values on the surface (visualization) nodes.
getVisualizationNodes()	Returns the nodes that should be visualized for this rigid body.
poseFromDofs(rpDisplacement, rpRotation)	The rigid body pose for given reference point DOF values.
querySurface(coords[, proximityDistance, ...])	Compute signed distances and outward face normals of the rigid surface, in its current configuration, for an array of query points.
referenceCoordinatesOfMovedNodes()	The reference coordinates of the surface nodes, which updateKinematics() moves.
referenceTriangles()	The triangles of the rigid surface and their outward unit normals, in the reference configuration (the body frame).
updateKinematics([timeStep])	Update the kinematics of the rigid body according to its prescribed or computed motion.

getAABB(kinematics=None)

Axis-aligned bounding box of the surface nodes’ current positions.

Parameters:

kinematics (tuple, optional) – See _currentAndReferenceSurfaceCoordinates(). If not given, uses the stacked array cached by the last updateKinematics() call (or, before the first such call, rebuilds it from the surface nodes’ own coordinates), as a fast path for the common case.

getCurrentKinematics()

Return the current rigid body motion.

Returns:

The total RP displacement, the 3x3 rotation matrix, and the reference (initial) RP coordinates. The current RP position is the sum of the first and third entries.

Return type:

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

getVisualizationElements()

Returns the geometric elements for visualization (e.g., facets). Each element is a dict with ‘type’ and ‘nodes’ keys.

Return type:

List[Dict[str, Any]]

getVisualizationField(fieldName)

Compute a field’s values on the surface (visualization) nodes.

The surface nodes carry no FieldVariables of their own – they are fully determined by the reference point’s already-solved kinematics – so this is computed directly rather than looked up from a NodeField.

Parameters:

fieldName (str) – The name of the field. Currently only "displacement" is supported.

Raises:

ValueError – If fieldName is not "displacement".

Return type:

ndarray

getVisualizationNodes()

Returns the nodes that should be visualized for this rigid body.

Return type:

List

poseFromDofs(rpDisplacement, rpRotation)

The rigid body pose for given reference point DOF values.

A contact constraint evaluates the pose from its own current trial solution, not from the last converged NodeField state that getCurrentKinematics() reads.

Parameters:

• rpDisplacement (ndarray) – The displacement of the reference point.

• rpRotation (ndarray) – The rotation pseudo-vector of the reference point.

Returns:

The current reference point position, the rotation matrix \(R\), and dPhysicalSpin_dTheta \(= R \, J_r(\boldsymbol{\theta})\), which maps a perturbation of the stored pseudo-vector DOF onto the spatial rotation it produces (the identity only at zero rotation), see rightJacobianSO3().

Return type:

tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]

querySurface(coords, proximityDistance=None, kinematics=None)

Compute signed distances and outward face normals of the rigid surface, in its current configuration, for an array of query points.

Parameters:

• coords (numpy.ndarray, shape (nPoints, 3)) – The query coordinates.

• proximityDistance (float, optional) – If given, a broadphase AABB check (the current AABB inflated by this distance) is performed first; points outside are assigned a distance of inf and a zero normal without querying the surface.

• kinematics (tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray], optional) – An explicit (u_rp, R, rp_initial) triple to use for the rigid body’s pose instead of getCurrentKinematics()’s (last-converged-increment) NodeField state. A caller that already holds a fresher/current RP kinematic state – e.g. a contact constraint mid-Newton-iteration, which has direct access to the trial solution vector – should pass it here so the query is evaluated for a pose consistent with the rest of that caller’s own current state (both this narrow-phase query and the broadphase AABB check use it consistently).

Returns:

• dists (numpy.ndarray, shape (nPoints,)) – The signed distances (negative = penetration).

• normals (numpy.ndarray, shape (nPoints, 3)) – The outward unit normals of the closest faces.

referenceCoordinatesOfMovedNodes()

The reference coordinates of the surface nodes, which updateKinematics() moves.

Returns:

The reference coordinates, keyed by surface node label.

Return type:

dict

referenceTriangles()

The triangles of the rigid surface and their outward unit normals, in the reference configuration (the body frame).

The rigid surface must be closed, so that outside and inside are well defined, and must not contain triangles of zero area, whose normal is undefined. Both are checked here, since a violation otherwise shows up only as unphysical contact forces much later in a run.

Returns:

The triangle vertex coordinates, of shape (nTriangles, 3, 3), and the outward unit normals, of shape (nTriangles, 3).

Return type:

tuple[numpy.ndarray, numpy.ndarray]

updateKinematics(timeStep=None)

Update the kinematics of the rigid body according to its prescribed or computed motion.

Polyhedron mass properties

Module edelweissfe.utils.polyhedronmassproperties

Exact mass properties (volume, mass, center of mass, and inertia tensor) of an arbitrary closed polyhedral surface mesh, computed via divergence-theorem integration over signed tetrahedra.

This module has no dependency on any particular mesh library (PyVista, VTK, …) or on EdelweissFE’s model/entity classes: it operates purely on raw vertex coordinates and face connectivity, so it can be reused wherever exact mass properties of a solid described by its boundary surface are needed (e.g. rigid bodies, indenters, generators reading mesh files).

class edelweissfe.utils.polyhedronmassproperties.PolyhedronMassProperties(volume, mass, centerOfMass, inertia)

Create new instance of PolyhedronMassProperties(volume, mass, centerOfMass, inertia)

Methods

count(value, /)	Return number of occurrences of value.
index(value[, start, stop])	Return first index of value.

Parameters:

• volume (float)

• mass (float)

• centerOfMass (ndarray)

• inertia (ndarray)

centerOfMass: ndarray

Alias for field number 2

inertia: ndarray

Alias for field number 3

mass: float

Alias for field number 1

volume: float

Alias for field number 0

edelweissfe.utils.polyhedronmassproperties.computePolyhedronMassProperties(vertices, faces, density=1.0)

Compute the exact mass properties of a closed polyhedral surface mesh.

Computes the volume, mass, center of mass, and inertia tensor of the homogeneous solid enclosed by a closed (watertight), consistently oriented polyhedral surface, via divergence-theorem integration over signed tetrahedra formed between the coordinate origin and each face [1]. The result is exact (up to floating-point round-off) for any polyhedron describable by its boundary faces – it is not limited to canonical primitives such as boxes, cylinders, or spheres.

Parameters:

• vertices (numpy.ndarray, shape (nVertices, 3)) – The Cartesian coordinates of all vertices referenced by faces.

• faces (sequence of sequence of int) – The faces of the closed surface, each given as the (0-based) indices of its vertices in vertices, in consistent winding order (all faces oriented outward, or all oriented inward – a single globally reversed mesh is detected and corrected automatically). Faces may be triangles, quadrilaterals, or other planar convex polygons.

• density (float, optional) – The (uniform) mass density of the solid. Default is 1.0, in which case mass numerically equals volume.

Returns:

A named tuple with fields volume, mass, centerOfMass, and inertia (the latter about the center of mass).

Return type:

PolyhedronMassProperties

Raises:

ValueError – If the computed volume is (numerically) zero, which indicates that faces does not describe a closed, watertight surface.

Notes

The surface must be closed (watertight) and consistently wound for the result to be meaningful; this is not verified here, as doing so requires a topological check of the mesh (e.g. that every edge is shared by exactly two oppositely-wound faces), which is outside the scope of this purely geometric utility.

References

[1]

B. Mirtich, “Fast and accurate computation of polyhedral mass properties,” Journal of Graphics Tools, 1(2), 1996.

Examples

>>> import numpy as np
>>> vertices = np.array([[0,0,0],[1,0,0],[1,1,0],[0,1,0],
...                       [0,0,1],[1,0,1],[1,1,1],[0,1,1]], dtype=float)
>>> faces = [[0,3,2,1],[4,5,6,7],[0,1,5,4],[1,2,6,5],[2,3,7,6],[3,0,4,7]]
>>> props = computePolyhedronMassProperties(vertices, faces)
>>> round(props.volume, 6)
1.0
>>> np.allclose(props.centerOfMass, [0.5, 0.5, 0.5])
True

Discrete surface query

Module edelweissfe.utils.discretesurfacequery

class edelweissfe.utils.discretesurfacequery.DiscreteSurfaceQuery(mesh)

Initializes the query engine for a given PyVista mesh.

Uses static VTK locators and evaluators to avoid memory leaks during repeated distance evaluations.

Parameters:

mesh (pyvista.PolyData) – The PyVista mesh to query.

Methods

query(coords[, translation, ...])

Batched query to compute signed distance and closest face normals for an array of points.

query(coords, translation=None, rotation_matrix=None, rotation_center=None)

Batched query to compute signed distance and closest face normals for an array of points.

The coordinate transform and result assembly are vectorized, but the actual VTK distance/closest-cell evaluation loops over points one at a time – this is a potential performance bottleneck for very large point counts.

Parameters:

• coords (numpy.ndarray) – Array of shape (N, 3) containing the query coordinates.

• translation (numpy.ndarray, optional) – A translation vector of the rigid body RP.

• rotation_matrix (numpy.ndarray, optional) – A 3x3 rotation matrix of the rigid body.

• rotation_center (numpy.ndarray, optional) – The center of rotation (initial position of the RP).

Returns:

• dists (numpy.ndarray) – Array of shape (N,) containing the signed distance (negative = inside/penetration).

• normals (numpy.ndarray) – Array of shape (N, 3) containing the outward normal vectors on the closest faces.

• .. note:: – Note on the rotation convention: both localCoords/normals here and e.g. getCurrentKinematics() use the same forward convention: a body-frame vector v is mapped to the world frame via v.dot(R.T). Because these are row vectors (shape (N, 3)), v.dot(M) computes M.T @ v per row, not M @ v – so v.dot(R.T) is R @ v (forward: body -> world), and v.dot(R) is R.T @ v (inverse: world -> body). This is the opposite of what a naive read of .dot(rotation_matrix) vs. .dot(rotation_matrix.T) suggests – verified numerically against a known rotated mesh before relying on it, since it is easy to get backwards (do not “fix” the transpose here without re-deriving/re-checking this convention first).

Return type:

tuple[ndarray, ndarray]