Module: shells

Shells and plates utility functions, e.g. for creation of plate / shell mesh.

  • Author: Johannes Gerstmayr

  • Date: 2026-01-11 (created)

  • Notes: For a list of plot colors useful for matplotlib, see also utilities.PlotLineCode(…) Extended ShellMesh class with symbolic slope computation via vertexMapping. Copied from ANCFThinPlatePrecurved.py; intended to replace exudyn.shells.ShellMesh once merged into the Exudyn library.

Function: _UnitVector

_UnitVector(v)

  • function description: Return the unit vector of v; raises ValueError if v is near-zero.

Function: AddNodeConstraints

AddNodeConstraints(mbs, nodeNumber, dofs, refNodeNumber = None)

  • function description: Constrain selected DOFs (0-8) of an ANCF plate node, to its reference value or to the same DOF of another node. Mirrors the GenerateStraightBeam pattern: one NodePointGround placed at the node’s reference position, one shared ground marker, one CoordinateConstraint per DOF. MarkerNodeCoordinate returns ODE2 displacement (without reference values), so default offset=0 pins each displacement to zero = node stays at its reference configuration.

  • input:

    • mbs: the MainSystem

    • nodeNumber: the NodePointSlope12 node

    • dofs: the coordinates of the node to constrain, 0-2 position, 3-5 slope x, 6-8 slope y

    • refNodeNumber: None - each selected coordinate keeps its reference value (ground); a node number - each selected coordinate equals the same coordinate of this node (C1 node-to-node connection)

  • output: list of CoordinateConstraint object indices.

Relevant Examples (Ex) and TestModels (TM) with weblink to github: ANCFThinPlateRevisionTest.py (TM)

Function: AddSphericalJointToGround

AddSphericalJointToGround(mbs, nodeNumber)

  • function description: Pin an ANCF node to ground: fix translations (DOF 0,1,2), slopes free.

Function: AddClampToGround

AddClampToGround(mbs, nodeNumber, fixedConstraints = None)

  • function description: Clamp an ANCF plate node to ground using a 9-element binary constraint vector. fixedConstraints: list/array of 9 values (0 or 1), one per DOF in the order [x, y, z, sx_x, sx_y, sx_z, sy_x, sy_y, sy_z]. A value of 1 fixes that DOF to its reference value; 0 leaves it free. Mirrors the fixedConstraintsNode0 convention of GenerateStraightLineANCFCable2D. Default (None): fix all 9 DOFs – equivalent to [1,1,1, 1,1,1, 1,1,1]. For a flat plate in the xy-plane (clamped left edge), the physically consistent choice that mirrors the 2D cable [1,1,0,1] convention is [1,1,1, 0,1,1, 1,0,1]: positions fixed, transverse slope components fixed (zero rotation), axial slope components (sx_x and sy_y, both = 1 in reference) left free to allow boundary stretching.

Function: _ThinPlateShapeFunctionDerivatives

_ThinPlateShapeFunctionDerivatives(xi, eta, scaleX, scaleY)

  • function description: First derivatives (d/dxi, d/deta) of the 12 ANCFThinPlate shape functions. Transcribed from CObjectANCFThinPlate::ComputeShapeFunctions_xy (C++), including the per-node slope scaling (slopesScalingX/Y parameters of the element). Node order: n0=(-1,-1), n1=(+1,-1), n2=(+1,+1), n3=(-1,+1); per node the SF order is [position, slopeX, slopeY]. Returns (sf_xi, sf_eta): two numpy arrays of length 12.

Function: _CreateMPCZeroMarker

_CreateMPCZeroMarker(mbs, nCoords = 36)

  • function description: Return a marker whose coordinate vector is clamped to zero. Used as marker0 in ObjectConnectorCoordinateVector MPCs so Exudyn does not warn about two identical markers (self-constraint via markerNumbers=[m, m]). One ground node can be shared by many element constraints in the same mbs.

Function: _AddSlopeDeviationConstraint

_AddSlopeDeviationConstraint(mbs, elem, rowFunc, groundMarker = None, nodalCoords = (-1.0, +1.0), components = (0, 1, 2))

  • function description: Add a linear MPC forcing the slope functional rowFunc(s) to equal the linear interpolation of its values at the edge ends, at the 2 interior Gauss points. rowFunc(s) must return the 12 shape-function-derivative coefficients of the constrained (transverse) slope at edge coordinate s in [-1, +1]. Returns the constraint object index.

Function: AddSlopeConformityConstraints

