"""Spring — 3D 2-node spring-damper (longitudinal mode only).
Default mode: longitudinal spring oriented along the I→J axis.
Real constants:
real[0] : K — spring stiffness (force / length)
real[1] : CV1 — linear damping (force · time / length, unused for K_e)
real[2] : CV2 — non-linear (cubic) damping (unused here)
real[3] : IL — initial length (unused; reference length is the nodal distance)
Stiffness::
K_e = K · [[ C, -C],
[-C, C]] with C_ij = d_i · d_j (3 × 3)
where ``(d_x, d_y, d_z)`` is the unit vector along I→J.
Spring is massless by convention (modal codes treat it as a
stiffness-only contribution); ``me`` returns a 6×6 zero matrix so
that assembly loops stay uniform.
Only the longitudinal mode is implemented; torsional and 2-D planar
variants need separate code paths and will be added alongside their
verification datasets.
References
----------
* Discrete spring / lumped connector element (3-D longitudinal
mode): Cook, Malkus, Plesha, Witt, *Concepts and Applications
of Finite Element Analysis*, 4th ed., Wiley, 2002, §2.3 and
§2.10 (direction-cosine rotation of a scalar spring between
two nodes in 3-D). Stiffness ``K · (d ⊗ d)`` on the I→J
axis is the standard result.
* Direction-cosine transformation for a bar / spring aligned
with an arbitrary 3-D axis: Zienkiewicz, O.C. and Taylor, R.L.,
*The Finite Element Method*, 7th ed., 2013, §2.4.4.
"""
from __future__ import annotations
from typing import ClassVar
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
def _axial_frame(coords: np.ndarray) -> tuple[float, np.ndarray]:
d = coords[1] - coords[0]
L = float(np.linalg.norm(d))
if L == 0.0:
# Coincident-node springs are physically legitimate — MAPDL's
# COMBIN14 / COMBIN40 routinely have ``E,n,n`` (same node twice)
# or two nodes at identical coordinates joined by a stiffness
# pinned to an external reference. Default the spring axis to
# global X (this matches MAPDL's KEYOPT(2)=1 default for
# coincident springs).
return 0.0, np.array([1.0, 0.0, 0.0], dtype=np.float64)
return L, d / L
[docs]
@register
class Spring(ElementBase):
name = "SPRING"
n_nodes = 2
n_dof_per_node = 3 # UX, UY, UZ — longitudinal mode only
vtk_cell_type = 3 # VTK_LINE
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
_, d = _axial_frame(np.asarray(coords, dtype=np.float64))
real = np.asarray(real, dtype=np.float64)
if real.size == 0:
raise ValueError("Spring requires REAL[0]=K (spring constant); got empty real set")
k_spring = float(real[0])
C = np.outer(d, d)
k = k_spring * np.block([[C, -C], [-C, C]])
return np.ascontiguousarray(k)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
return np.zeros((6, 6), dtype=np.float64)
[docs]
@register
class SpringTorsional(ElementBase):
"""3-D torsional spring (COMBIN14 with KEYOPT(3)=1).
Connects rotational DOFs (ROTX, ROTY, ROTZ) at the two end nodes
with a stiffness aligned along the I->J axis. Stiffness::
K_e = K_t * [[d⊗d, -d⊗d],
[-d⊗d, d⊗d]]
where ``d`` is the unit vector from node I to node J. Coincident
nodes default to a Z-aligned axis.
Real constants
--------------
real[0] : K_t - torsional spring constant (moment / radian)
Massless (mass / rotary inertia comes from companion MASS21 and
similar elements).
"""
name = "SPRING_TORSIONAL"
n_nodes = 2
n_dof_per_node = 3
vtk_cell_type = 3
dof_labels: ClassVar[tuple[str, ...]] = ("ROTX", "ROTY", "ROTZ")
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
_, d = _axial_frame(np.asarray(coords, dtype=np.float64))
real = np.asarray(real, dtype=np.float64)
if real.size == 0:
raise ValueError(
"SpringTorsional requires REAL[0]=K_t (torsional stiffness); got empty real set"
)
# Coincident-node torsional springs default to a Z-aligned axis
# rather than X (the longitudinal default) - matches MAPDL's
# COMBIN14 KEYOPT(3)=1 KEYOPT(2)=0 behaviour where the spring
# twist axis is the global Z when the two nodes overlap.
if np.allclose(coords[0], coords[1]):
d = np.array([0.0, 0.0, 1.0], dtype=np.float64)
k_spring = float(real[0])
C = np.outer(d, d)
k = k_spring * np.block([[C, -C], [-C, C]])
return np.ascontiguousarray(k)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
return np.zeros((6, 6), dtype=np.float64)