visualdynamics.core.fem

fem

A beam finite element model, and the eigensolution it gives.

visualdynamics is an analysis toolset, so this is deliberately the smallest modeling capability that produces something worth analyzing: three-dimensional two-node beams, flat-shell plates — four-node rectangles (MITC4) and three-node triangles (MITC3) — and lumped masses, assembled into mass and stiffness matrices and solved for real normal modes. It exists because a demonstration needs a truth model — a dense analytical answer the measured one can be compared against — and because building one should not require reaching for another package.

There are three ways in. Build a structure member by member, as the example below does. Give a geometry's element groups their properties — a material and a thickness for an element group of plates, a material and a section for an element group of beams (GroupProperties) — and let from_geometry build every element as the element it is, which is how an exodus file means a structure and the way a user builds a model of a simple article. Or draw a shape as a surface mesh and let from_geometry, given one material and one section, put a member along every edge of it and share the mass over its nodes: the shorter road from any geometry to some set of modes, and how visualdynamics.demo.drone is built throughout.

What it is not: a general finite element code. The plate is rectangular only — a skewed or warped quad is refused rather than solved badly — and there are no curved shells, no solids, no constraints beyond fixing degrees of freedom, and no static solution. Wiring a surface mesh's edges with from_geometry makes a grillage of beams, not plates, which is a real modeling choice with a known cost rather than an approximation hidden inside an element: it answers what order of mode density and what mode families a shape has, and it does not pretend to be shell theory.

Everything here is SI, because the mass and stiffness matrices are the one place in visualdynamics where several dimensions have to be consistent with each other at once — a length in millimeters beside a modulus in pascals is not wrong in any single entry, it is wrong only in the answer. The objects that come out (Geometry, ShapeSet) carry their units declared, so the conversion happens once, at the boundary, as it does everywhere else.

from visualdynamics import fem

aluminum = fem.Material('aluminum', youngs_modulus=70e9, density=2700)
tube = fem.Section.round_tube('16 mm tube', outer=0.016, wall=0.001)

model = fem.Model('cantilever')
for i in range(11):
    model.add_node(100 + i, i * 0.1, 0.0, 0.0)
model.add_chain(range(100, 111), aluminum, tube)
shapes = model.eigensolution(maximum_frequency=2000,
                             fixed=['100'], damping=0.01)

The element formulation is Euler-Bernoulli with a consistent mass matrix: no shear flexibility and no section rotary inertia, which is right while a member is slender (length more than about ten times its depth) and progressively optimistic when it is not. A drone arm at 180 mm long and 16 mm deep is comfortably inside that; a stubby mounting stub is not, and frequencies there read high.

References

The Euler-Bernoulli beam element and its consistent mass matrix are textbook; these are where they are set out.

1. Przemieniecki, J. S. (1968). Theory of Matrix Structural Analysis. McGraw-Hill.
2. Cook, R. D., Malkus, D. S., Plesha, M. E., & Witt, R. J. (2002). Concepts and Applications of Finite Element Analysis, 4th ed. Wiley. The cubic Hermite shape functions the consistent mass follows from, and why consistent rather than lumped changes the frequencies returned.
3. Craig, R. R., & Kurdila, A. J. (2006). Fundamentals of Structural Dynamics, 2nd ed. Wiley. The generalized symmetric eigenproblem this solves, and normalization to unit modal mass.

Classes:

Name	Description
Material	An isotropic elastic material.
LibraryMaterial	A material the library offers, with where its numbers came from.
Section	A beam cross section, as the four numbers the element needs.
GroupProperties	What an element group of a geometry is made of, for Model.from_geometry.
Beam	One two-node beam element.
Plate	One four-node rectangular plate-bending element.
Triangle	One three-node triangular plate-bending element.
Solid	One solid element: a hexahedron, a wedge or a tetrahedron, told
RigidLink	Two nodes held rigidly together, with no mass of their own: the
Spring	A discrete spring between two degrees of freedom, or from one to
LumpedMass	A rigid item carried at a node: a motor, a battery, a camera.
Face	A cosmetic face: shading, not stiffness.
Model	Nodes, beams and lumped masses, and the modes they imply.

Functions:

Name	Description
ground_axes	A ground's directions as AXES names, in their order: True for
material	A library material by name — material('6061-T6') — the same
with_article	'a channel', 'an I-beam', 'an angle' — a shape named in a
angle_major_axis	Where an angle's major principal axis lies — the direction its
connected_pieces	The connected components of an adjacency map, largest first.
polish_modes	Refine a block of approximate modes until they are converged.

Classes

Material dataclass

Material(name: str, youngs_modulus: float, density: float, poissons_ratio: float = 0.3, modulus_of_rigidity: float | None = None)

An isotropic elastic material.

Shear modulus is derived from the modulus and Poisson's ratio unless it is given: for carbon fiber laminates the isotropic relation is a poor guess and the torsional stiffness it implies can be out by a factor of two, so the door is left open to state it.

Attributes:

Name	Type	Description
is_rigid	bool	Whether this is RIGID: a link, not a material. A modulus

Attributes

is_rigid property

is_rigid: bool

Whether this is RIGID: a link, not a material. A modulus left blank in the Element Groups table (None) is a material not yet finished, not a link.

LibraryMaterial dataclass

LibraryMaterial(material: Material, note: str)

A material the library offers, with where its numbers came from.

The note names the kind of source and what the values are typical of, because a handbook's typical modulus and density are what a modal solution wants and also not what a particular heat of a particular alloy measures: a person with the part's certification in hand types those numbers over the library's.

Section dataclass

Section(name: str, area: float, iy: float, iz: float, j: float, shape: str = '', dimensions: tuple[float, ...] = ())

A beam cross section, as the four numbers the element needs.

iy and iz are second moments of area about the element's own local y and z axes, and j is the torsion constant — St Venant's, which equals the polar second moment only for a circular section. For anything else it is smaller, and using the polar value overstates torsional stiffness; the constructors below carry the right formula for the shapes they build.

The shapes (SHAPES): a section built from one remembers it — shape and its dimensions in meters, in the order the constructor takes them — so a table can show the dimensions and a file or a session script rebuilds the section from them. One rule for which way a shape faces: its width or flanges run along local y and its depth or height along local z, so the orientation vector (which names local y) points across the flanges, and iy is the strong-axis bending of an I-beam or a channel. An angle is the exception, because its leg axes are not principal and the element has no product term: its iy and iz are its principal moments, and the orientation vector points along the major principal axis (angle_major_axis says where that is relative to the legs).

The element assumes the shear center is the centroid. For the symmetric shapes it is; a channel's and an angle's are not, and the twisting a load through their centroid causes is not in the model.

polar is the inertia term, always the true polar second moment (iy + iz), because rotary inertia about the axis is a property of where the material is and has nothing to do with warping.

Methods:

Name	Description
of_shape	The section of a named shape (SHAPES), given its dimensions
round_tube	A circular tube, given its outside diameter and wall thickness.
rod	A solid circular rod.
rectangle	A solid rectangle, width along local y and height along z.
rectangular_tube	A rectangular tube of one wall thickness, width along local y
square_tube	A square tube, given its outside width and wall thickness — a
i_beam	A doubly symmetric I-beam: overall depth along local z,
channel	A channel: overall depth along local z, the flanges
angle	An angle of one thickness. Its leg axes are not principal and

Methods:

of_shape classmethod

of_shape(name: str, shape: str, dimensions: Sequence[float]) -> Section

The section of a named shape (SHAPES), given its dimensions in meters in the order the shape's constructor takes them.

Parameters:

Name	Type	Description	Default
name	str	What to call the section.	required
shape	str	A key of SHAPES.	required
dimensions	sequence of float	The shape's dimensions, in meters.	required

Returns:

Type	Description
Section

Source code in src/visualdynamics/core/fem.py

@classmethod
def of_shape(cls, name: str, shape: str,
             dimensions: Sequence[float]) -> Section:
    """The section of a named shape (`SHAPES`), given its dimensions
    in meters in the order the shape's constructor takes them.

    Parameters
    ----------
    name : str
        What to call the section.
    shape : str
        A key of `SHAPES`.
    dimensions : sequence of float
        The shape's dimensions, in meters.

    Returns
    -------
    Section
    """
    if shape not in SHAPES:
        raise ValueError(f'{shape!r} is not a section shape: '
                         + ', '.join(SHAPES))
    method, names = SHAPES[shape]
    if len(dimensions) != len(names):
        raise ValueError(f'{with_article(shape)} takes {len(names)} '
                         'dimensions — ' + ', '.join(names))
    return getattr(cls, method)(name, *(float(v) for v in dimensions))

round_tube classmethod

round_tube(name: str, outer: float, wall: float) -> Section

A circular tube, given its outside diameter and wall thickness.

Source code in src/visualdynamics/core/fem.py

@classmethod
def round_tube(cls, name: str, outer: float, wall: float) -> Section:
    """A circular tube, given its outside diameter and wall thickness."""
    ro, ri = outer / 2.0, outer / 2.0 - wall
    if ri < 0 or wall <= 0 or outer <= 0:
        raise ValueError(f'{name}: a round tube needs a wall between 0 '
                         'and the radius')
    area = math.pi * (ro ** 2 - ri ** 2)
    i = math.pi * (ro ** 4 - ri ** 4) / 4.0
    # a closed circular section is the one case where torsion is the
    # polar moment exactly: it does not warp
    return cls(name, area, i, i, 2.0 * i, 'round tube',
               (float(outer), float(wall)))

rod classmethod

rod(name: str, diameter: float) -> Section

A solid circular rod.

Source code in src/visualdynamics/core/fem.py

@classmethod
def rod(cls, name: str, diameter: float) -> Section:
    """A solid circular rod."""
    section = cls.round_tube(name, diameter, diameter / 2.0)
    return cls(name, section.area, section.iy, section.iz, section.j,
               'rod', (float(diameter),))

rectangle classmethod

rectangle(name: str, width: float, height: float) -> Section

A solid rectangle, width along local y and height along z.

Source code in src/visualdynamics/core/fem.py

@classmethod
def rectangle(cls, name: str, width: float, height: float) -> Section:
    """A solid rectangle, `width` along local y and `height` along z."""
    if width <= 0 or height <= 0:
        raise ValueError(f'{name}: a rectangle needs a width and a height')
    area = width * height
    iz = width ** 3 * height / 12.0     # bending in the local x-y plane
    iy = width * height ** 3 / 12.0     # bending in the local x-z plane
    long, short = max(width, height), min(width, height)
    # St Venant's constant for a solid rectangle, from the exact series
    # solution of the Prandtl stress function (Timoshenko & Goodier,
    # Theory of Elasticity, sec. 109). The odd terms fall as n^-5, so a
    # hundred of them leave under a part in 1e10. This replaced Roark's
    # closed approximation (2026-09-23), which had been described as
    # good to 0.1% and is 0.45% off at an aspect ratio of 1.2.
    tail = sum(math.tanh(n * math.pi * long / (2.0 * short)) / n ** 5
               for n in range(1, 200, 2))
    j = long * short ** 3 / 3.0 * (
        1.0 - 192.0 / math.pi ** 5 * (short / long) * tail)
    return cls(name, area, iy, iz, j, 'rectangle',
               (float(width), float(height)))

rectangular_tube classmethod

rectangular_tube(name: str, width: float, height: float, wall: float) -> Section

A rectangular tube of one wall thickness, width along local y and height along z, both outside dimensions.

Source code in src/visualdynamics/core/fem.py