AddSlopeConformityConstraints(mbs, plateMesh, direction = 'x', components = (0, 1, 2))

  • function description: Suppress the non-conforming transverse-slope modes of ALL plate elements. The blended (ACM-type) ANCFThinPlate basis interpolates the transverse slope r_,x along edge lines NON-conformingly: between the nodes it deviates from the linear nodal interpolation by contributions of the far-node DOFs. This deviation is CONSTANT in xi (the eta-nonlinear terms of sf_,xi carry only constant-in-xi factors), so constraining it on ONE face per element removes it from the whole element; analogously for r_,y in eta. Effect: with direction=’x’, the membrane strain field of a beam-like strip (ny=1) reproduces ANCFCable2D EXACTLY (verified: identical L2/max membrane error vs fine cable reference), eliminating the spurious self-equilibrated strain oscillations across the width at clamps and inter-element lines. CAUTION: this stiffens the element (removes the non-conforming modes that contribute to 2D bending/twist softness). Intended for beam-like verification/benchmark setups, not for general shell meshes. Do NOT combine with AddDistributedClampToEdge on the same elements (duplicate equations -> singular Jacobian). LIMITATION: direction=’both’ creates REDUNDANT constraints (singular Jacobian) on meshes with ny>=2 (rank check: the y-direction rows of stacked element rows become dependent); use ‘both’ only for ny=1 (and ‘y’+nx>=2 analogously untested) – for the beam-like strip use the default ‘x’.

  • input:

    • mbs: the MainSystem

    • plateMesh: the ShellMesh whose elements are constrained

    • direction: ‘x’ (constrain r_,x along eta), ‘y’ (r_,y along xi), or ‘both’

    • components: slope vector components to constrain (default all)

  • output: list of constraint object indices.

Function: AddDistributedClampToEdge

AddDistributedClampToEdge(mbs, plateMesh, edgeKey, components = (0, 1, 2))

  • function description: Constrain the TRANSVERSE slope field along a clamped edge at the 2 interior Gauss points of every edge element (‘distributed clamp’). Background: nodal clamps (AddClampToGround) fix the edge POSITION field completely (it is conforming), but the transverse slope r_,n along the edge is the classic ACM-type NON-CONFORMING field: between the nodes it receives contributions from the interior/far-node DOFs of the edge element, so the clamp condition is violated mid-edge. This produces a spurious self-equilibrated strain oscillation across the clamped edge (boundary layer). Constraint formulation (per edge element, per Gauss point +-1/sqrt(3), per component): the transverse slope at the interior edge point must equal the LINEAR interpolation of the two nodal transverse slopes, r_,n(s*) - [(1-s*)/2 r_,n(-1) + (1+s*)/2 r_,n(+1)] = 0. This kills exactly the non-conforming deviation while leaving the nodal slope DOFs (e.g. free axial stretch sx_x of a clamp like [1,1,1, 0,1,1, 1,0,1]) untouched – mirroring an ideal 1D (cable-like) clamp. Since the transverse slope is cubic along the edge, 2 interior points reduce it exactly to the linear nodal interpolation. Each constraint is a LINEAR multipoint constraint on the element’s 36 ODE2 coordinates (MarkerObjectODE2Coordinates + ObjectConnectorCoordinateVector with exact constant Jacobian). NOTE: this makes the strain AT the clamp face exact (matches an ideal 1D clamp), but the suppressed non-conforming mode reappears at the next inter-element line (overall L2 error barely changes). To suppress the non-conformity in the WHOLE mesh (e.g. to reproduce ANCFCable2D exactly with a beam-like strip), use AddSlopeConformityConstraints instead – but do not combine both on the same elements (duplicate equations -> singular Jacobian).

  • input:

    • mbs: the MainSystem

    • plateMesh: the ShellMesh of the edge

    • edgeKey: ‘left’/’right’/’bottom’/’top’ – should also be nodally clamped

    • components: vector components (0=x, 1=y, 2=z) of the transverse slope to constrain; default all three.

  • output: list of constraint object indices.

Function: ApplyEdgeLoad

