"""PointMass — concentrated point mass.
Two single-node element variants:
* :class:`PointMass` — translational-only (3 DOFs/node ``UX, UY, UZ``).
* :class:`PointMassRotary` — translational + rotary inertia (6 DOFs/node
``UX, UY, UZ, ROTX, ROTY, ROTZ``), used at nodes that already carry
rotational DOFs (e.g. attached to a beam or shell).
The MAPDL ``MASS21`` 3-D KEYOPT(3) variants map onto these as:
* ``KEYOPT(3) ∈ {0, unspecified}`` — 3-D mass *with* rotary inertia → :class:`PointMassRotary`.
* ``KEYOPT(3) = 2`` — 3-D mass *without* rotary inertia → :class:`PointMass`.
* ``KEYOPT(3) = 4`` — 2-D mass *without* rotary inertia → :class:`PointMass`.
* ``KEYOPT(3) = 3`` — 2-D mass *with* rotary inertia (``UX, UY, ROTZ``)
→ currently routed to :class:`PointMassRotary`; the unused
out-of-plane DOFs receive zero mass and are pinned by the deck's
symmetry constraints in 2-D usage.
Real constants
--------------
:class:`PointMass`::
real[0]: MASSX — translational mass along x (mandatory)
real[1]: MASSY — translational mass along y (defaults to MASSX)
real[2]: MASSZ — translational mass along z (defaults to MASSX)
:class:`PointMassRotary` (matches MAPDL ``MASS21`` real-constant order)::
real[0]: MASSX — translational mass along x (defaults to 0.0)
real[1]: MASSY — translational mass along y (defaults to MASSX)
real[2]: MASSZ — translational mass along z (defaults to MASSX)
real[3]: IXX — rotary inertia about x (defaults to 0.0)
real[4]: IYY — rotary inertia about y (defaults to IXX)
real[5]: IZZ — rotary inertia about z (defaults to IXX)
A ``MASS21`` deck that omits leading slots leaves zeros — e.g.
``R, 2, , , , 31E-3`` is ``(MASSX=0, MASSY=0, MASSZ=0, IXX=31E-3,
IYY=0, IZZ=0)``: pure rotary inertia about the X axis with no
translational mass.
Matrices
--------
:class:`PointMass`::
K_e = 0₃×₃ (a point mass has no stiffness contribution)
M_e = diag(MASSX, MASSY, MASSZ)
:class:`PointMassRotary`::
K_e = 0₆×₆
M_e = diag(MASSX, MASSY, MASSZ, IXX, IYY, IZZ)
Lumped and consistent formulations coincide for a one-node element.
References
----------
* Concentrated / point-mass element (lumped nodal mass, no stiffness
contribution): Cook, Malkus, Plesha, Witt, *Concepts and
Applications of Finite Element Analysis*, 4th ed., Wiley, 2002,
§16.3 (discussion of lumped nodal masses). Bathe, K.J., *Finite
Element Procedures*, 2nd ed., Prentice Hall, 2014, §4.2.4.
* Anisotropic diagonal mass tensor (``diag(MASSX, MASSY, MASSZ)``,
supports different masses on each translational axis): standard
generalisation of the isotropic ``m·I₃`` lumped mass.
* Rotary-inertia tensor as a diagonal ``diag(IXX, IYY, IZZ)``: Meirovitch,
L., *Methods of Analytical Dynamics*, McGraw-Hill, 1970, §4.7
(rigid-body kinetic energy decomposition into translational and
rotational parts on the principal axes).
"""
from __future__ import annotations
import numpy as np
from femorph_solver.elements._base import ElementBase
from femorph_solver.elements._registry import register
def _mass_vector(real: np.ndarray) -> np.ndarray:
r = np.asarray(real, dtype=np.float64)
if r.size == 0:
raise ValueError("PointMass requires REAL[0]=MASSX (point mass); got empty real set")
mx = float(r[0])
my = float(r[1]) if r.size > 1 else mx
mz = float(r[2]) if r.size > 2 else mx
return np.array([mx, my, mz], dtype=np.float64)
def _rotary_diagonal(real: np.ndarray) -> np.ndarray:
"""Unpack a ``MASS21``-style 6-slot real-constant set.
Slots default to zero so that decks like ``R, 2, , , , 31E-3``
(rotary inertia about X only, no translational mass) round-trip
without raising. Trailing rotary slots fall back to ``IXX`` to
match MAPDL's "unset → MASSX/IXX" defaulting on the leading
translational triple.
"""
r = np.asarray(real, dtype=np.float64)
n = r.size
mx = float(r[0]) if n > 0 else 0.0
my = float(r[1]) if n > 1 else mx
mz = float(r[2]) if n > 2 else mx
ixx = float(r[3]) if n > 3 else 0.0
iyy = float(r[4]) if n > 4 else ixx
izz = float(r[5]) if n > 5 else ixx
return np.array([mx, my, mz, ixx, iyy, izz], dtype=np.float64)
[docs]
@register
class PointMass(ElementBase):
name = "POINT_MASS"
n_nodes = 1
n_dof_per_node = 3 # UX, UY, UZ — translational-only (no rotary inertia)
vtk_cell_type = 1 # VTK_VERTEX
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
return np.zeros((3, 3), dtype=np.float64)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
return np.diag(_mass_vector(real))
[docs]
@register
class PointMassRotary(ElementBase):
name = "POINT_MASS_ROTARY"
n_nodes = 1
n_dof_per_node = 6 # UX, UY, UZ, ROTX, ROTY, ROTZ
vtk_cell_type = 1 # VTK_VERTEX
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
return np.zeros((6, 6), dtype=np.float64)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
return np.diag(_rotary_diagonal(real))
[docs]
@register
class PointMass2DRotary(ElementBase):
"""MAPDL ``MASS21`` KEYOPT(3) = 3 - 2-D mass with rotary inertia.
The 2-D ``MASS21`` real-constant layout differs from the 3-D form::
real[0] : MASS (translational, applied to both UX and UY)
real[1] : IZZ (rotary inertia about Z)
so the standard 6-slot ``_rotary_diagonal`` unpack misreads
``R, 2, MASS, IZZ`` as ``MASSX = MASS, MASSY = IZZ`` and dumps
the IZZ value into the wrong DOF. vm47 (Thomson torsional
pendulum) needs the dedicated layout to match the textbook
``f = sqrt(K_t / IZZ) / (2 pi)``.
"""
name = "POINT_MASS_2D_ROTARY"
n_nodes = 1
n_dof_per_node = 6
vtk_cell_type = 1
[docs]
@staticmethod
def ke(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
) -> np.ndarray:
return np.zeros((6, 6), dtype=np.float64)
[docs]
@staticmethod
def me(
coords: np.ndarray,
material: dict[str, float],
real: np.ndarray,
lumped: bool = False,
) -> np.ndarray:
r = np.asarray(real, dtype=np.float64)
mass = float(r[0]) if r.size > 0 else 0.0
izz = float(r[1]) if r.size > 1 else 0.0
# Layout matches PointMassRotary (UX, UY, UZ, ROTX, ROTY, ROTZ):
# 2-D in-plane mass on UX + UY only; IZZ on ROTZ only.
return np.diag([mass, mass, 0.0, 0.0, 0.0, izz])