@classmethod
def rectangular_tube(cls, name: str, width: float, height: float,
                     wall: float) -> Section:
    """A rectangular tube of one wall thickness, `width` along local y
    and `height` along z, both outside dimensions."""
    inner_w, inner_h = width - 2.0 * wall, height - 2.0 * wall
    if wall <= 0 or inner_w <= 0 or inner_h <= 0:
        raise ValueError(f'{name}: wall {wall} closes the section')
    area = width * height - inner_w * inner_h
    iy = (width * height ** 3 - inner_w * inner_h ** 3) / 12.0
    iz = (width ** 3 * height - inner_w ** 3 * inner_h) / 12.0
    # Bredt's thin-wall formula, J = 4 A_m^2 t / s: A_m the area inside
    # the wall's centreline, s that centreline's length
    mean_w, mean_h = width - wall, height - wall
    j = 4.0 * (mean_w * mean_h) ** 2 * wall / (2.0 * (mean_w + mean_h))
    return cls(name, area, iy, iz, j, 'rectangular tube',
               (float(width), float(height), float(wall)))

square_tube classmethod

square_tube(name: str, width: float, wall: float) -> Section

A square tube, given its outside width and wall thickness — a rectangular tube of equal sides.

Source code in src/visualdynamics/core/fem.py

@classmethod
def square_tube(cls, name: str, width: float, wall: float) -> Section:
    """A square tube, given its outside width and wall thickness — a
    rectangular tube of equal sides."""
    return cls.rectangular_tube(name, width, width, wall)

i_beam classmethod

i_beam(name: str, depth: float, flange_width: float, flange_thickness: float, web_thickness: float) -> Section

A doubly symmetric I-beam: overall depth along local z, flanges flange_width wide along y. Fillets are left out, which puts area and torsion a few percent under a rolled shape's table values (a W8x31: 8.99 in² against 9.13, J 0.50 in⁴ against 0.54).

Source code in src/visualdynamics/core/fem.py

@classmethod
def i_beam(cls, name: str, depth: float, flange_width: float,
           flange_thickness: float, web_thickness: float) -> Section:
    """A doubly symmetric I-beam: overall `depth` along local z,
    flanges `flange_width` wide along y. Fillets are left out, which
    puts area and torsion a few percent under a rolled shape's table
    values (a W8x31: 8.99 in² against 9.13, J 0.50 in⁴ against 0.54)."""
    d, bf, tf, tw = depth, flange_width, flange_thickness, web_thickness
    web = d - 2.0 * tf
    if min(d, bf, tf, tw) <= 0 or web <= 0 or tw > bf:
        raise ValueError(f'{name}: the flanges and web do not make an I')
    area = 2.0 * bf * tf + web * tw
    iy = (bf * d ** 3 - (bf - tw) * web ** 3) / 12.0
    iz = 2.0 * tf * bf ** 3 / 12.0 + web * tw ** 3 / 12.0
    # an open thin-walled section: the sum of b t^3 / 3 over its
    # plates, the web taken between the flanges' mid-planes
    j = (2.0 * bf * tf ** 3 + (d - tf) * tw ** 3) / 3.0
    return cls(name, area, iy, iz, j, 'I-beam',
               (float(d), float(bf), float(tf), float(tw)))

channel classmethod

channel(name: str, depth: float, flange_width: float, flange_thickness: float, web_thickness: float) -> Section

A channel: overall depth along local z, the flanges flange_width wide (web included) running along +y from the web. iz is about the centroid, which sits off the web; the shear center sits further off it the other way, and the element does not know. The flanges are of one thickness: a rolled channel's taper toward their tips, and its table flange thickness is an average, so this iz runs high for one (a C6x10.5: 1.06 in⁴ against the table's 0.86) — a bent-plate channel it gives exactly.

Source code in src/visualdynamics/core/fem.py

@classmethod
def channel(cls, name: str, depth: float, flange_width: float,
            flange_thickness: float, web_thickness: float) -> Section:
    """A channel: overall `depth` along local z, the flanges
    `flange_width` wide (web included) running along +y from the web.
    `iz` is about the centroid, which sits off the web; the shear
    center sits further off it the other way, and the element does
    not know. The flanges are of one thickness: a rolled channel's
    taper toward their tips, and its table flange thickness is an
    average, so this `iz` runs high for one (a C6x10.5: 1.06 in⁴
    against the table's 0.86) — a bent-plate channel it gives
    exactly."""
    d, bf, tf, tw = depth, flange_width, flange_thickness, web_thickness
    web = d - 2.0 * tf
    if min(d, bf, tf, tw) <= 0 or web <= 0 or tw > bf:
        raise ValueError(f'{name}: the flanges and web do not make a '
                         'channel')
    flange_area, web_area = bf * tf, web * tw
    area = 2.0 * flange_area + web_area
    iy = (bf * d ** 3 - (bf - tw) * web ** 3) / 12.0
    # centroid along y, from the back of the web
    centroid = (2.0 * flange_area * bf / 2.0 + web_area * tw / 2.0) / area
    iz = (2.0 * (tf * bf ** 3 / 12.0
                 + flange_area * (bf / 2.0 - centroid) ** 2)
          + web * tw ** 3 / 12.0 + web_area * (tw / 2.0 - centroid) ** 2)
    j = (2.0 * bf * tf ** 3 + (d - tf) * tw ** 3) / 3.0
    return cls(name, area, iy, iz, j, 'channel',
               (float(d), float(bf), float(tf), float(tw)))

angle classmethod

angle(name: str, long_leg: float, short_leg: float, thickness: float) -> Section

An angle of one thickness. Its leg axes are not principal and the element has no product term, so iy and iz are the principal moments — major and minor — and the orientation vector is to point along the major principal axis (angle_major_axis).

Source code in src/visualdynamics/core/fem.py

@classmethod
def angle(cls, name: str, long_leg: float, short_leg: float,
          thickness: float) -> Section:
    """An angle of one thickness. Its leg axes are not principal and
    the element has no product term, so `iy` and `iz` are the
    principal moments — major and minor — and the orientation vector
    is to point along the major principal axis (`angle_major_axis`)."""
    a, b, t = long_leg, short_leg, thickness
    if min(a, b, t) <= 0 or t >= min(a, b) or b > a:
        raise ValueError(f'{name}: an angle needs a long leg, a short '
                         'leg no longer, and a thickness under both')
    area, (i_major, i_minor), _theta = _angle_properties(a, b, t)
    j = (a + b - t) * t ** 3 / 3.0
    return cls(name, area, i_major, i_minor, j, 'angle',
               (float(a), float(b), float(t)))

GroupProperties dataclass

GroupProperties(material: Material | None = None, thickness: float | None = None, section: Section | None = None, orientation: tuple[float, float, float] | None = None, mass: float | None = None, stiffness: tuple[float | None, ...] | None = None, ground: tuple[str, ...] | bool = ())

What an element group of a geometry is made of, for Model.from_geometry.

One property set per element group of one element type, the way every finite element format states a structure: a material for any element group, plus a thickness for an element group of plates (triangles or quads) or a section — and an orientation vector for the roll, as Model.add_beam takes it — for an element group of beams; an element group of solids takes the material alone. An element group given both, or plates or beams given neither, is refused when the model is built, by name.

An element group of point elements is a set of lumped masses and takes a mass alone, no material (2026-10-07): GroupProperties(mass=m) puts m kilograms at the node of every element in the element group — a bolt, a sensor, a fitting too small to mesh — which is how a finite element deck states one, a CONM2 or a point-mass element group.

Springs and ground (2026-10-08, Brandon: a spring is a kind of line element, and ground a kind of point). A group given a stiffness, six numbers in the global X, Y, Z, RX, RY and RZ directions with None for free, is a group of springs: every two-node line in it a spring between its two nodes in each direction given (Model.add_spring), the nodes free to coincide, and every point element a spring from its node to ground — a Nastran CBUSH, or CELAS cards one direction at a time. A group given ground holds points, and each one's node is held in each direction given (Model.add_ground): all six with ground=True, the translations alone with ground=('X', 'Y', 'Z') — a support, or the far end of a spring to ground drawn as a line. A Nastran SPC1 reads as one, its components the directions.

Attributes:

Name	Type	Description
kind	str	'plate', 'beam', 'solid', 'rigid', 'mass', or what is wrong

Attributes

kind property

kind: str

'plate', 'beam', 'solid', 'rigid', 'mass', or what is wrong with it. A rigid element group takes no thickness and no section, and any left from before the material was picked are ignored; a material alone is an element group of solids (2026-09-30), which take nothing else — an element group of plates or beams given only a material is refused where the elements are built, by what they are. A mass makes an element group of point masses whatever else is set: the mass is the one number such an element group can use; ground and a stiffness likewise make a group of supports and of springs.

Beam dataclass

Beam(node_a: int, node_b: int, material: Material, section: Section, orientation: tuple[float, float, float] | None = None, color: int = 1, group: str = '')

One two-node beam element.

Plate dataclass

Plate(nodes: tuple[int, int, int, int], material: Material, thickness: float, color: int = 1, group: str = '')

One four-node rectangular plate-bending element.

A flat shell: plane-stress membrane action in its own plane and Mindlin bending out of it, with the transverse shear tied at the edge midpoints (MITC4). The tying is not optional finesse — a plain bilinear Mindlin element locks in shear as the plate gets thin, and a 12x12 mesh of the locked element puts the first elastic mode of a thin free plate several times too high.

Rectangles only, and refused otherwise rather than silently mis-integrated: the tying directions assume the natural axes align with the sides, which is exactly true for a rectangle and only approximately for anything else. The models this module exists to build mesh rectangular panels; a skewed general quad earns its place when something needs it, with the covariant transforms and the validation that come with it.

Nodes run around the perimeter: 1-2 is the first edge, 1-4 the second, corner 3 opposite corner 1.

Triangle dataclass

Triangle(nodes: tuple[int, int, int], material: Material, thickness: float, color: int = 1, group: str = '')

One three-node triangular plate-bending element.

The quad's sibling, so a mesh of triangles and a mesh of rectangles are one theory: a flat shell with plane-stress membrane action in its own plane and Mindlin bending out of it, the transverse shear tied along the three edges (MITC3, Lee & Bathe 2004). The tying is what keeps a linear triangle from locking in shear as the plate gets thin — a plain linear Mindlin triangle is the worst locker there is. Any flat triangle is a valid element; only a degenerate one (zero area) is refused.

Nodes run around the perimeter, counterclockwise about the normal the element takes as its own +z.

Solid dataclass

Solid(nodes: tuple[int, ...], material: Material, color: int = 1, group: str = '')

One solid element: a hexahedron, a wedge or a tetrahedron, told apart by how many nodes it names (8, 6 or 4).

Three translations per node and no rotations — a solid has no rotational stiffness, and the rotations of a node only solids touch are grounded by the eigensolution rather than left as degrees of freedom with nothing on them (Model.dangling_rotations).

The hexahedron is trilinear with Wilson's incompatible bending modes, Taylor's form (the extra modes' strains taken from the centroid's Jacobian, so a distorted brick still passes the patch test): a plain trilinear brick is far too stiff in bending, and a part meshed a few elements through its thickness would come out a third high. With the modes, one layer of bricks bends like a beam. The wedge is the linear six-node element and the tetrahedron the constant-strain one — transition and imported shapes, not what a part is meshed with here (mesh.block makes bricks); a linear tetrahedron locks in bending and a mesh of them is trusted only where it is fine.