ApplyEdgeLoad(mbs, plateMesh, edgeKey, loadVector = None, torqueVector = None, distributed = False)

  • function description: Apply a uniformly distributed force/torque along one plate edge. Accepts the TOTAL load/torque; the force is distributed in one of two ways: distributed=False (default, nodal): Trapezoidal rule on the edge nodes: corner nodes receive 1x the per-element base, interior nodes receive 2x, so the sum equals the total exactly. NOTE: nodal point forces are NOT work-equivalent to a uniform line load on the cubic (Hermite) edge interpolation – the slope-conjugate load components are missing, which excites a spurious Saint-Venant boundary layer (local membrane/bending oscillation across the width) at the loaded edge. distributed=True (consistent line load): Work-equivalent uniform line load: per edge element, the force is applied at the 2 Gauss points of the edge (local coordinate +-1/sqrt(3)) via MarkerBodyPosition with weight 1/2 each. Since the edge field is cubic and the load direction is constant, 2-point Gauss integration of f = int S^T(edge) q ds is EXACT – including the Hermite slope-conjugate components that nodal forces miss. Use this to reproduce ideal beam-like load introduction (e.g. matching an ANCFCable2D tip load on a beam-like strip). Torques are always applied nodally (trapezoidal rule).

  • input:

    • mbs: the MainSystem

    • plateMesh: the ShellMesh of the edge

    • edgeKey: the edge, ‘left’, ‘right’, ‘bottom’ or ‘top’ of plateMesh.boundaryNodeNumbers

    • loadVector: 3-component total force [N] or None

    • torqueVector: 3-component total torque [N·m] or None

    • distributed: if True, apply force as consistent line load (see above)

Function: AddRotationalSpringDamper

AddRotationalSpringDamper(mbs, nodeA, nodeB, dofs, stiffness, damping = 0.0)

  • function description: Add a CoordinateSpringDamper on each selected DOF between nodeA and nodeB. Used to add a rotational spring/damper at an ANCF plate hinge: dofs=[3,4,5] resists relative rotation about the y-edge-axis (sx slopes). dofs=[6,7,8] resists relative twist about the x-edge-axis (sy slopes).

  • input:

    • mbs: the MainSystem

    • nodeA: the first node

    • nodeB: the second node

    • dofs: the node coordinates to couple

    • stiffness: the stiffness of each CoordinateSpringDamper

    • damping: the damping of each CoordinateSpringDamper

  • output: list of CoordinateSpringDamper object indices.

Function: AddEdgeSpringDamper

AddEdgeSpringDamper(mbs, nodesA, nodesB, dofs, totalStiffness, totalDamping = 0.0)

  • function description: Distribute a total rotational spring-damper along matching edge node lists. Mirrors ApplyEdgeLoad: distributes totalStiffness via the trapezoidal rule so that the effective stiffness per unit length is uniform regardless of mesh refinement. Corner nodes (first/last) receive 1x the base weight; interior nodes receive 2x, giving the correct trapezoidal integration.

  • input:

    • mbs: the MainSystem

    • nodesA: list of node numbers on side A of the hinge (e.g. plate1 right edge)

    • nodesB: list of node numbers on side B of the hinge (e.g. plate2 left edge)

    • dofs: DOF indices to couple, e.g. [3,4,5] for x-slope (bending hinge)

    • totalStiffness: total rotational stiffness [Nm/rad] of the whole hinge

    • totalDamping: total rotational damping [Nm·s/rad] of the whole hinge

  • output: list of all CoordinateSpringDamper object indices.

CLASS ShellMesh (in module shells)

class description: class for generation, representation of plate and shell meshes; creaton of Exudyn elements

Class function: init

__init__(self, vertices = [[-1,-1,0],[ 1,-1,0],[ 1, 1,0],[-1, 1,0]], numberOfElementsX = 1, numberOfElementsY = 1, youngsModulus = None, poissonsRatio = None, density = None, thickness = None, massProportionalDamping = 0., thicknessAtNodes = None, thicknessFunction = None, stiffnessProportionalDamping = 0., surfaceMap = None)

  • class function description: initialize rectangular shell mesh with geometry, discretization and physics parameters

  • input:

    • vertices: list of four 3D vectors (numpy array or list), sorted [bottom-left, bottom-right, top-right, top-left]; defining the reference positions of the corner nodes; if further transformations are added, use unit coordinates!

    • numberOfElementsX: number of elements in x-direction

    • numberOfElementsY: number of elements in y-direction

    • youngsModulus: Young’s modulus; used for calculation of membrane and bending stiffness

    • poissonsRatio: Poisson’s ratio for inplane shear deformation

    • density: average density of plate/shell

    • thickness: thickness of plate/shell

    • massProportionalDamping: damping parameter which introduces damping proportional to distributed mass

    • thicknessAtNodes: optional thickness at each node, in the order of the nodes; interpolated bilinearly in each element

    • thicknessFunction: optional function f(x, y) of the thickness at the global x and y of a point, which must accept exu.symbolic.Real; its value and its gradients along the slopes of each node give 12 thickness values per element

    • stiffnessProportionalDamping: Kelvin-Voigt damping coefficient [s] of the membrane and the bending stiffness

    • surfaceMap: optional function f(x, y) -> [X, Y, Z] that maps the positions (x, y) of the nodes in the rectangle of the vertices onto a surface, e.g. SurfaceMap.Cylinder(…); the slopes of the nodes follow from its derivatives

  • notes: x-axis is aligned with bottom (y=min) and top (y=max); y-axis is aligned with left (x=min) and right (x=max)

