Source code for femorph_solver.elements.spring

"""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

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:
        raise ValueError("SPRING: coincident nodes, element length is zero")
    return L, d / L




@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



    @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)




    @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)