Nodes run around the bottom face and then the top, the same way round (UFF 2412 and Nastran's CHEXA/CPENTA/CTETRA order).

RigidLink dataclass

RigidLink(node_a: int, node_b: int, group: str = '')

Two nodes held rigidly together, with no mass of their own: the second moves as the first does, translated by the first's rotation about it (RIGID).

Spring dataclass

Spring(node_a: int, direction_a: int, node_b: int | None, direction_b: int, stiffness: float, name: str = '')

A discrete spring between two degrees of freedom, or from one to ground: the flexible counterpart of a RigidLink, and what a Nastran CELAS states (2026-10-08, proposed for two-beam substructuring cases, where a test article meets its fixture at coincident nodes with no length between them for a beam).

Each end is a node and a signed direction code (direction_code: 1-3 the translations, 4-6 the rotations); the spring resists the difference of the two motions, each taken along its own sign, so 'X+' to 'X+' is the ordinary axial spring. node_b None grounds it.

LumpedMass dataclass

LumpedMass(node: int, mass: float, inertia: tuple[float, float, float] = (0.0, 0.0, 0.0), name: str = '')

A rigid item carried at a node: a motor, a battery, a camera.

The inertias are about the global axes through the node. A point mass leaves them zero, which is honest for something small against the members carrying it and wrong for a battery the size of the structure — a mass with no inertia cannot rock, so a rocking mode simply will not appear.

Face dataclass

Face(nodes: tuple[int, ...], color: int = 1, group: str = '')

A cosmetic face: shading, not stiffness.

A grillage of beams reads as a wireframe, and a wireframe of a drone deck reads as nothing much. Faces spanning nodes that are already there give the renderer something to shade and the animation something to deform, while contributing nothing to the matrices — which is exactly the truth about them, and is why they are a separate kind rather than a zero-stiffness element.

Model

Model(name: str = '', length_unit: str = 'm')

Nodes, beams and lumped masses, and the modes they imply.

Two assemblies of one set of element matrices. Dense (matrices): (6 x nodes) square, every mode solved whole — exact, and the path for a model up to SPARSE_ABOVE degrees of freedom. Sparse (sparse_matrices): only the nonzeros, the lowest modes by shift-invert Lanczos. The dense path was the only one, and its ceiling (~1000 nodes) deliberate, until a plate model of a small real structure met it: the BARC at a quarter inch, 1,500 nodes, is ~13 GB dense; sparse, its eighth-inch mesh of 6,150 nodes solves in three seconds in under half a gigabyte (2026-09-26).

Attributes: name: What the model is called; it becomes the geometry's name and the shape set's comment. length_unit: What the coordinates are in. Everything inside is SI; this is what the Geometry is told on the way out. beams: Every member, as Beam records naming two nodes, a material and a section. masses: Lumped masses, as LumpedMass records at a node. springs: Discrete springs between two degrees of freedom or to ground, as Spring records (add_spring). grounds: The nodes held, and in which directions, {node: axes 0-5} (add_ground). faces: Surfaces, as Face records naming three or four nodes. They carry no stiffness — a face is drawn, and its edges are what carry members.

Methods:

Name	Description
add_chain	A run of beams through consecutive nodes — a member, in one call.
add_plate	One rectangular plate element over four existing nodes.
add_triangle	One triangular plate element over three existing nodes.
add_solid	One solid element over eight, six or four existing nodes: a
add_rigid_link	Join two nodes rigidly, adding no mass.
add_spring	A spring between two degrees of freedom, or from one to ground.
add_ground	Hold a node in some or all of the six directions: a support,
rigid_bodies	The groups of nodes the rigid links join, each in the model's
from_geometry	A structure from a drawn shape: members on its edges, mass at its
pieces	The structure's disconnected parts, largest first.
wire_faces	Put a beam along every edge of every face, and return how many.
distribute_mass	Share total over the nodes by the members meeting at each.
group	What this node was added as part of — 'arm 2', 'top deck'.
dof_strings	'101X+', '101Y+', … in the matrices' own order.
matrices	Assemble the global mass and stiffness matrices, dense.
sparse_matrices	The same mass and stiffness matrices as matrices, stored
rigid_body_vectors	The six rigid-body motions, as columns over the model's DOFs.
eigensolution	Real normal modes, mass-normalized, as a ShapeSet.
scaled_system	The sparse eigenproblem as the sparse solver poses it.
constraint_transform	T, with u = T q: every degree of freedom of the model written
dangling_rotations	The nodes whose rotations nothing acts on: touched by solids or
idle_lead_rotations	(lead node, axis 0-2) for each rotation of a rigid body's lead
loose_nodes	The nodes nothing touches: no element, no rigid link, no
geometry	The model as a Geometry: nodes, beams as elements, faces as faces.

Attributes:

Name	Type	Description
node_ids	list[int]	In insertion order: the matrices' row order follows this.
structural_mass	float	What the members weigh, before anything is hung on them.

Source code in src/visualdynamics/core/fem.py

def __init__(self, name: str = '', length_unit: str = 'm') -> None:
    self.name: str = name
    #: everything here is SI; this is what the Geometry is told
    self.length_unit: str = length_unit
    self._nodes: dict[int, np.ndarray] = {}
    self._node_group: dict[int, str] = {}
    self.beams: list[Beam] = []
    self.plates: list[Plate] = []
    self.triangles: list[Triangle] = []
    self.solids: list[Solid] = []
    self.rigid_links: list[RigidLink] = []
    self.springs: list[Spring] = []
    #: {node: the axes held, 0-5}
    self.grounds: dict[int, tuple[int, ...]] = {}
    self.masses: list[LumpedMass] = []
    self.faces: list[Face] = []

Attributes

node_ids property

node_ids: list[int]

In insertion order: the matrices' row order follows this.

structural_mass property

structural_mass: float

What the members weigh, before anything is hung on them.

Methods:

add_chain

add_chain(nodes: Sequence[int], material: Material, section: Section, orientation: Sequence[float] | None = None, color: int = 1, group: str = '') -> list[Beam]

A run of beams through consecutive nodes — a member, in one call.

Source code in src/visualdynamics/core/fem.py

def add_chain(self, nodes: Sequence[int], material: Material,
              section: Section,
              orientation: Sequence[float] | None = None,
              color: int = 1, group: str = '') -> list[Beam]:
    """A run of beams through consecutive nodes — a member, in one call."""
    nodes = [int(n) for n in nodes]
    return [self.add_beam(a, b, material, section, orientation, color, group)
            for a, b in pairwise(nodes)]

add_plate

add_plate(nodes: Sequence[int], material: Material, thickness: float, color: int = 1, group: str = '') -> Plate

One rectangular plate element over four existing nodes.

Source code in src/visualdynamics/core/fem.py

def add_plate(self, nodes: Sequence[int], material: Material,
              thickness: float, color: int = 1,
              group: str = '') -> Plate:
    """One rectangular plate element over four existing nodes."""
    nodes = tuple(int(n) for n in nodes)
    if len(nodes) != 4 or len(set(nodes)) != 4:
        raise ValueError('a plate spans four distinct nodes')
    for node in nodes:
        if node not in self._nodes:
            raise ValueError(
                f'plate names node {node}, which is not in the model')
    if float(thickness) <= 0.0:
        raise ValueError('a plate needs a positive thickness')
    # validate the rectangle at build time, not at solve time: the
    # person holding the bad corner coordinate is the one adding it
    _plate_frame(*[self._nodes[n] for n in nodes])
    plate = Plate(nodes, material, float(thickness), color, group)
    self.plates.append(plate)
    return plate

add_triangle

add_triangle(nodes: Sequence[int], material: Material, thickness: float, color: int = 1, group: str = '') -> Triangle

One triangular plate element over three existing nodes.

Source code in src/visualdynamics/core/fem.py

def add_triangle(self, nodes: Sequence[int], material: Material,
                 thickness: float, color: int = 1,
                 group: str = '') -> Triangle:
    """One triangular plate element over three existing nodes."""
    nodes = tuple(int(n) for n in nodes)
    if len(nodes) != 3 or len(set(nodes)) != 3:
        raise ValueError('a triangle spans three distinct nodes')
    for node in nodes:
        if node not in self._nodes:
            raise ValueError(
                f'triangle names node {node}, which is not in the model')
    if float(thickness) <= 0.0:
        raise ValueError('a triangle needs a positive thickness')
    _triangle_frame(*[self._nodes[n] for n in nodes])
    triangle = Triangle(nodes, material, float(thickness), color, group)
    self.triangles.append(triangle)
    return triangle

add_solid

add_solid(nodes: Sequence[int], material: Material, color: int = 1, group: str = '') -> Solid

One solid element over eight, six or four existing nodes: a hexahedron, a wedge or a tetrahedron.

Parameters:

Name	Type	Description	Default
nodes	sequence of int	The corners, bottom face then top, the same way round.	required
material	Material	What it is made of.	required
color	optional	As for a plate.	1
group	optional	As for a plate.	1

Returns:

Type	Description
Solid

Source code in src/visualdynamics/core/fem.py

def add_solid(self, nodes: Sequence[int], material: Material,
              color: int = 1, group: str = '') -> Solid:
    """One solid element over eight, six or four existing nodes: a
    hexahedron, a wedge or a tetrahedron.

    Parameters
    ----------
    nodes : sequence of int
        The corners, bottom face then top, the same way round.
    material : Material
        What it is made of.
    color, group : optional
        As for a plate.

    Returns
    -------
    Solid
    """
    nodes = tuple(int(n) for n in nodes)
    if len(nodes) not in (8, 6, 4) or len(set(nodes)) != len(nodes):
        raise ValueError('a solid spans eight, six or four distinct nodes')
    for node in nodes:
        if node not in self._nodes:
            raise ValueError(
                f'solid names node {node}, which is not in the model')
    if material.is_rigid:
        raise ValueError('a solid is made of a material, not rigid')
    # a flat element is caught here, by the person holding the bad
    # coordinate, not by the eigensolver. One numbered the other way
    # round (its volume negative) is the same element and is taken:
    # meshers disagree on the handedness, and Linderholt's frame
    # mesh arrived inside out to this convention (2026-09-30)
    xyz = np.array([self._nodes[n] for n in nodes])
    if abs(_solid_volume(xyz)) <= 1e-12 * float(np.ptp(xyz)) ** 3:
        raise ValueError(f'solid {nodes} has no volume')
    solid = Solid(nodes, material, color, group)
    self.solids.append(solid)
    return solid

add_rigid_link

add_rigid_link(node_a: int, node_b: int, group: str = '') -> RigidLink

Join two nodes rigidly, adding no mass.

Links that share nodes join into one rigid body, however they are chained; each body moves as its first node does (the one added to the model first), and the others follow it exactly. The eigensolution eliminates the followers' degrees of freedom rather than stiffening anything, so the answer is the limit of an infinitely stiff member and the matrices stay well conditioned.

Parameters:

Name	Type	Description	Default
node_a	int	The nodes, both already in the model, and different.	required
node_b	int	The nodes, both already in the model, and different.	required
group	str	The part the link belongs to — its element group, from a geometry.	''

Returns:

Type	Description
RigidLink

Source code in src/visualdynamics/core/fem.py

def add_rigid_link(self, node_a: int, node_b: int,
                   group: str = '') -> RigidLink:
    """Join two nodes rigidly, adding no mass.

    Links that share nodes join into one rigid body, however they are
    chained; each body moves as its first node does (the one added to
    the model first), and the others follow it exactly. The
    eigensolution eliminates the followers' degrees of freedom rather
    than stiffening anything, so the answer is the limit of an
    infinitely stiff member and the matrices stay well conditioned.

    Parameters
    ----------
    node_a, node_b : int
        The nodes, both already in the model, and different.
    group : str, optional
        The part the link belongs to — its element group, from a geometry.

    Returns
    -------
    RigidLink
    """
    for node in (node_a, node_b):
        if int(node) not in self._nodes:
            raise ValueError(f'rigid link names node {node}, which is not '
                             'in the model')
    if int(node_a) == int(node_b):
        raise ValueError(f'a rigid link joins two nodes; it names node '
                         f'{node_a} twice')
    link = RigidLink(int(node_a), int(node_b), group)
    self.rigid_links.append(link)
    return link

add_spring

add_spring(dof_a: str, dof_b: str | None, stiffness: float, name: str = '') -> Spring

A spring between two degrees of freedom, or from one to ground.

The ends are written the way fixed writes them, a node and a direction: model.add_spring('101Z+', '201Z+', 5e5) joins two nodes' Z translations, model.add_spring('101RY+', '201RY+', 2e3) their rotations about Y, and model.add_spring('1Z+', None, 1e4) grounds node 1 in Z. The nodes may coincide — a joint between two parts meshed to the same point, which no beam can be — and a joint stiff in more than one direction is one spring per direction. It adds no mass. A spring to a node nothing else touches is a spring to ground, since that node is grounded whole (loose_nodes).

A spring far stiffer than the structure costs the solve digits: two beams joined by springs ten orders over their EI/L read 0.2 % low on the dense solver (2026-10-08), where five orders were within 3e-5 of the one beam. A joint meant to be rigid is a rigid link (add_rigid_link), which costs nothing.

Parameters:

Name	Type	Description	Default
dof_a	str	The first end, as '': 'X+' to 'Z+' a translation, 'RX+' to 'RZ+' a rotation; a minus sign turns the end around.	required
dof_b	str or None	The second end, the same way, or None for ground. Both ends are translations or both rotations.	required
stiffness	float	Positive: N/m between translations, N m/rad between rotations.	required
name	str	What the spring is called — its element group, from a geometry.	''

Returns:

Type	Description
Spring

Source code in src/visualdynamics/core/fem.py

def add_spring(self, dof_a: str, dof_b: str | None, stiffness: float,
               name: str = '') -> Spring:
    """A spring between two degrees of freedom, or from one to ground.

    The ends are written the way `fixed` writes them, a node and a
    direction: ``model.add_spring('101Z+', '201Z+', 5e5)`` joins two
    nodes' Z translations, ``model.add_spring('101RY+', '201RY+',
    2e3)`` their rotations about Y, and ``model.add_spring('1Z+',
    None, 1e4)`` grounds node 1 in Z. The nodes may coincide — a
    joint between two parts meshed to the same point, which no beam
    can be — and a joint stiff in more than one direction is one
    spring per direction. It adds no mass. A spring to a node
    nothing else touches is a spring to ground, since that node is
    grounded whole (`loose_nodes`).

    A spring far stiffer than the structure costs the solve digits:
    two beams joined by springs ten orders over their EI/L read
    0.2 % low on the dense solver (2026-10-08), where five orders
    were within 3e-5 of the one beam. A joint meant to be rigid is
    a rigid link (`add_rigid_link`), which costs nothing.

    Parameters
    ----------
    dof_a : str
        The first end, as '<node><direction>': 'X+' to 'Z+' a
        translation, 'RX+' to 'RZ+' a rotation; a minus sign turns
        the end around.
    dof_b : str or None
        The second end, the same way, or None for ground. Both ends
        are translations or both rotations.
    stiffness : float
        Positive: N/m between translations, N m/rad between
        rotations.
    name : str, optional
        What the spring is called — its element group, from a geometry.

    Returns
    -------
    Spring
    """
    ends = [self._spring_end(text) for text in (dof_a, dof_b)
            if text is not None]
    if not float(stiffness) > 0.0:
        raise ValueError(f'a spring\'s stiffness is positive; it was '
                         f'given {stiffness!r}')
    if len(ends) == 2:
        (node_a, code_a), (node_b, code_b) = ends
        if (abs(code_a) > 3) != (abs(code_b) > 3):
            raise ValueError(f'{dof_a!r} and {dof_b!r}: a spring joins two '
                             'translations or two rotations')
        if node_a == node_b and abs(code_a) == abs(code_b):
            raise ValueError(f'{dof_a!r} and {dof_b!r} are one degree of '
                             'freedom; a spring joins two, or one to '
                             'ground')
    else:
        (node_a, code_a), (node_b, code_b) = ends[0], (None, 0)
    spring = Spring(node_a, code_a, node_b, code_b, float(stiffness), name)
    self.springs.append(spring)
    return spring

add_ground

add_ground(node: int, directions: Any = True) -> int

Hold a node in some or all of the six directions: a support, solved for as fixed would hold it, but carried by the model — what a ground point in a geometry builds. A node nothing else touches is held whole already (loose_nodes); this holds one the structure touches too. Held twice, a node is held in both sets.

Parameters:

Name	Type	Description	Default
node	int	The node, already in the model.	required
directions	bool or sequence of str	The directions held, as AXES names ('X', 'RY'); True for all six.	True

Returns:

Type	Description
int	The node.

Source code in src/visualdynamics/core/fem.py

def add_ground(self, node: int, directions: Any = True) -> int:
    """Hold a node in some or all of the six directions: a support,
    solved for as `fixed` would hold it, but carried by the model —
    what a ground point in a geometry builds. A node nothing else
    touches is held whole already (`loose_nodes`); this holds one the
    structure touches too. Held twice, a node is held in both sets.

    Parameters
    ----------
    node : int
        The node, already in the model.
    directions : bool or sequence of str, default True
        The directions held, as `AXES` names ('X', 'RY'); True for
        all six.

    Returns
    -------
    int
        The node.
    """
    if int(node) not in self._nodes:
        raise ValueError(f'ground names node {node}, which is not in '
                         'the model')
    axes = {AXES.index(axis) for axis in ground_axes(directions)}
    if not axes:
        raise ValueError(f'a ground at node {node} holds no direction')
    held = set(self.grounds.get(int(node), ())) | axes
    self.grounds[int(node)] = tuple(sorted(held))
    return int(node)

rigid_bodies

rigid_bodies() -> list[list[int]]

The groups of nodes the rigid links join, each in the model's node order — its first node is the one the others follow.

Source code in src/visualdynamics/core/fem.py

def rigid_bodies(self) -> list[list[int]]:
    """The groups of nodes the rigid links join, each in the model's
    node order — its first node is the one the others follow."""
    neighbors: dict[int, set[int]] = {}
    for link in self.rigid_links:
        neighbors.setdefault(link.node_a, set()).add(link.node_b)
        neighbors.setdefault(link.node_b, set()).add(link.node_a)
    order = {node: i for i, node in enumerate(self._nodes)}
    return [sorted(piece, key=order.__getitem__)
            for piece in connected_pieces(neighbors)]

from_geometry classmethod

from_geometry(geometry: Geometry, material: Material | None = None, section: Section | None = None, total_mass: float | None = None, name: str = '', groups: dict[int, str] | None = None, sections: dict[str, Section] | None = None) -> Model

A structure from a drawn shape: members on its edges, mass at its nodes.

The shortest route from any geometry to a set of modes. Every element contributes its own edges as beams — a face gives its perimeter, a line element gives its run — and the mass is shared equally over the nodes. What comes back is a model that can be solved like any other.

This is a sanity-check tool, not a mesher. The members are a stand-in for whatever the real structure is: one section for the whole model, chosen to put the modes where they are wanted, and the answer scales as its square root. Shell bending, membrane action and any real thickness are simply not represented. What it is good for is looking at a shape and asking what order of mode density and what mode families it has — which is the question a display model usually raises first.

A geometry whose element groups carry properties builds itself. When geometry.group_properties names what each element group is made of (GroupProperties), every element becomes the element it is: a quad a plate, a triangle a triangle, a two-node line a beam, each with its element group's material and thickness or section, and a point element in an element group given a mass a lumped mass at its node. That is the model an exodus file means, one property set per element group of one element type (Brandon, 2026-09-25), and the demonstration plate rebuilt from its own geometry this way is the same model to the last digit (tests/test_block_model.py). An element group with no properties, or an element type the solver has no element for, is refused by name; material and section are not consulted. Without element group properties the grillage below is built, and material and section are required for it.

A drawn line — an element group of two-node line elements with no properties, what a traceline was — is an element like any other here and gives its run, which is what a wireframe geometry needs to hold together at all. groups labels the nodes by the part they belong to; a Geometry does not carry that, and the first question asked of any result is which part of the structure a mode lives in. sections then gives one part a section of its own — {'prop': stiffer} — applied where both ends of an edge belong to it, which is how a part is made stiffer or softer than the rest without redrawing anything.

Source code in src/visualdynamics/core/fem.py

@classmethod
def from_geometry(cls, geometry: Geometry, material: Material | None = None,
                  section: Section | None = None,
                  total_mass: float | None = None,
                  name: str = '',
                  groups: dict[int, str] | None = None,
                  sections: dict[str, Section] | None = None) -> Model:
    """A structure from a drawn shape: members on its edges, mass at its
    nodes.

    The shortest route from *any* geometry to a set of modes. Every
    element contributes its own edges as beams — a face gives its
    perimeter, a line element gives its run — and the mass is shared
    equally over the nodes. What comes back is a model that can be
    solved like any other.

    This is a sanity-check tool, not a mesher. The members are a stand-in
    for whatever the real structure is: one section for the whole model,
    chosen to put the modes where they are wanted, and the answer scales
    as its square root. Shell bending, membrane action and any real
    thickness are simply not represented. What it *is* good for is
    looking at a shape and asking what order of mode density and what
    mode families it has — which is the question a display model usually
    raises first.

    **A geometry whose element groups carry properties builds itself.** When
    `geometry.group_properties` names what each element group is made of
    (`GroupProperties`), every element becomes the element it is:
    a quad a plate, a triangle a triangle, a two-node line a beam,
    each with its element group's material and thickness or section, and a
    point element in an element group given a mass a lumped mass at its
    node. That
    is the model an exodus file means, one property set per element group
    of one element type (Brandon, 2026-09-25), and the demonstration
    plate rebuilt from its own geometry this way is the same model
    to the last digit (`tests/test_block_model.py`). An element group with no
    properties, or an element type the solver has no element for,
    is refused by name; `material` and `section` are not consulted.
    Without element group properties the grillage below is built, and
    `material` and `section` are required for it.

    A drawn line — an element group of two-node line elements with no
    properties, what a traceline was — is an element like any other
    here and gives its run, which is what a wireframe geometry
    needs to hold together at all.
    `groups` labels the nodes by the part they belong to; a Geometry
    does not carry that, and the first question asked of any result is
    which part of the structure a mode lives in. `sections` then gives
    one part a section of its own — {'prop': stiffer} — applied where
    both ends of an edge belong to it, which is how a part is made
    stiffer or softer than the rest without redrawing anything.
    """
    model = cls(name or getattr(geometry, 'name', '') or 'geometry',
                length_unit=geometry.length_unit or 'm')
    properties = dict(getattr(geometry, 'group_properties', {}) or {})
    if properties:
        return _from_groups(model, geometry, properties, total_mass,
                            groups)
    if material is None or section is None:
        raise ValueError(
            'the geometry carries no element group properties, so a material '
            'and a section are needed to make a grillage of it — or '
            'give each element group its properties (fem.GroupProperties) and '
            'the elements build themselves')
    labels = dict(groups or {})
    if not labels:
        # A geometry that carries element groups says for itself which
        # part each node belongs to, so nothing has to be passed
        # alongside it. That matters because a side-channel does not
        # survive being saved: a geometry written to a file and read
        # back could not reproduce the model it came from.
        labels = _labels_from_groups(geometry)
    for node, xyz in zip(geometry.node_id, geometry.node_xyz):
        model.add_node(int(node), *[float(v) for v in xyz],
                       group=labels.get(int(node), ''))

    # A member takes its section from the element group of the element it came
    # from, not from labels on its end nodes. A node on a seam belongs
    # to two parts and can only answer for one, which left 24 of the
    # drone's 1248 blade members reading as ordinary frame; an element
    # belongs to exactly one element group and is never ambiguous.
    parts = _group_labels(geometry)
    runs: list[tuple[list[int], str]] = []
    for index, (kind, conn) in enumerate(zip(geometry.elem_type,
                                             geometry.elem_conn)):
        nodes = [int(n) for n in conn]
        label = parts[index] if index < len(parts) else ''
        shape = ELEMENT_TYPES.get(int(kind), (None, 0, 'line'))[2]
        if shape == 'line':
            runs.append((nodes, label))
        else:
            # a face closes on itself; three or four of them also
            # become a Face, so the shape can be drawn back
            runs.append((nodes + nodes[:1], label))
            if len(nodes) in (3, 4):
                model.add_face(nodes, group=label)
    seen: set[frozenset[int]] = set()
    for run, label in runs:
        chosen = section
        for part, alternative in (sections or {}).items():
            if label.startswith(part):
                chosen = alternative
                break
        for a, b in pairwise(run):
            edge = frozenset((a, b))
            if a == b or edge in seen:
                continue
            seen.add(edge)
            model.add_beam(a, b, material, chosen, group=label)
    if not seen:
        raise ValueError(
            'the geometry has no elements to make members '
            'from, so there is nothing to connect its nodes')

    loose = sorted(set(model.node_ids)
                   - {n for edge in seen for n in edge})
    if loose:
        raise ValueError(
            f'{len(loose)} nodes are connected to nothing, so they would '
            'carry mass with no stiffness and the solution would not '
            'factorize: ' + ', '.join(str(n) for n in loose[:8]))
    pieces = model.pieces()
    if len(pieces) > 1:
        sizes = ', '.join(str(len(p)) for p in pieces[:6])
        raise ValueError(
            f'the geometry is {len(pieces)} disconnected pieces ({sizes} '
            'nodes), which would solve as that many free bodies and give '
            f'{6 * len(pieces)} zero-frequency modes rather than 6. Its '
            'elements do not join them: either they are '
            'meant to be separate structures, or the connectivity is '
            'incomplete')
    if total_mass is not None:
        model.distribute_mass(total_mass)
    return model

pieces

pieces() -> list[list[int]]

The structure's disconnected parts, largest first.

One piece is a structure; more than one is that many free bodies, each bringing its own six zero-frequency modes. It is the first thing to ask of any model that comes back too floppy, and the answer is almost never what was intended — the old airplane fixture, meshed and wired through its own drawn lines, turned out to be three: the fuselage, a wing and the tail, none joined.

Source code in src/visualdynamics/core/fem.py

def pieces(self) -> list[list[int]]:
    """The structure's disconnected parts, largest first.

    One piece is a structure; more than one is that many free bodies,
    each bringing its own six zero-frequency modes. It is the first
    thing to ask of any model that comes back too floppy, and the
    answer is almost never what was intended — the old airplane
    fixture, meshed and wired through its own drawn lines, turned out
    to be three: the fuselage, a wing and the tail, none joined.
    """
    neighbors: dict[int, set[int]] = {n: set() for n in self.node_ids}
    for beam in self.beams:
        if beam.node_a in neighbors and beam.node_b in neighbors:
            neighbors[beam.node_a].add(beam.node_b)
            neighbors[beam.node_b].add(beam.node_a)
    for plate in self.plates:
        for k, node in enumerate(plate.nodes):
            other = plate.nodes[(k + 1) % 4]
            neighbors[node].add(other)
            neighbors[other].add(node)
    for triangle in self.triangles:
        for k, node in enumerate(triangle.nodes):
            other = triangle.nodes[(k + 1) % 3]
            neighbors[node].add(other)
            neighbors[other].add(node)
    for solid in self.solids:
        first = solid.nodes[0]
        for node in solid.nodes[1:]:
            neighbors[first].add(node)
            neighbors[node].add(first)
    for link in self.rigid_links:
        neighbors[link.node_a].add(link.node_b)
        neighbors[link.node_b].add(link.node_a)
    # a spring joins what it connects: a test article on springs to
    # its fixture is one structure, not two free bodies
    for spring in self.springs:
        if spring.node_b is not None:
            neighbors[spring.node_a].add(spring.node_b)
            neighbors[spring.node_b].add(spring.node_a)
    return connected_pieces(neighbors)

wire_faces

wire_faces(material: Material, section: Section, color: int = 1, group: str = '') -> int

Put a beam along every edge of every face, and return how many.

This is the shortest road from a shape to a structure: draw the thing as a surface mesh, and let its own edges be its members. The geometry then is the model, with nothing derived, nothing tied, and no second set of nodes that only exist to be looked at.

Beams, not axial springs. A spring on each edge leaves a quad free to shear and a flat sheet free to fold — the edges never change length, so nothing resists it — and the model comes back a mechanism with as many zero-frequency modes as it has panels. A beam carries moment, so a wireframe of them is a space frame and stands up. It is the same element the rest of this module uses.

Shared edges are wired once. An edge that already has a beam is left alone, so explicit members (a truss strut, a standoff) can be placed first and keep their own section.

Source code in src/visualdynamics/core/fem.py

def wire_faces(self, material: Material, section: Section,
               color: int = 1, group: str = '') -> int:
    """Put a beam along every edge of every face, and return how many.

    This is the shortest road from a shape to a structure: draw the
    thing as a surface mesh, and let its own edges be its members. The
    geometry then *is* the model, with nothing derived, nothing tied,
    and no second set of nodes that only exist to be looked at.

    Beams, not axial springs. A spring on each edge leaves a quad free
    to shear and a flat sheet free to fold — the edges never change
    length, so nothing resists it — and the model comes back a
    mechanism with as many zero-frequency modes as it has panels. A
    beam carries moment, so a wireframe of them is a space frame and
    stands up. It is the same element the rest of this module uses.

    Shared edges are wired once. An edge that already has a beam is
    left alone, so explicit members (a truss strut, a standoff) can be
    placed first and keep their own section.
    """
    seen = {frozenset((beam.node_a, beam.node_b)) for beam in self.beams}
    added = 0
    for face in self.faces:
        nodes = face.nodes
        for k, node in enumerate(nodes):
            other = nodes[(k + 1) % len(nodes)]
            edge = frozenset((node, other))
            if edge in seen or node == other:
                continue
            seen.add(edge)
            self.add_beam(node, other, material, section, color=color,
                          group=group or f'edge {len(self.beams)}')
            added += 1
    return added

distribute_mass

distribute_mass(total: float, spin: float = 1.0 / 12.0) -> dict[int, float]

Share total over the nodes by the members meeting at each.

A node's share is proportional to the length of member it carries — half of every member that reaches it — so mass follows the material rather than the mesh. Returns {node: mass}.

Sharing it equally is the obvious thing and it is a trap, because a node count is a statement about how finely something was drawn rather than about how much of it there is. On the demonstration airframe that put 42% of the mass into the propellers — they need the finest mesh to look right, so they collect the most nodes — and 117 g on each rotor buried every airframe mode below 1.8 kHz under blade motion, while the battery, the heaviest real item on the aircraft, was left with 51 g.

Each node also gets a rotary inertia, without which the three rotational degrees of freedom carry nothing, the mass matrix is singular and the Cholesky factorization fails outright. It is taken as spin * m * L^2 with L the mean length of the members meeting there — the patch of structure the node stands for, spun about its own middle. The default 1/12 is a uniform rod's.

Source code in src/visualdynamics/core/fem.py

def distribute_mass(self, total: float, spin: float = 1.0 / 12.0
                    ) -> dict[int, float]:
    """Share `total` over the nodes by the members meeting at each.

    A node's share is proportional to the length of member it carries —
    half of every member that reaches it — so mass follows the material
    rather than the mesh. Returns {node: mass}.

    Sharing it *equally* is the obvious thing and it is a trap, because
    a node count is a statement about how finely something was drawn
    rather than about how much of it there is. On the demonstration
    airframe that put 42% of the mass into the propellers — they need
    the finest mesh to look right, so they collect the most nodes — and
    117 g on each rotor buried every airframe mode below 1.8 kHz under
    blade motion, while the battery, the heaviest real item on the
    aircraft, was left with 51 g.

    Each node also gets a rotary inertia, without which the three
    rotational degrees of freedom carry nothing, the mass matrix is
    singular and the Cholesky factorization fails outright. It is
    taken as `spin * m * L^2` with L the mean length of the members
    meeting there — the patch of structure the node stands for, spun
    about its own middle. The default 1/12 is a uniform rod's.
    """
    nodes = self.node_ids
    if not nodes:
        raise ValueError('there are no nodes to share the mass over')
    reach: dict[int, list[float]] = {node: [] for node in nodes}
    for beam in self.beams:
        length = self._length(beam)
        for node in (beam.node_a, beam.node_b):
            if node in reach:
                reach[node].append(length)
    carried = {node: sum(lengths) / 2.0 for node, lengths in reach.items()}
    span = sum(carried.values())
    if span <= 0.0:
        raise ValueError(
            'no node carries any member, so there is nothing to share '
            'the mass out in proportion to')

    self.masses = [item for item in self.masses if item.name != 'shared']
    shares = {}
    for node in nodes:
        share = float(total) * carried[node] / span
        lengths = reach[node]
        scale = float(np.mean(lengths)) if lengths else 0.0
        inertia = spin * share * scale ** 2
        self.add_mass(node, share, (inertia, inertia, inertia),
                      name='shared')
        shares[node] = share
    return shares

group

group(node: int) -> str

What this node was added as part of — 'arm 2', 'top deck'.

A label for the modeler's own use: nothing here reads it, but working out which part of a structure a mode lives in is the first question asked of any result, and reconstructing it from coordinates afterwards is guesswork.

Source code in src/visualdynamics/core/fem.py

def group(self, node: int) -> str:
    """What this node was added as part of — 'arm 2', 'top deck'.

    A label for the modeler's own use: nothing here reads it, but
    working out which part of a structure a mode lives in is the first
    question asked of any result, and reconstructing it from
    coordinates afterwards is guesswork.
    """
    return self._node_group[int(node)]

dof_strings

dof_strings() -> list[str]

'101X+', '101Y+', … in the matrices' own order.

Source code in src/visualdynamics/core/fem.py

def dof_strings(self) -> list[str]:
    """'101X+', '101Y+', … in the matrices' own order."""
    return [f'{node}{direction}'
            for node in self._nodes for direction in DIRECTIONS]

matrices

matrices() -> tuple[ndarray, ndarray]

Assemble the global mass and stiffness matrices, dense.

Rows and columns run structural node by structural node in insertion order, six per node in DIRECTIONS order, which is what dof_strings() spells out. Display nodes are absent: they carry nothing, so there is nothing of theirs to assemble.

Source code in src/visualdynamics/core/fem.py

def matrices(self) -> tuple[np.ndarray, np.ndarray]:
    """Assemble the global mass and stiffness matrices, dense.

    Rows and columns run structural node by structural node in
    insertion order, six per node in `DIRECTIONS` order, which is what
    `dof_strings()` spells out. Display nodes are absent: they carry
    nothing, so there is nothing of theirs to assemble.
    """
    n = self.num_dof
    mass = np.zeros((n, n), dtype=np.float64)
    stiffness = np.zeros((n, n), dtype=np.float64)
    for rows, k, m in self._contributions():
        grid = np.ix_(rows, rows)
        if k is not None:
            stiffness[grid] += k
        mass[grid] += m
    # assembly is symmetric by construction, but floating point addition
    # is not associative and the halves drift apart in the last bits;
    # eigh reads only one triangle, so an asymmetry here is silent
    return (mass + mass.T) / 2.0, (stiffness + stiffness.T) / 2.0

sparse_matrices

sparse_matrices(ticker: Any = None) -> tuple[Any, Any]

The same mass and stiffness matrices as matrices, stored sparse (scipy CSR): each node couples only to the nodes of the elements it touches, so a row holds a few dozen entries of thousands, and the storage grows with the nodes rather than their square.

Returns:

Type	Description
tuple of scipy.sparse.csr_matrix	(mass, stiffness), symmetric.

Source code in src/visualdynamics/core/fem.py

def sparse_matrices(self, ticker: Any = None) -> tuple[Any, Any]:
    """The same mass and stiffness matrices as `matrices`, stored
    sparse (scipy CSR): each node couples only to the nodes of the
    elements it touches, so a row holds a few dozen entries of
    thousands, and the storage grows with the nodes rather than
    their square.

    Returns
    -------
    tuple of scipy.sparse.csr_matrix
        (mass, stiffness), symmetric.
    """
    from scipy import sparse

    n = self.num_dof
    rows_k, cols_k, vals_k, rows_m, cols_m, vals_m = [], [], [], [], [], []
    if ticker is not None:
        ticker.add(len(self.beams) + len(self.plates)
                   + len(self.triangles) + len(self.solids)
                   + len(self.springs))
    for rows, k, m in self._batches():
        count, r = rows.shape
        i = np.broadcast_to(rows[:, :, None], (count, r, r)).ravel()
        j = np.broadcast_to(rows[:, None, :], (count, r, r)).ravel()
        if k is not None:
            rows_k.append(i)
            cols_k.append(j)
            vals_k.append(k.ravel())
            if ticker is not None:
                ticker.tick(count)
        rows_m.append(i)
        cols_m.append(j)
        vals_m.append(m.ravel())

    def gather(rows, cols, vals):
        if not rows:
            return sparse.csr_matrix((n, n))
        matrix = sparse.coo_matrix(
            (np.concatenate(vals), (np.concatenate(rows),
                                    np.concatenate(cols))),
            shape=(n, n)).tocsr()
        return ((matrix + matrix.T) / 2.0).tocsr()

    return gather(rows_m, cols_m, vals_m), gather(rows_k, cols_k, vals_k)

rigid_body_vectors

rigid_body_vectors() -> ndarray

The six rigid-body motions, as columns over the model's DOFs.

Written down from the node positions rather than found from the matrices: they are what the null space of an unconstrained stiffness matrix is, and knowing them in advance is what lets a rigid mode be recognized as rigid rather than as a very soft one.

Source code in src/visualdynamics/core/fem.py

def rigid_body_vectors(self) -> np.ndarray:
    """The six rigid-body motions, as columns over the model's DOFs.

    Written down from the node positions rather than found from the
    matrices: they are what the null space of an unconstrained
    stiffness matrix *is*, and knowing them in advance is what lets a
    rigid mode be recognized as rigid rather than as a very soft one.
    """
    n = self.num_dof
    solved = self.node_ids
    vectors = np.zeros((n, 6), dtype=np.float64)
    center = np.mean(np.array([self._nodes[k] for k in solved]), axis=0)
    for i, node in enumerate(solved):
        offset = self._nodes[node] - center
        for axis in range(3):
            vectors[6 * i + axis, axis] = 1.0            # translations
            # a small rotation about `axis` moves a point by the cross
            # product of the axis with its offset, and rotates it by
            # the axis itself
            spin = np.zeros(3)
            spin[axis] = 1.0
            vectors[6 * i:6 * i + 3, 3 + axis] = np.cross(spin, offset)
            vectors[6 * i + 3 + axis, 3 + axis] = 1.0
    return vectors

eigensolution

eigensolution(maximum_frequency: float | None = None, num_modes: int | None = None, damping: float = 0.0, fixed: Sequence[str] = (), solver: str = 'auto', progress: Callable[[int, int], None] | None = None) -> ShapeSet

Real normal modes, mass-normalized, as a ShapeSet.

fixed names degrees of freedom to ground: '101X+' fixes one, '101' fixes all six of that node. Rigid links are eliminated exactly (constraint_transform).

Two solvers, one answer. Dense (a model of up to SPARSE_ABOVE degrees of freedom): the symmetric generalized problem K phi = lambda M phi solved whole, by factoring M (Cholesky), reducing to a standard symmetric problem and transforming back — every mode, mass-normalized to machine precision. Sparse (larger models): the matrices stored as their nonzeros (sparse_matrices) and the lowest modes found by shift-invert Lanczos (ARPACK, scipy.sparse.linalg.eigsh), the family the large finite element codes use; it finds the lowest num_modes, or every mode up to maximum_frequency, and one of the two must be said. Memory grows with the nodes instead of their square: a plate model of 1,500 nodes, ~13 GB dense, is tens of megabytes sparse (Brandon, 2026-09-26, for the BARC example — the dense ceiling, deliberate until then, measured and met).

damping is a fraction of critical, applied uniformly. A model has no damping of its own; it is stated so the modes can synthesize an FRF that looks like a measurement.

Parameters:

Name	Type	Description	Default
maximum_frequency	float	Keep every mode up to this frequency, in Hz.	None
num_modes	int	Keep this many, the lowest, rigid ones included.	None
damping	float	The fraction of critical damping every mode is given.	0.0
fixed	sequence of str	Degrees of freedom to ground.	()
solver	('auto', 'dense', 'sparse')	Which solver; 'auto' is dense up to SPARSE_ABOVE degrees of freedom and sparse beyond.	'auto'
progress	callable	Told (done, total) as the solve advances: the assembly per element, the factorization, the eigen iterations (an estimated length, extended while they run), the polish. The window's strip bar reads it; anything it raises stops the solve.	None

Returns:

Type	Description
ShapeSet

Source code in src/visualdynamics/core/fem.py

def eigensolution(self, maximum_frequency: float | None = None,
                  num_modes: int | None = None, damping: float = 0.0,
                  fixed: Sequence[str] = (),
                  solver: str = 'auto',
                  progress: Callable[[int, int], None] | None = None
                  ) -> ShapeSet:
    """Real normal modes, mass-normalized, as a ShapeSet.

    `fixed` names degrees of freedom to ground: '101X+' fixes one,
    '101' fixes all six of that node. Rigid links are eliminated
    exactly (`constraint_transform`).

    Two solvers, one answer. **Dense** (a model of up to
    `SPARSE_ABOVE` degrees of freedom): the symmetric generalized
    problem K phi = lambda M phi solved whole, by factoring M
    (Cholesky), reducing to a standard symmetric problem and
    transforming back — every mode, mass-normalized to machine
    precision. **Sparse** (larger models): the matrices stored as
    their nonzeros (`sparse_matrices`) and the lowest modes found by
    shift-invert Lanczos (ARPACK, `scipy.sparse.linalg.eigsh`), the
    family the large finite element codes use; it finds the lowest
    `num_modes`, or every mode up to `maximum_frequency`, and one of
    the two must be said. Memory grows with the nodes instead of
    their square: a plate model of 1,500 nodes, ~13 GB dense, is tens
    of megabytes sparse (Brandon, 2026-09-26, for the BARC example —
    the dense ceiling, deliberate until then, measured and met).

    `damping` is a fraction of critical, applied uniformly. A model
    has no damping of its own; it is stated so the modes can
    synthesize an FRF that looks like a measurement.

    Parameters
    ----------
    maximum_frequency : float, optional
        Keep every mode up to this frequency, in Hz.
    num_modes : int, optional
        Keep this many, the lowest, rigid ones included.
    damping : float, default 0.0
        The fraction of critical damping every mode is given.
    fixed : sequence of str
        Degrees of freedom to ground.
    solver : {'auto', 'dense', 'sparse'}, default 'auto'
        Which solver; 'auto' is dense up to `SPARSE_ABOVE` degrees
        of freedom and sparse beyond.
    progress : callable, optional
        Told ``(done, total)`` as the solve advances: the assembly
        per element, the factorization, the eigen iterations (an
        estimated length, extended while they run), the polish.
        The window's strip bar reads it; anything it raises stops
        the solve.

    Returns
    -------
    ShapeSet
    """
    if solver not in ('auto', 'dense', 'sparse'):
        raise ValueError(f'{solver!r} is not a solver: auto, dense, sparse')
    ticker = Ticker(progress)
    if solver == 'sparse' or (solver == 'auto'
                              and self.num_dof > SPARSE_ABOVE):
        eigenvalues, full, stiffness = self._sparse_modes(
            fixed, maximum_frequency, num_modes, ticker)
    else:
        eigenvalues, full, stiffness = self._dense_modes(fixed, ticker)
    if ticker.total:
        ticker.tick(ticker.total - ticker.done)
    # a rigid-body eigenvalue is zero plus round-off, and comes out
    # either side of it; the negative ones are not oscillations
    frequency = np.sqrt(np.clip(eigenvalues, 0.0, None)) / (2.0 * np.pi)
    frequency[_strains_nothing(full, stiffness)] = 0.0

    order = np.argsort(frequency, kind='stable')
    frequency, full = frequency[order], full[:, order]
    keep = np.ones(len(frequency), dtype=bool)
    if maximum_frequency is not None:
        keep &= frequency <= float(maximum_frequency)
    keep = np.flatnonzero(keep)
    if num_modes is not None:
        keep = keep[:int(num_modes)]

    return ShapeSet(frequency=frequency[keep],
                    damping=np.full(len(keep), float(damping)),
                    coordinate=self.dof_strings(),
                    shape_matrix=full[:, keep].T,
                    mass_unit='kg',
                    comment=self.name or '')

scaled_system

scaled_system(fixed: Sequence[str] = (), ticker: Any = None) -> tuple[Any, ...]

The sparse eigenproblem as the sparse solver poses it.

Returns (K, M, T, S, sigma, stiffness): the constrained, symmetrically scaled stiffness and mass in CSC form, the constraint transform T and the scaling S that carry a solution back to every degree of freedom as T (S phi), the shift sigma, and the unconstrained sparse stiffness. Public so a test can hand the polish a start of its own choosing.

The shift is just below zero: the rigid-body modes (eigenvalue 0) are nearest it, so they come first, and K - sigma M = K + |sigma| M is positive definite even when K alone is singular (free-free). The scaling is symmetric and diagonal, S K S and S M S with S = 1/sqrt of the shifted diagonal: the same eigenvalues, the modes S phi'. A model of plates mixes meters with radians, and rigid links fold lever arms into the rotations — the rigid-link test model's shifted matrix had a condition number of 4e12, and Lanczos, which converges no better than its linear solves, left residuals of 1e-2. Scaled it is 4e9 (2026-09-26).

Source code in src/visualdynamics/core/fem.py

def scaled_system(self, fixed: Sequence[str] = (),
                  ticker: Any = None) -> tuple[Any, ...]:
    """The sparse eigenproblem as the sparse solver poses it.

    Returns (K, M, T, S, sigma, stiffness): the constrained,
    symmetrically scaled stiffness and mass in CSC form, the
    constraint transform T and the scaling S that carry a solution
    back to every degree of freedom as T (S phi), the shift sigma, and
    the unconstrained sparse stiffness. Public so a test can hand the
    polish a start of its own choosing.

    The shift is just below zero: the rigid-body modes (eigenvalue 0)
    are nearest it, so they come first, and K - sigma M = K + |sigma| M
    is positive definite even when K alone is singular (free-free).
    The scaling is symmetric and diagonal, S K S and S M S with S =
    1/sqrt of the shifted diagonal: the same eigenvalues, the modes
    S phi'. A model of plates mixes meters with radians, and rigid
    links fold lever arms into the rotations — the rigid-link test
    model's shifted matrix had a condition number of 4e12, and
    Lanczos, which converges no better than its linear solves, left
    residuals of 1e-2. Scaled it is 4e9 (2026-09-26).
    """
    from scipy import sparse

    mass, stiffness = self.sparse_matrices(ticker)
    transform = self.constraint_transform(fixed, sparse=True)
    reduced_m = (transform.T @ mass @ transform).tocsc()
    reduced_k = (transform.T @ stiffness @ transform).tocsc()
    if not reduced_m.shape[0]:
        raise ValueError('every degree of freedom is fixed')
    sigma = -(2.0 * np.pi) ** 2
    scale = sparse.diags(1.0 / np.sqrt(
        (reduced_k - sigma * reduced_m).diagonal()))
    reduced_k = (scale @ reduced_k @ scale).tocsc()
    reduced_m = (scale @ reduced_m @ scale).tocsc()
    return reduced_k, reduced_m, transform, scale, sigma, stiffness

constraint_transform

constraint_transform(fixed: Sequence[str] = (), sparse: bool = False) -> Any

T, with u = T q: every degree of freedom of the model written in terms of those that remain free — grounded ones gone, and each rigid body's followers written through its first node, u_b = u_a + θ_a × (x_b − x_a) and θ_b = θ_a. The identity's columns when there is nothing to constrain.

Parameters:

Name	Type	Description	Default
fixed	sequence of str	Degrees of freedom to ground, as eigensolution takes them. A rigid body is grounded through its first node; naming one of its followers is refused, since fixing a follower alone would fix part of a rigid body and not the rest.	()
sparse	bool	Return it as a scipy CSR matrix, for the sparse solver.	False

Returns:

Type	Description
ndarray or csr_matrix	(num_dof, remaining) and real.

Source code in src/visualdynamics/core/fem.py

def constraint_transform(self, fixed: Sequence[str] = (),
                         sparse: bool = False) -> Any:
    """T, with u = T q: every degree of freedom of the model written
    in terms of those that remain free — grounded ones gone, and each
    rigid body's followers written through its first node, u_b = u_a +
    θ_a × (x_b − x_a) and θ_b = θ_a. The identity's columns when there
    is nothing to constrain.

    Parameters
    ----------
    fixed : sequence of str
        Degrees of freedom to ground, as `eigensolution` takes them. A
        rigid body is grounded through its first node; naming one of
        its followers is refused, since fixing a follower alone would
        fix part of a rigid body and not the rest.

    sparse : bool, default False
        Return it as a scipy CSR matrix, for the sparse solver.

    Returns
    -------
    numpy.ndarray or scipy.sparse.csr_matrix
        (num_dof, remaining) and real.
    """
    index = {node: 6 * i for i, node in enumerate(self.node_ids)}
    follows: dict[int, int] = {}
    for body in self.rigid_bodies():
        for node in body[1:]:
            follows[node] = body[0]
    free = self._free_dofs(fixed)
    followers = {index[node] + k for node in follows for k in range(6)}
    held = sorted(set(range(self.num_dof)) - {int(i) for i in free})
    clash = sorted({self.node_ids[i // 6] for i in held if i in followers})
    if clash:
        raise ValueError(
            'node ' + ', '.join(str(n) for n in clash) + ' follows a '
            'rigid link; ground the node it follows instead')
    kept = [int(i) for i in free if int(i) not in followers]
    column = {dof: k for k, dof in enumerate(kept)}
    entries = [(dof, k, 1.0) for dof, k in column.items()]
    for node, lead in follows.items():
        r = self._nodes[node] - self._nodes[lead]
        # θ × r = -[r]× θ: the lead's rotation moves the follower
        arm = np.array([[0.0, r[2], -r[1]],
                        [-r[2], 0.0, r[0]],
                        [r[1], -r[0], 0.0]])
        rows = index[node]
        lead_row = index[lead]
        for k in range(6):
            source = lead_row + k
            if source not in column:
                continue                  # the lead is grounded there
            entries.append((rows + k, column[source], 1.0))
            if k >= 3:
                entries.extend((rows + axis, column[source],
                                float(arm[axis, k - 3]))
                               for axis in range(3) if arm[axis, k - 3])
    shape = (self.num_dof, len(kept))
    if sparse:
        from scipy import sparse as sp

        if not entries:
            return sp.csr_matrix(shape)
        i, j, v = zip(*entries, strict=True)
        return sp.coo_matrix((v, (i, j)), shape=shape).tocsr()
    transform = np.zeros(shape, dtype=np.float64)
    for i, j, v in entries:
        transform[i, j] = v
    return transform

dangling_rotations

dangling_rotations() -> list[int]

The nodes whose rotations nothing acts on: touched by solids or translational springs and by nothing that carries a rotation — no beam, plate or triangle, and no rigid link, whose lead's rotation moves its followers. Neither a solid nor a spring along an axis gives a rotation stiffness, so these rotations are grounded by the eigensolution; left free they would be degrees of freedom with neither mass nor stiffness, which no factorization survives. A mass on springs is the newer case (2026-10-08): a tuned mass hung from a structure by spring lines, its node no beam's.

Returns:

Type	Description
list of int	The nodes, in the model's order.

Source code in src/visualdynamics/core/fem.py

def dangling_rotations(self) -> list[int]:
    """The nodes whose rotations nothing acts on: touched by solids or
    translational springs and by nothing that carries a rotation — no
    beam, plate or triangle, and no rigid link, whose lead's rotation
    moves its followers. Neither a solid nor a spring along an axis
    gives a rotation stiffness, so these rotations are grounded by the
    eigensolution; left free they would be degrees of freedom with
    neither mass nor stiffness, which no factorization survives. A
    mass on springs is the newer case (2026-10-08): a tuned mass hung
    from a structure by spring lines, its node no beam's.

    Returns
    -------
    list of int
        The nodes, in the model's order.
    """
    rotating = self._rotating_nodes()
    touched = {n for solid in self.solids for n in solid.nodes}
    touched |= {node for spring in self.springs
                for node, code in ((spring.node_a, spring.direction_a),
                                   (spring.node_b, spring.direction_b))
                if node is not None and abs(code) <= 3}
    # a held node's rotations are held already
    held = set(self.loose_nodes()) | {
        node for node, axes in self.grounds.items()
        if {3, 4, 5} <= set(axes)}
    return [n for n in self.node_ids
            if n in touched and n not in rotating and n not in held]

idle_lead_rotations

idle_lead_rotations() -> list[tuple[int, int]]

(lead node, axis 0-2) for each rotation of a rigid body's lead that nothing acts on: a body only solids touch, whose nodes lie on one line, cannot feel a rotation about that line — it moves no follower, and a solid has no rotational stiffness — so that rotation is a degree of freedom with neither mass nor stiffness.

Found tying a beam of bricks to the channels under it (the BARC in hexes, 2026-10-07): the two meshes' columns line up, so each tied channel node leads only the beam nodes straight above it, and the app's Tie on any two conforming brick meshes does the same. Grounded like the rotations of a node only solids touch (dangling_rotations): the rotation about the line where the line runs along a global axis, all three where the body's nodes coincide. A line at a slant is left to the factorization's own refusal, which names it; grounding an oblique axis would need a rotated basis, and no model here has needed one.

Returns:

Type	Description
list of (int, int)	The lead's node and the rotation's axis, in the model's order.

Source code in src/visualdynamics/core/fem.py

def idle_lead_rotations(self) -> list[tuple[int, int]]:
    """(lead node, axis 0-2) for each rotation of a rigid body's lead
    that nothing acts on: a body only solids touch, whose nodes lie
    on one line, cannot feel a rotation about that line — it moves
    no follower, and a solid has no rotational stiffness — so that
    rotation is a degree of freedom with neither mass nor stiffness.

    Found tying a beam of bricks to the channels under it (the BARC
    in hexes, 2026-10-07): the two meshes' columns line up, so each
    tied channel node leads only the beam nodes straight above it,
    and the app's Tie on any two conforming brick meshes does the
    same. Grounded like the rotations of a node only solids touch
    (`dangling_rotations`): the rotation about the line where the
    line runs along a global axis, all three where the body's nodes
    coincide. A line at a slant is left to the factorization's own
    refusal, which names it; grounding an oblique axis would need a
    rotated basis, and no model here has needed one.

    Returns
    -------
    list of (int, int)
        The lead's node and the rotation's axis, in the model's
        order.
    """
    if not self.solids:
        return []
    carried = {n for beam in self.beams for n in (beam.node_a, beam.node_b)}
    carried |= {n for plate in self.plates for n in plate.nodes}
    carried |= {n for triangle in self.triangles for n in triangle.nodes}
    carried |= {m.node for m in self.masses if any(m.inertia)}
    idle = []
    for body in self.rigid_bodies():
        if any(node in carried for node in body):
            continue
        lead = body[0]
        arms = np.array([self._nodes[node] - self._nodes[lead]
                         for node in body[1:]], dtype=float)
        size = float(np.abs(arms).max()) if arms.size else 0.0
        if size == 0.0:
            idle.extend((lead, axis) for axis in range(3))
            continue
        _u, spread, direction = np.linalg.svd(arms / size)
        if spread.size > 1 and spread[1] > 1e-9 * spread[0]:
            continue                      # not on one line: all stiff
        axis = int(np.argmax(np.abs(direction[0])))
        if abs(direction[0][axis]) > 1.0 - 1e-9:
            idle.append((lead, axis))
    return idle

loose_nodes

loose_nodes() -> list[int]

The nodes nothing touches: no element, no rigid link, no lumped mass. A finite element deck carries them routinely — a reference point, a constraint's own grid — and a model built from its element groups grounds them whole rather than refusing the deck (2026-09-30); the grillage path still refuses, since there a loose node is a drawing that was never wired.

Returns:

Type	Description
list of int	The nodes, in the model's order.

Source code in src/visualdynamics/core/fem.py

def loose_nodes(self) -> list[int]:
    """The nodes nothing touches: no element, no rigid link, no
    lumped mass. A finite element deck carries them routinely — a
    reference point, a constraint's own grid — and a model built
    from its element groups grounds them whole rather than refusing the
    deck (2026-09-30); the grillage path still refuses, since there
    a loose node is a drawing that was never wired.

    Returns
    -------
    list of int
        The nodes, in the model's order.
    """
    touched = self._rotating_nodes()
    touched |= {n for solid in self.solids for n in solid.nodes}
    touched |= {m.node for m in self.masses}
    return [n for n in self.node_ids if n not in touched]

geometry

geometry(beams: bool = True) -> Geometry

The model as a Geometry: nodes, beams as elements, faces as faces.

Beams become element type 21 (beam2) rather than tracelines, because they are elements — a traceline is a line drawn through nodes to make a display readable, and confusing the two would make the model's own connectivity indistinguishable from a drawing aid the moment anything edited it.

beams=False leaves them out, for a model whose members are all wrapped in surfaces: there the beams run inside the shells, so drawing them puts a wireframe over the thing it is the skeleton of.

Source code in src/visualdynamics/core/fem.py

def geometry(self, beams: bool = True) -> Geometry:
    """The model as a Geometry: nodes, beams as elements, faces as faces.

    Beams become element type 21 (beam2) rather than tracelines,
    because they *are* elements — a traceline is a line drawn through
    nodes to make a display readable, and confusing the two would make
    the model's own connectivity indistinguishable from a drawing aid
    the moment anything edited it.

    `beams=False` leaves them out, for a model whose members are all
    wrapped in surfaces: there the beams run *inside* the shells, so
    drawing them puts a wireframe over the thing it is the skeleton of.
    """
    node_ids = self.node_ids
    connectivity, types, colors = [], [], []
    for beam in self.beams if beams else ():
        connectivity.append([beam.node_a, beam.node_b])
        types.append(21)
        colors.append(beam.color)
    # rigid links are two-node lines too, in element groups of their own — a
    # geometry given `RIGID` on those element groups rebuilds them
    for link in self.rigid_links:
        connectivity.append([link.node_a, link.node_b])
        types.append(21)
        colors.append(1)
    # plates are structural quads and appear as such — unlike faces,
    # which are drawings; both shade, only one carries stiffness
    for plate in self.plates:
        connectivity.append(list(plate.nodes))
        types.append(44)
        colors.append(plate.color)
    for triangle in self.triangles:
        connectivity.append(list(triangle.nodes))
        types.append(41)
        colors.append(triangle.color)
    for solid in self.solids:
        connectivity.append(list(solid.nodes))
        types.append(SOLID_CODES[len(solid.nodes)])
        colors.append(solid.color)
    for face in self.faces:
        connectivity.append(list(face.nodes))
        types.append(44 if len(face.nodes) == 4 else 41)
        colors.append(face.color)
    # Elements go into the element group of the part they belong to, so the
    # geometry states its own regions and a saved file can rebuild the
    # structure. Without this the element groups live only on the drawing, and
    # the model's own geometry — which is what gets saved — carries
    # nothing: the file comes back with every member the same section.
    # The part comes off the element, never off its first node: a node
    # on a seam belongs to two parts and answers for one, which put
    # five of the drone's blade members into the frame element group.
    parts, groups = {}, []
    for part in ([b.group for b in (self.beams if beams else ())]
                 + [link.group or 'rigid links'
                    for link in self.rigid_links]
                 + [p.group for p in self.plates]
                 + [t.group for t in self.triangles]
                 + [s.group for s in self.solids]
                 + [f.group for f in self.faces]):
        groups.append(parts.setdefault(part or 'body', len(parts) + 1))
    return Geometry(
        node_id=node_ids,
        node_xyz=np.array([self._nodes[node] for node in node_ids]),
        elem_conn=connectivity or None,
        elem_type=types or None,
        elem_color=colors or None,
        elem_group=groups or None,
        group_id=list(parts.values()) or None,
        group_name=list(parts) or None,
        length_unit=self.length_unit)

Functions:

ground_axes

ground_axes(directions: Any) -> tuple[str, ...]

A ground's directions as AXES names, in their order: True for all six, False or None or nothing for none, else names ('X', 'RY'), a sign or case making no difference — 'x+' is 'X'.

Source code in src/visualdynamics/core/fem.py

def ground_axes(directions: Any) -> tuple[str, ...]:
    """A ground's directions as `AXES` names, in their order: True for
    all six, False or None or nothing for none, else names ('X', 'RY'),
    a sign or case making no difference — 'x+' is 'X'."""
    if directions is True:
        return AXES
    if not directions:
        return ()
    if isinstance(directions, str):
        directions = directions.replace(',', ' ').split()
    wanted = set()
    for name in directions:
        axis = str(name).strip().upper().rstrip('+-')
        if axis not in AXES:
            raise ValueError(f'{name!r} is not a direction to hold: '
                             + ', '.join(AXES))
        wanted.add(axis)
    return tuple(axis for axis in AXES if axis in wanted)

material

material(name: str) -> Material

A library material by name — material('6061-T6') — the same entry the Element Groups table's Material drop-down fills a row from.

Raises KeyError naming the library when the name is not in it, so a typo reads as one rather than as a missing material.

Source code in src/visualdynamics/core/fem.py

def material(name: str) -> Material:
    """A library material by name — `material('6061-T6')` — the same
    entry the Element Groups table's Material drop-down fills a row from.

    Raises `KeyError` naming the library when the name is not in it,
    so a typo reads as one rather than as a missing material.
    """
    try:
        return MATERIALS[name]
    except KeyError:
        raise KeyError(f'{name!r} is not in the material library; it '
                       f'has {", ".join(MATERIALS)}') from None

with_article

with_article(shape: str) -> str

'a channel', 'an I-beam', 'an angle' — a shape named in a sentence.

Source code in src/visualdynamics/core/fem.py

def with_article(shape: str) -> str:
    """'a channel', 'an I-beam', 'an angle' — a shape named in a
    sentence."""
    return ('an ' if shape[:1].lower() in 'aeio' or shape.startswith('I-')
            else 'a ') + shape

angle_major_axis

angle_major_axis(long_leg: float, short_leg: float, thickness: float) -> float

Where an angle's major principal axis lies — the direction its orientation vector points — in degrees from the long leg, turning toward the short leg. 45 for an equal angle.

Parameters:

Name	Type	Description	Default
long_leg	float	The angle's dimensions, in any one unit.	required
short_leg	float	The angle's dimensions, in any one unit.	required
thickness	float	The angle's dimensions, in any one unit.	required

Returns:

Type	Description
float	Degrees from the long leg, 0 to 180.

Source code in src/visualdynamics/core/fem.py

def angle_major_axis(long_leg: float, short_leg: float,
                     thickness: float) -> float:
    """Where an angle's major principal axis lies — the direction its
    orientation vector points — in degrees from the long leg, turning
    toward the short leg. 45 for an equal angle.

    Parameters
    ----------
    long_leg, short_leg, thickness : float
        The angle's dimensions, in any one unit.

    Returns
    -------
    float
        Degrees from the long leg, 0 to 180.
    """
    return _angle_properties(long_leg, short_leg, thickness)[2]

connected_pieces

connected_pieces(neighbors: dict[int, set[int]]) -> list[list[int]]

The connected components of an adjacency map, largest first.

One piece is a structure; more than one is that many free bodies, each bringing its own six zero-frequency modes. The flood fill is shared by the model (asking over its beams and plates) and the drone's drawing (asking over its face edges, before a model exists), so the two cannot disagree about what "joined" means.

Source code in src/visualdynamics/core/fem.py

def connected_pieces(neighbors: dict[int, set[int]]) -> list[list[int]]:
    """The connected components of an adjacency map, largest first.

    One piece is a structure; more than one is that many free bodies,
    each bringing its own six zero-frequency modes. The flood fill is
    shared by the model (asking over its beams and plates) and the
    drone's drawing (asking over its face edges, before a model
    exists), so the two cannot disagree about what "joined" means.
    """
    seen: set[int] = set()
    found: list[list[int]] = []
    for start in neighbors:
        if start in seen:
            continue
        stack, piece = [start], []
        while stack:
            node = stack.pop()
            if node in seen:
                continue
            seen.add(node)
            piece.append(node)
            stack.extend(neighbors[node] - seen)
        found.append(sorted(piece))
    return sorted(found, key=len, reverse=True)

polish_modes

polish_modes(stiffness: Any, mass: Any, sigma: float, vectors: ndarray, wanted: int, factor: Any = None, ticker: Any = None) -> tuple[ndarray, ndarray, ndarray]

Refine a block of approximate modes until they are converged.

Block inverse iteration with the shifted factor (K - sigma M), each step followed by Rayleigh-Ritz on the block — the small projected generalized problem solved densely — repeated until every one of the first wanted modes has a relative residual below POLISH_RESIDUAL, or a step no longer halves the worst of them (the floor the conditioning sets), or POLISH_STEPS are spent. The residual of a mode is |K v - lambda M v| over |(K - sigma M) v|, which is defined for a rigid-body mode too (its numerator and |K v| are both zero). The block is polished whole: the Ritz step resolves the modes inside it exactly, and iteration removes what leaks in from outside, at a rate set by the gap between the last wanted mode and the first beyond the block, which is why the caller hands over a guard band.

Returns the eigenvalues, the mass-normal vectors, both for the whole block in ascending order, and the residual of each after the last step. Sparse stiffness and mass in CSC form.

Source code in src/visualdynamics/core/fem.py

def polish_modes(stiffness: Any, mass: Any, sigma: float, vectors: np.ndarray,
                 wanted: int, factor: Any = None, ticker: Any = None
                 ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Refine a block of approximate modes until they are converged.

    Block inverse iteration with the shifted factor (K - sigma M), each
    step followed by Rayleigh-Ritz on the block — the small projected
    generalized problem solved densely — repeated until every one of the
    first `wanted` modes has a relative residual below `POLISH_RESIDUAL`,
    or a step no longer halves the worst of them (the floor the
    conditioning sets), or `POLISH_STEPS` are spent. The residual of a mode is
    |K v - lambda M v| over |(K - sigma M) v|, which is defined for a
    rigid-body mode too (its numerator and |K v| are both zero). The
    block is polished whole: the Ritz step resolves the modes *inside*
    it exactly, and iteration removes what leaks in from *outside*, at
    a rate set by the gap between the last wanted mode and the first
    beyond the block, which is why the caller hands over a guard band.

    Returns the eigenvalues, the mass-normal vectors, both for the whole
    block in ascending order, and the residual of each after the last
    step. Sparse `stiffness` and `mass` in CSC form.
    """
    from scipy.linalg import eigh as scipy_eigh
    from scipy.sparse.linalg import splu

    if factor is None:
        factor = splu((stiffness - sigma * mass).tocsc())
    wanted = min(int(wanted), vectors.shape[1])
    worst = math.inf
    if ticker is not None:
        ticker.add(POLISH_STEPS)
    for step in range(POLISH_STEPS):
        if ticker is not None:
            ticker.tick()
        if step:
            vectors = factor.solve(mass @ vectors)
        eigenvalues, ritz = scipy_eigh(vectors.T @ (stiffness @ vectors),
                                       vectors.T @ (mass @ vectors))
        vectors = vectors @ ritz
        k_v = stiffness @ vectors
        m_v = mass @ vectors
        residuals = (np.linalg.norm(k_v - m_v * eigenvalues, axis=0)
                     / np.linalg.norm(k_v - sigma * m_v, axis=0))
        before, worst = worst, float(residuals[:wanted].max())
        if worst < POLISH_RESIDUAL or worst > 0.5 * before:
            break
    return eigenvalues, vectors, residuals