Relevant Examples (Ex) and TestModels (TM) with weblink to github: ANCFThinPlateRevisionTest.py (TM), ANCFThinPlateTests.py (TM), energiesFlexibleBodiesTest.py (TM)

CLASS SurfaceMap (in module shells)

class description: Maps of a rectangle onto a surface, for ShellMesh(surfaceMap=…): each static method returns a function f(x, y) -> [X, Y, Z] of the position (x, y) in the rectangle of the vertices of the mesh. A map written by the user has the same form. It is called with exu.symbolic.Real values for x and y, so that the slopes of the nodes follow exactly from its derivatives: use exu.symbolic.sin, exu.symbolic.cos and the other functions of exu.symbolic, which take floats as well; a map that computes with numpy (BezierStrip, ConeFrustum, Cylinder) gets its slopes by central differences.

Class function: SkewParallelogram

SkewParallelogram(skewY = 0.0, skewX = 0.0)

  • class function description: Skew a rectangle into a parallelogram. skewY: adds x += skewYy (shear x with y). skewX: adds y += skewXx (shear y with x).

Class function: Trapezoid

Trapezoid(Lx, Ly, topScale = 1.0, aboutMidline = True)

  • class function description: Trapezoid by scaling x-span linearly with y. topScale<1: narrower at top; topScale>1: wider at top. aboutMidline=True keeps the plate centreline fixed.

Class function: CurvedEdge

CurvedEdge(Lx, Ly, amp = 0.0, mode = 'top', shape = 'sin', aboutMidline = True)

  • class function description: Curve one edge by an x-shift that varies with y (and x). mode: ‘top’ or ‘bottom’ — where the curvature reaches full amplitude. shape: ‘sin’ (sinusoidal in x) or ‘parabola’. amp: maximum x-shift [m].

Class function: OutOfPlaneWarp

OutOfPlaneWarp(Lx, Ly, amp = 0.0)

  • class function description: Initial out-of-plane warp vanishing at all four boundaries: z = amp * sin(pix/Lx) * sin(piy/Ly).

Class function: BezierStrip

BezierStrip(P0, P1, P2, P3, Lx, Ly, up = (0, 0, 1), useArcLength = True, nArc = 400)

  • class function description: Map a plate onto a cubic Bezier centerline strip. x in [0, Lx] follows the curve (optionally by arc-length). y in [0, Ly] offsets across the strip width.

Class function: ConeFrustum

ConeFrustum(Lx, Ly, r0, r1, phiDeg = 20.0, thetaCenterDeg = 0.0, centered = True, yMode = "arclength")

  • class function description: Map a rectangle onto a conical frustum mantle. x in [0,Lx] along the axis; y in [0,Ly] around the circumference. r0, r1: radii at x=0 and x=Lx. phiDeg: angular extent [deg].

Class function: Cylinder

Cylinder(Lx, Ly, R, phiDeg = 90.0, thetaCenterDeg = 0.0, centered = False, yMode = 'arclength')

  • class function description: Map a rectangle onto a cylindrical mantle (special case of ConeFrustum with r0=r1=R).

Class function: HemisphericalShell

HemisphericalShell(Lx, Ly, R, alpha, center = (0.0,0.0,0.0), thetaOffset = 0.0, thetaRange = None)

  • class function description: Map a rectangle onto a hemispherical shell with a circular apex cutout. R: radius. alpha: cutout angle from +z axis [rad]. x maps to azimuthal angle; y maps from cutout edge (y=0) to equator (y=Ly).

Class function: ToroidalPanel

ToroidalPanel(Lx, Ly, r1 = 1.5, r2 = 0.5)

  • class function description: Map a rectangle onto a toroidal panel. x -> meridional angle q1 in [-pi/2, pi/2]; y -> azimuthal angle q2 in [0, pi/2].