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 modelling capability that produces something worth analysing: three-dimensional two-node beams 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 two ways in. Build a structure member by member, as the example below does; or draw a shape as a surface mesh and let from_geometry put a member along every edge of it and share the mass over its nodes. The second is the shorter road from any geometry — imported or built — to a set of modes, and visualdynamics.demo.drone is built that way throughout.

What it is not: a general finite element code. There are no shells, no solids, no constraints beyond fixing degrees of freedom, and no static solution. A plate is modelled the way a frame is, as a grillage of beams, which is a real modelling choice with a known cost rather than an approximation hidden inside an element. Wiring a surface mesh's edges is the same bargain: 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 millimetres 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

aluminium = fem.Material('aluminium', 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), aluminium, 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.

Classes:

Name	Description
Material	An isotropic elastic material.
Section	A beam cross section, as the four numbers the element needs.
Beam	One two-node beam element.
Plate	One four-node rectangular plate-bending element.
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
connected_pieces	The connected components of an adjacency map, largest first.

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 fibre 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.

Section dataclass

Section(
    name: str, area: float, iy: float, iz: float, j: 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.

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
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.
square_tube	A square tube, given its outside width and wall thickness.

Methods:

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:
        raise ValueError(f'{name}: wall {wall} is thicker than 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)

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."""
    return cls.round_tube(name, diameter, diameter / 2.0)

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."""
    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
    a, b = max(width, height) / 2.0, min(width, height) / 2.0
    # St Venant's constant for a solid rectangle, to better than 0.1%
    # over the whole aspect range (Roark): the series in b/a truncated
    # where the next term is under a part in a thousand
    j = a * b ** 3 * (16 / 3 - 3.36 * (b / a) * (1 - b ** 4 / (12 * a ** 4)))
    return cls(name, area, iy, iz, j)

square_tube classmethod

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

A square tube, given its outside width and wall thickness.

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."""
    inner = width - 2.0 * wall
    if inner <= 0:
        raise ValueError(f'{name}: wall {wall} closes the section')
    area = width ** 2 - inner ** 2
    i = (width ** 4 - inner ** 4) / 12.0
    # Bredt's thin-wall formula: J = 4 A_m^2 t / s, with A_m the area
    # inside the wall centreline and s that centreline's length
    mean = width - wall
    j = 4.0 * (mean ** 2) ** 2 * wall / (4.0 * mean)
    return cls(name, area, i, i, j)

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.

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.

Assembly is dense. The matrices are (6 x nodes) square, so a three-hundred-node model is a 1800 x 1800 pair — 26 MB and a couple of seconds to solve — and a thousand nodes is 290 MB and a minute or two. That is the working ceiling, and it is a deliberate one: a sparse assembly and a subspace solver would raise it by an order of magnitude and cost more code than the models this is here to build need. Measure before assuming it is the problem.

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. 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.
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.
rigid_body_vectors	The six rigid-body motions, as columns over the model's DOFs.
eigensolution	Real normal modes, mass-normalized, as a ShapeSet.
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.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

from_geometry classmethod

from_geometry(
    geometry: Geometry,
    material: Material,
    section: Section,
    total_mass: float | None = None,
    tracelines: bool = False,
    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.

tracelines wires the display polylines too, which is what a wireframe geometry with no elements 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,
                  section: Section,
                  total_mass: float | None = None,
                  tracelines: bool = False, 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.

    `tracelines` wires the display polylines too, which is what a
    wireframe geometry with no elements 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')
    labels = dict(groups or {})
    if not labels:
        # A geometry that carries element blocks 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_blocks(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 block 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 block and is never ambiguous.
    parts = _block_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)
    if tracelines:
        runs.extend(([int(n) for n in conn], '')
                    for conn in geometry.traceline_conn)

    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 or tracelines 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 and tracelines 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 tracelines, 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 tracelines, turned out
    to be three: the fuselage, a wing and the tail, none joined.
    """
    neighbours: dict[int, set[int]] = {n: set() for n in self.node_ids}
    for beam in self.beams:
        if beam.node_a in neighbours and beam.node_b in neighbours:
            neighbours[beam.node_a].add(beam.node_b)
            neighbours[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]
            neighbours[node].add(other)
            neighbours[other].add(node)
    return connected_pieces(neighbours)

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 modeller'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 modeller'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.

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.

    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)
    index = {node: 6 * i for i, node in enumerate(self.node_ids)}

    for beam in self.beams:
        length = self._length(beam)
        if length == 0.0:
            raise ValueError(f'beam {beam.node_a}-{beam.node_b} has zero length')
        rotation = _element_axes(self._nodes[beam.node_a],
                                 self._nodes[beam.node_b], beam.orientation)
        transform = _block_diagonal(rotation, 4)
        # local -> global: k_g = T^T k_l T, with T mapping global
        # displacements onto local ones
        k = transform.T @ _beam_stiffness(beam.material, beam.section, length) @ transform
        m = transform.T @ _beam_mass(beam.material, beam.section, length) @ transform
        rows = np.r_[index[beam.node_a]:index[beam.node_a] + 6,
                     index[beam.node_b]:index[beam.node_b] + 6]
        grid = np.ix_(rows, rows)
        stiffness[grid] += k
        mass[grid] += m

    for plate in self.plates:
        corners = [self._nodes[n] for n in plate.nodes]
        rotation, a, b = _plate_frame(*corners)
        transform = _block_diagonal(rotation, 8)
        k_local, m_local = _plate_matrices(plate.material,
                                           plate.thickness, a, b)
        k = transform.T @ k_local @ transform
        m = transform.T @ m_local @ transform
        rows = np.concatenate([np.arange(index[n], index[n] + 6)
                               for n in plate.nodes])
        grid = np.ix_(rows, rows)
        stiffness[grid] += k
        mass[grid] += m

    for item in self.masses:
        start = index[item.node]
        for offset in range(3):
            mass[start + offset, start + offset] += item.mass
        for offset, inertia in enumerate(item.inertia):
            mass[start + 3 + offset, start + 3 + offset] += inertia

    # 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

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 recognised 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 recognised 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)
    centre = np.mean(np.array([self._nodes[k] for k in solved]), axis=0)
    for i, node in enumerate(solved):
        offset = self._nodes[node] - centre
        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] = (),
) -> 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. The remaining problem is the symmetric generalized one, K phi = lambda M phi, solved by factoring M (Cholesky), reducing to a standard symmetric problem and transforming back — which is what makes the shapes come out mass-normalized to machine precision rather than normalized and then rescaled.

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.

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] = ()) -> 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. The remaining problem is the
    symmetric generalized one, K phi = lambda M phi, solved by
    factoring M (Cholesky), reducing to a standard symmetric problem
    and transforming back — which is what makes the shapes come out
    mass-normalized to machine precision rather than normalized and
    then rescaled.

    `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.
    """
    mass, stiffness = self.matrices()
    free = self._free_dofs(fixed)
    if not len(free):
        raise ValueError('every degree of freedom is fixed')
    reduced_m = mass[np.ix_(free, free)]
    reduced_k = stiffness[np.ix_(free, free)]

    try:
        factor = np.linalg.cholesky(reduced_m)
    except np.linalg.LinAlgError:
        starved = [self.dof_strings()[free[i]]
                   for i in np.flatnonzero(np.diag(reduced_m) <= 0.0)]
        raise ValueError(
            'the mass matrix is not positive definite, so these degrees '
            'of freedom carry no mass and no rotary inertia: '
            + (', '.join(starved[:6]) or 'none on the diagonal, so the '
               'model is nearly a mechanism')) from None

    # M = L L^T, so K phi = lambda M phi becomes A y = lambda y with
    # A = L^-1 K L^-T and phi = L^-T y. y orthonormal then gives
    # phi^T M phi = y^T y = I exactly, which is the normalization we
    # want and not something applied afterwards.
    temporary = np.linalg.solve(factor, reduced_k)
    standard = np.linalg.solve(factor, temporary.T).T
    eigenvalues, vectors = np.linalg.eigh((standard + standard.T) / 2.0)
    shapes = np.linalg.solve(factor.T, vectors)

    full = np.zeros((self.num_dof, shapes.shape[1]), dtype=np.float64)
    full[free] = shapes
    # 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 '')

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)
    # 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 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 block of the part they belong to, so the
    # geometry states its own regions and a saved file can rebuild the
    # structure. Without this the blocks 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 block.
    parts, blocks = {}, []
    for part in ([b.group for b in (self.beams if beams else ())]
                 + [p.group for p in self.plates]
                 + [f.group for f in self.faces]):
        blocks.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_block=blocks or None,
        block_id=list(parts.values()) or None,
        block_name=list(parts) or None,
        length_unit=self.length_unit)

Functions:

connected_pieces

connected_pieces(
    neighbours: 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(neighbours: 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 neighbours:
        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(neighbours[node] - seen)
        found.append(sorted(piece))
    return sorted(found, key=len, reverse=True)