visualdynamics.demo.drone

drone

The demonstration airframe: a quadcopter, built with visualdynamics.

This is the model every drone fixture comes from — the truth set, the reduced test geometry, and the plant the Rattlesnake runs fly against. It is built entirely with visualdynamics.fem — like every fixture model since the sdynpy demo articles left the repository: a demonstration of the toolset should be made by the toolset.

It builds a model and returns it; it writes no files. testdata/'s generate_drone.py, in the visualdynamics-generators repository, is what turns it into fixtures on disk.

from visualdynamics.demo import drone

model = drone.build()
shapes = model.eigensolution(maximum_frequency=400, damping=0.01)

Run it directly to print what it weighs and where its modes are:

python3 -m visualdynamics.demo.drone

The geometry is the model. Every surface is swept from a profile along a path, every edge of every face becomes a beam, and the mass is shared equally over the nodes. There is nothing derived and nothing tied: what you see is what is solved, node for node.

That is the second design. The first was a beam skeleton with surfaces hung off it on rigid ties — nodes that were drawn but not solved, their motion extrapolated from the beam they rode. It worked, and it was cheap, and it was a picture of a fidelity the model did not have: the arm walls looked meshed and were carrying nothing. Letting the drawing be the structure is simpler to write, simpler to explain, and honest about what it is.

Beams on the edges, not springs. A spring on each edge leaves a quad free to shear and a flat panel free to fold, since neither changes any edge length, and the model comes back a mechanism with a zero-frequency mode per panel. A beam carries moment, so the same wireframe stands up. Proved on a cube before anything else was built: six rigid-body modes, then 121.7 Hz.

The members are massless and the nodes carry it all. Equal shares of the total, so the mass follows the model rather than a section table — which also means it follows the mesh, and the mesh is kept roughly even for that reason. Each node also gets a rotary inertia, without which the rotational degrees of freedom carry nothing and the mass matrix will not factorize.

The numbers are plausible, not surveyed. One section size sets the whole frame's stiffness, chosen to put the elastic modes in a band a shaker test would use. This is a demonstration article, not a drone anyone has weighed.

Classes:

Name	Description
Shape	A drawing: nodes and faces, and nothing else.

Functions:

Name	Description
ring	Points round a circle, in a sweep's own cross-section plane.
sweep	Sweep cross-sections along a path, as rings of nodes joined by quads.
bar	A flat rectangular member between two points, in a stated plane.
blade_profile	A blade section: a thin cambered shape, rotated to its pitch.
bridge	Quads between two rings that already exist — a shared joint.
girder	A hollow box swept along a path, with windows cut through it.
build	The quadcopter, at whatever mesh density is asked for.
draw	The airframe as a drawing — nodes and faces, no structure at all.
part_of	Which part of the airframe a node belongs to.
instrumented	The nodes a modal survey of this airframe would put sensors on.
describe	What it weighs and where its modes are, by what moves in each.

Classes

Shape

Shape(first: int = 1000, tolerance: float = 1e-06)

A drawing: nodes and faces, and nothing else.

No members, no mass, no elements of any kind — the whole airframe is built into one of these, and fem.Model.from_geometry is what turns it into a structure afterwards. Keeping the two apart is the point: the geometry is authored once and the physics is derived from it, so there is no second description to disagree with the first.

Anything landing on an existing node is welded to it. Parts are drawn independently — a strut between two ring nodes, a nacelle standing on an arm — and where two of them meet they put a node in the same place twice. Drawn, that looks joined; solved, it is not, and each loose piece brings six zero-frequency modes. The airframe came back with 210 against the six it should have before this was here: 24 truss struts, 4 nacelles, 4 legs and 2 camera mounts, every one a free body.

Methods:

Name	Description
drop_face	Remove a face by the nodes it spans, if it is there.
pieces	The drawing's disconnected parts, largest first.
orient	Wind every face the same way round, and turn the lot outward.
prune	Drop nodes no face uses, and say how many there were.
geometry	What was drawn, as a Geometry: quads and triangles, no lines.

Source code in src/visualdynamics/demo/drone.py

def __init__(self, first: int = 1000,
             tolerance: float = 1e-6) -> None:
    self.tolerance: float = tolerance
    self.xyz: dict[int, np.ndarray] = {}
    self.group: dict[int, str] = {}
    self.faces: list[tuple[list[int], int]] = []
    self.face_group: list[str] = []
    self._at: dict[tuple, int] = {}
    self._next = first

Methods:

drop_face

drop_face(nodes: Sequence[int]) -> bool

Remove a face by the nodes it spans, if it is there.

Source code in src/visualdynamics/demo/drone.py

def drop_face(self, nodes: Sequence[int]) -> bool:
    """Remove a face by the nodes it spans, if it is there."""
    wanted = set(nodes)
    for index, (existing, _colour) in enumerate(self.faces):
        if set(existing) == wanted:
            self.faces.pop(index)
            return True
    return False

pieces

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

The drawing's disconnected parts, largest first.

Asked of the drawing rather than of a solved model, because that is when it can still be fixed cheaply. Two surfaces that touch are not joined unless they share nodes or something spans them, and the cost of finding out later is six zero-frequency modes per loose piece with nothing on screen to say which.

Source code in src/visualdynamics/demo/drone.py

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

    Asked of the drawing rather than of a solved model, because that is
    when it can still be fixed cheaply. Two surfaces that touch are not
    joined unless they share nodes or something spans them, and the
    cost of finding out later is six zero-frequency modes per loose
    piece with nothing on screen to say which.
    """
    neighbours: dict[int, set[int]] = {n: set() for n in self.xyz}
    for nodes, _color in self.faces:
        for k, node in enumerate(nodes):
            other = nodes[(k + 1) % len(nodes)]
            neighbours[node].add(other)
            neighbours[other].add(node)
    return fem.connected_pieces(neighbours)

orient

orient() -> int

Wind every face the same way round, and turn the lot outward.

Windings were being corrected girder by girder — reverse if the frame was flipped, reverse again if the rows ran the other way — and it never came right, because a part's winding depends on how it was built and the answer wanted is a property of the finished surface. So this asks the surface instead: walk face to face over shared edges, and where two neighbours traverse their shared edge the same way round, one of them is inside out. Then check the signed volume and turn everything over if the whole shell ended up pointing in.

Faces whose normals point into the solid show as dark patches, and would colour by displacement from the wrong side.

Source code in src/visualdynamics/demo/drone.py

def orient(self) -> int:
    """Wind every face the same way round, and turn the lot outward.

    Windings were being corrected girder by girder — reverse if the
    frame was flipped, reverse again if the rows ran the other way —
    and it never came right, because a part's winding depends on how
    it was *built* and the answer wanted is a property of the finished
    surface. So this asks the surface instead: walk face to face over
    shared edges, and where two neighbours traverse their shared edge
    the same way round, one of them is inside out. Then check the
    signed volume and turn everything over if the whole shell ended up
    pointing in.

    Faces whose normals point into the solid show as dark patches, and
    would colour by displacement from the wrong side.
    """
    edges: dict[frozenset, list[int]] = {}
    for index, (nodes, _colour) in enumerate(self.faces):
        for k, node in enumerate(nodes):
            edges.setdefault(
                frozenset((node, nodes[(k + 1) % len(nodes)])),
                []).append(index)

    def runs(index: int, a: int, b: int) -> bool:
        """Does face `index` traverse this edge from a to b?"""
        nodes = self.faces[index][0]
        for k, node in enumerate(nodes):
            if node == a and nodes[(k + 1) % len(nodes)] == b:
                return True
        return False

    seen, flipped = set(), 0
    for start in range(len(self.faces)):
        if start in seen:
            continue
        seen.add(start)
        stack = [start]
        while stack:
            index = stack.pop()
            nodes = self.faces[index][0]
            for k, node in enumerate(nodes):
                other = nodes[(k + 1) % len(nodes)]
                for neighbour in edges[frozenset((node, other))]:
                    if neighbour in seen:
                        continue
                    seen.add(neighbour)
                    # neighbours agree when they cross the shared edge
                    # in opposite directions; agreeing means one is
                    # wound the wrong way round
                    if runs(neighbour, node, other):
                        face, colour = self.faces[neighbour]
                        self.faces[neighbour] = (face[::-1], colour)
                        flipped += 1
                    stack.append(neighbour)

    volume = 0.0
    for nodes, _colour in self.faces:
        points = [self.xyz[n] for n in nodes]
        for k in range(1, len(points) - 1):
            volume += float(np.dot(points[0],
                                   np.cross(points[k], points[k + 1])))
    if volume < 0.0:
        self.faces = [(nodes[::-1], colour) for nodes, colour in self.faces]
        flipped += len(self.faces)
    return flipped

prune

prune() -> int

Drop nodes no face uses, and say how many there were.

A window more than one cell tall leaves the nodes strictly inside it belonging to nothing: the edges between two hole cells get no wall, and there is no face on either side. Rather than forbid tall windows — which is what keeps a truss looking like a slotted plate — the drawing simply forgets the nodes nothing referred to.

Source code in src/visualdynamics/demo/drone.py

def prune(self) -> int:
    """Drop nodes no face uses, and say how many there were.

    A window more than one cell tall leaves the nodes strictly inside
    it belonging to nothing: the edges between two hole cells get no
    wall, and there is no face on either side. Rather than forbid tall
    windows — which is what keeps a truss looking like a slotted plate
    — the drawing simply forgets the nodes nothing referred to.
    """
    used = {n for nodes, _color in self.faces for n in nodes}
    loose = [n for n in self.xyz if n not in used]
    for node in loose:
        key = tuple(round(float(v) / self.tolerance)
                    for v in self.xyz[node])
        self._at.pop(key, None)
        self.xyz.pop(node)
        self.group.pop(node, None)
    return len(loose)

geometry

geometry(length_unit: str = 'm') -> Geometry

What was drawn, as a Geometry: quads and triangles, no lines.

Each face goes into the block of the part that drew it, so the geometry says by itself which region is which — and a saved file can rebuild the structure without anything passed alongside it.

Blocks are named for the whole part, sides and all ('arm front left', not 'arm'). Collapsing them to the part alone reads more tidily and loses the one thing the sensor set is picked by: which of the four arms a node is on.

Source code in src/visualdynamics/demo/drone.py

def geometry(self, length_unit: str = 'm') -> Geometry:
    """What was drawn, as a Geometry: quads and triangles, no lines.

    Each face goes into the block of the part that drew it, so the
    geometry says by itself which region is which — and a saved file
    can rebuild the structure without anything passed alongside it.

    Blocks are named for the whole part, sides and all ('arm front
    left', not 'arm'). Collapsing them to the part alone reads more
    tidily and loses the one thing the sensor set is picked by: which
    of the four arms a node is on.
    """
    ids = sorted(self.xyz)
    parts, blocks = {}, []
    for group in self.face_group:
        part = group or 'body'
        blocks.append(parts.setdefault(part, len(parts) + 1))
    return Geometry(
        node_id=ids,
        node_xyz=np.array([self.xyz[n] for n in ids]),
        elem_conn=[nodes for nodes, _ in self.faces],
        elem_type=[44 if len(nodes) == 4 else 41
                   for nodes, _ in self.faces],
        elem_color=[color for _, color in self.faces],
        elem_block=blocks,
        block_id=list(parts.values()),
        block_name=list(parts),
        length_unit=length_unit)

Functions:

ring

ring(
    radius: float, sides: int, phase: float = 0.0
) -> list[tuple[float, float]]

Points round a circle, in a sweep's own cross-section plane.

Source code in src/visualdynamics/demo/drone.py

def ring(radius: float, sides: int,
         phase: float = 0.0) -> list[tuple[float, float]]:
    """Points round a circle, in a sweep's own cross-section plane."""
    return [(radius * math.cos(phase + 2.0 * math.pi * s / sides),
             radius * math.sin(phase + 2.0 * math.pi * s / sides))
            for s in range(sides)]

sweep

sweep(
    shape: Shape,
    path: Sequence[ArrayLike],
    profiles: Profile | Sequence[Profile],
    group: str,
    color: int,
    cap_start: bool = False,
    cap_end: bool = False,
) -> list[list[int]]

Sweep cross-sections along a path, as rings of nodes joined by quads.

profiles is one profile per station, so a tube can taper, swell or change shape along its length; pass a single profile to keep it uniform. Returns the rings.

Source code in src/visualdynamics/demo/drone.py

def sweep(shape: Shape, path: Sequence[ArrayLike],
          profiles: Profile | Sequence[Profile], group: str, color: int,
          cap_start: bool = False,
          cap_end: bool = False) -> list[list[int]]:
    """Sweep cross-sections along a path, as rings of nodes joined by quads.

    `profiles` is one profile per station, so a tube can taper, swell or
    change shape along its length; pass a single profile to keep it
    uniform. Returns the rings.
    """
    if not isinstance(profiles[0], (list, tuple)) or isinstance(
            profiles[0][0], (int, float)):
        stations: Sequence[Profile] = [profiles] * len(path)  # type: ignore[list-item]
    else:
        stations = profiles                                 # type: ignore[assignment]
    rings = []
    for (centre, side, up), profile in zip(_frames(path), stations):
        row = []
        for a, b in profile:
            point = centre + a * side + b * up
            row.append(shape.node(point, group))
        rings.append(row)
    for lower, upper in pairwise(rings):
        for s in range(len(lower)):
            t = (s + 1) % len(lower)
            shape.face([lower[s], lower[t], upper[t], upper[s]], color, group)
    for wanted, row, centre in ((cap_start, rings[0], path[0]),
                                (cap_end, rings[-1], path[-1])):
        if not wanted:
            continue
        hub = shape.node(np.asarray(centre, dtype=float), group)
        for s in range(len(row)):
            shape.face([hub, row[s], row[(s + 1) % len(row)]], color, group)
    return rings

bar

bar(
    shape: Shape,
    start: ArrayLike,
    end: ArrayLike,
    thick: float,
    depth: float,
    normal: ArrayLike,
    group: str,
    color: int,
) -> list[list[int]] | None

A flat rectangular member between two points, in a stated plane.

The frame is given rather than worked out, because a truss has members at every angle and letting each pick its own reference direction lets the sections twist relative to one another. normal is the plane the truss lies in, so every bar in one arm is the same slab of material seen from a different angle — which is what a plate truss is.

Source code in src/visualdynamics/demo/drone.py

def bar(shape: Shape, start: ArrayLike, end: ArrayLike, thick: float, depth: float,
        normal: ArrayLike, group: str,
        color: int) -> list[list[int]] | None:
    """A flat rectangular member between two points, in a stated plane.

    The frame is given rather than worked out, because a truss has members
    at every angle and letting each pick its own reference direction lets
    the sections twist relative to one another. `normal` is the plane the
    truss lies in, so every bar in one arm is the same slab of material
    seen from a different angle — which is what a plate truss is.
    """
    start = np.asarray(start, dtype=float)
    end = np.asarray(end, dtype=float)
    tangent = end - start
    length = float(np.linalg.norm(tangent))
    if length < 1e-9:
        return None
    tangent = tangent / length
    across = np.asarray(normal, dtype=float)
    across = across / np.linalg.norm(across)
    within = np.cross(tangent, across)
    within = within / np.linalg.norm(within)
    rings = []
    for point in (start, end):
        rings.append([shape.node(point + a * across + b * within, group)
                      for a, b in ((-thick / 2, -depth / 2),
                                   (thick / 2, -depth / 2),
                                   (thick / 2, depth / 2),
                                   (-thick / 2, depth / 2))])
    for k in range(4):
        t = (k + 1) % 4
        shape.face([rings[0][k], rings[0][t], rings[1][t], rings[1][k]], color)
    # the ends are closed onto the joint itself, which is a node the
    # neighbouring bars also close onto — so the joint is shared, not
    # stitched, and the truss is one surface rather than a pile of sticks
    for row, point in ((rings[0], start), (rings[1], end)):
        hub = shape.node(point, group)
        for k in range(4):
            shape.face([hub, row[k], row[(k + 1) % 4]], color)
    return rings

blade_profile

blade_profile(
    chord: float, thick: float, twist: float
) -> list[tuple[float, float]]

A blade section: a thin cambered shape, rotated to its pitch.

Twist is applied to the profile rather than to the sweep, because the path is a straight radial line and its frame does not turn — a blade that is flat at the tip and coarse at the root is the whole reason a propeller looks like one rather than like a paddle.

Source code in src/visualdynamics/demo/drone.py

def blade_profile(chord: float, thick: float,
                  twist: float) -> list[tuple[float, float]]:
    """A blade section: a thin cambered shape, rotated to its pitch.

    Twist is applied to the profile rather than to the sweep, because the
    path is a straight radial line and its frame does not turn — a blade
    that is flat at the tip and coarse at the root is the whole reason a
    propeller looks like one rather than like a paddle.
    """
    points = ((-0.50, 0.00), (-0.30, -0.42), (0.10, -0.50), (0.42, -0.26),
              (0.50, 0.00), (0.30, 0.40), (-0.10, 0.46), (-0.38, 0.28))
    angle = math.radians(twist)
    cos, sin = math.cos(angle), math.sin(angle)
    out = []
    for a, b in points:
        u, v = a * chord, b * thick
        out.append((u * cos - v * sin, u * sin + v * cos))
    return out

bridge

bridge(
    shape: Shape,
    lower: Sequence[int],
    upper: Sequence[int],
    color: int,
) -> None

Quads between two rings that already exist — a shared joint.

Source code in src/visualdynamics/demo/drone.py

def bridge(shape: Shape, lower: Sequence[int], upper: Sequence[int],
           color: int) -> None:
    """Quads between two rings that already exist — a shared joint."""
    for s in range(len(lower)):
        t = (s + 1) % len(lower)
        shape.face([lower[s], lower[t], upper[t], upper[s]], color)

girder

girder(
    shape: Shape,
    left: Sequence[int],
    right: Sequence[int],
    path: Sequence[ArrayLike],
    half_height: float,
    half_thick: float,
    windows: Callable[[int, int], bool] | None,
    group: str,
    color: int,
    cap_end: bool = True,
    shear: float = 0.0,
    end_left: Sequence[int] | None = None,
    end_right: Sequence[int] | None = None,
) -> tuple[list[list[int]], list[list[int]]]

A hollow box swept along a path, with windows cut through it.

One closed surface: two flat sides, a strip along the top and bottom, a cap at the far end, and a wall round every window joining one side to the other. Nothing overlaps anything and there is no face inside the solid — which is the whole point. Built as separate bars welded at their centres, as this was first, a truss is a heap of interpenetrating boxes with their end caps buried in the joints, and colouring it by displacement shows the insides through the skin.

left and right are the rows of nodes the root starts from, so a girder can grow out of a surface that already exists rather than being parked against it; end_left and end_right do the same at the far end, so a member can grow into something as well as out of it.

shear leans the interior node rows alternately along the run, which turns the openings from upright rectangles into a zigzag — the difference between a slotted plate and a truss.

Source code in src/visualdynamics/demo/drone.py

def girder(shape: Shape, left: Sequence[int], right: Sequence[int],
           path: Sequence[ArrayLike], half_height: float, half_thick: float,
           windows: Callable[[int, int], bool] | None, group: str,
           color: int, cap_end: bool = True, shear: float = 0.0,
           end_left: Sequence[int] | None = None,
           end_right: Sequence[int] | None = None
           ) -> tuple[list[list[int]], list[list[int]]]:
    """A hollow box swept along a path, with windows cut through it.

    One closed surface: two flat sides, a strip along the top and bottom,
    a cap at the far end, and a wall round every window joining one side to
    the other. Nothing overlaps anything and there is no face inside the
    solid — which is the whole point. Built as separate bars welded at
    their centres, as this was first, a truss is a heap of interpenetrating
    boxes with their end caps buried in the joints, and colouring it by
    displacement shows the insides through the skin.

    `left` and `right` are the rows of nodes the root starts from, so a
    girder can grow out of a surface that already exists rather than being
    parked against it; `end_left` and `end_right` do the same at the far
    end, so a member can grow *into* something as well as out of it.

    `shear` leans the interior node rows alternately along the run, which
    turns the openings from upright rectangles into a zigzag — the
    difference between a slotted plate and a truss.
    """
    rows = len(left) - 1
    lo = [list(left)]
    hi = [list(right)]
    # One reference for the whole run, chosen from where it ends up rather
    # than from each step. Picking it per station lets it swap axes the
    # moment a path tips past vertical, and the section spins a quarter
    # turn between one station and the next — which is what put a visible
    # twist in the legs, where they steepen towards the foot.
    overall = np.asarray(path[-1], dtype=float) - np.asarray(path[0],
                                                             dtype=float)
    overall = overall / np.linalg.norm(overall)
    # The cross-section direction comes from the panel this girder starts
    # on, not from a global axis. Taken from a fixed reference it is
    # whatever that axis gives, and any panel not lined up with it starts
    # rotated -- an arm at sixty degrees, a strap under a curved belly.
    # A sign flip cannot correct a rotation, which is why correcting the
    # handedness left the nacelle roots and the payload straps twisted.
    seam0 = (np.asarray(shape.position(right[0]), dtype=float)
             - np.asarray(shape.position(left[0]), dtype=float))
    seam0 = seam0 - float(np.dot(seam0, overall)) * overall
    if float(np.linalg.norm(seam0)) < 1e-9:
        seam0 = np.cross(overall, np.array([0.0, 0.0, 1.0]))
    seam0 = seam0 / np.linalg.norm(seam0)
    # Which way round the root panel already runs. The frame's own `across`
    # is whatever the cross product gives, and when it points from right to
    # left the first bay is built inside out: its quads cross over between
    # the panel and the next station, which is the twist that shows where
    # an arm meets the body. Ask the panel rather than assume.
    # And which way the rows run. Flipping `across` flips `up` with it, so
    # correcting left-for-right silently turns the rows upside down: the
    # girder counts j upward while the panel it started from counts it
    # downward, and every one of them is built inverted. That is what put
    # the legs on top of the arms -- row 0 of an arm is meant to be its
    # underside, because row 0 of the waist panel is.
    rise = (np.asarray(shape.position(left[-1]), dtype=float)
            - np.asarray(shape.position(left[0]), dtype=float))
    climbs = (1.0 if float(np.dot(np.cross(seam0, overall), rise)) >= 0.0
              else -1.0)
    for k in range(1, len(path)):
        centre = np.asarray(path[k], dtype=float)
        tangent = np.asarray(path[k], dtype=float) - np.asarray(path[k - 1],
                                                                dtype=float)
        tangent = tangent / np.linalg.norm(tangent)
        across = seam0 - float(np.dot(seam0, tangent)) * tangent
        across = across / np.linalg.norm(across)
        up = np.cross(across, tangent)
        h, t = half_height[k], half_thick[k]
        if k == len(path) - 1 and end_left is not None:
            # Sort the panel it lands on into this girder's own frame
            # rather than trusting the order it was handed in. The two
            # panels a graft joins were built by different parts, for
            # different reasons, and their rows have no reason to run the
            # same way: taken as given, the last bay connects row j of one
            # to row j of the other and crosses over. That is the twist in
            # the camera's mounts.
            ends = list(end_left) + list(end_right)
            middle = np.mean([shape.position(n) for n in ends], axis=0)

            def placed(node: int, middle: np.ndarray = middle,
                       across: np.ndarray = across,
                       up: np.ndarray = up) -> tuple[float, float]:
                offset = np.asarray(shape.position(node),
                                    dtype=float) - middle
                return (float(np.dot(offset, across)),
                        float(np.dot(offset, up * climbs)))

            near = sorted(ends, key=lambda n: placed(n)[0])
            half = len(ends) // 2
            lo.append(sorted(near[:half], key=lambda n: placed(n)[1]))
            hi.append(sorted(near[half:], key=lambda n: placed(n)[1]))
            continue
        lean = shear * (1.0 if k % 2 else -1.0)
        lo.append([shape.node(centre - across * t
                              + up * climbs * (-h + 2 * h * j / rows)
                              + tangent * lean * _bulge(j, rows), group)
                   for j in range(rows + 1)])
        hi.append([shape.node(centre + across * t
                              + up * climbs * (-h + 2 * h * j / rows)
                              + tangent * lean * _bulge(j, rows), group)
                   for j in range(rows + 1)])

    # Flipping the frame reverses which way round a quad runs, so the
    # faces that follow it have to be reversed too or their normals point
    # into the solid — which renders as a dark band at the join and would
    # colour by displacement from the inside.
    def face(nodes: Sequence[int], colour: int | None = None) -> None:
        # Wound however it comes out; Shape.orient() settles the whole
        # surface afterwards, which is the only level at which the answer
        # is well defined.
        shape.face(nodes, color if colour is None else colour, group)

    holes = {(k, j) for k in range(len(path) - 1) for j in range(rows)
             if windows(k, j)}
    for k in range(len(path) - 1):
        for j in range(rows):
            if (k, j) in holes:
                continue
            face([lo[k][j], lo[k][j + 1], lo[k + 1][j + 1], lo[k + 1][j]])
            face([hi[k][j], hi[k + 1][j], hi[k + 1][j + 1], hi[k][j + 1]])
        # the narrow strips along the top and the bottom of the box
        face([lo[k][0], lo[k + 1][0], hi[k + 1][0], hi[k][0]])
        face([lo[k][rows], hi[k][rows], hi[k + 1][rows], lo[k + 1][rows]])
    # every window is walled from one side to the other, so the hole is a
    # hole through a solid rather than two gaps in two skins
    for k, j in holes:
        for a, b, neighbour in (((k, j), (k, j + 1), (k - 1, j)),
                                ((k + 1, j + 1), (k + 1, j), (k + 1, j)),
                                ((k + 1, j), (k, j), (k, j - 1)),
                                ((k, j + 1), (k + 1, j + 1), (k, j + 1))):
            if neighbour in holes:
                continue
            face([lo[a[0]][a[1]], lo[b[0]][b[1]],
                  hi[b[0]][b[1]], hi[a[0]][a[1]]])
    if cap_end:
        for j in range(rows):
            face([lo[-1][j], hi[-1][j], hi[-1][j + 1], lo[-1][j + 1]])
    return lo, hi

build

build(
    sides: int = 12,
    arm_stations: int = 9,
    leg_stations: int = 5,
    body_rings: int = 3,
    total_mass: float = TOTAL_MASS,
) -> Model

The quadcopter, at whatever mesh density is asked for.

sides is how many facets go round every swept tube and round the body, and is what most of the node count comes from. The defaults are about a thousand nodes, which the dense eigensolver clears in half a minute; tests build it coarser, since what they check is the airframe and not the mesh.

Source code in src/visualdynamics/demo/drone.py

def build(sides: int = 12, arm_stations: int = 9, leg_stations: int = 5,
          body_rings: int = 3, total_mass: float = TOTAL_MASS) -> fem.Model:
    """The quadcopter, at whatever mesh density is asked for.

    `sides` is how many facets go round every swept tube and round the body,
    and is what most of the node count comes from. The defaults are about a
    thousand nodes, which the dense eigensolver clears in half a minute;
    tests build it coarser, since what they check is the airframe and not
    the mesh.
    """
    shape = draw(sides, arm_stations, leg_stations, body_rings)
    shape.prune()
    shape.orient()
    # The drawing becomes the structure here and only here: a member along
    # every edge of every face, and the mass shared over the nodes. There
    # is one description of this aircraft, and the physics is derived from
    # it rather than written beside it.
    # No groups passed: the geometry's own blocks say which part each face
    # belongs to, so a model built here and a model rebuilt from a saved
    # file are the same model. They were not while the parts arrived
    # alongside the drawing — 16 blade members read as frame.
    return fem.Model.from_geometry(
        shape.geometry(), FRAME, MEMBER, total_mass=total_mass,
        name='quadcopter', sections={'prop': BLADE})

draw

draw(
    sides: int = 12,
    arm_stations: int = 9,
    leg_stations: int = 5,
    body_rings: int = 3,
) -> Shape

The airframe as a drawing — nodes and faces, no structure at all.

Source code in src/visualdynamics/demo/drone.py

def draw(sides: int = 12, arm_stations: int = 9, leg_stations: int = 5,
         body_rings: int = 3) -> Shape:
    """The airframe as a drawing — nodes and faces, no structure at all."""
    shape = Shape()
    rows = 5
    waist, _apex, floor, belly = _body(shape, sides, body_rings, rows)
    # each arm takes one panel of the waist band as its root, so the panel
    # itself is not skinned over -- the arm is what closes it
    taken = {round(sides * angle / 360.0) % sides for _n, angle, _r in ARMS}
    _skin_waist(shape, waist, sides, rows, taken)
    for (name, angle, radius), index in zip(ARMS, sorted(taken)):
        lo, hi, along = _arm(shape, waist, sides, rows, index, radius,
                             arm_stations, name)
        _leg(shape, lo, hi, along, leg_stations,
             LEG_DROP * (0.85 if 'front' in name else 1.0), name)
    _payload(shape, floor, belly, sides)
    return shape

part_of

part_of(model: Model, node: int) -> str

Which part of the airframe a node belongs to.

Source code in src/visualdynamics/demo/drone.py

def part_of(model: fem.Model, node: int) -> str:
    """Which part of the airframe a node belongs to."""
    group = model.group(node)
    head = group.split()[0] if group else ''
    return PARTS.get(head, 'body')

instrumented

instrumented(model: Model) -> dict[str, list[int]]

The nodes a modal survey of this airframe would put sensors on.

Found from the model rather than written down, because writing them down has been wrong twice: the numbering moves whenever the mesh does, and a generator holding stale ids produces a test of nowhere.

Source code in src/visualdynamics/demo/drone.py

def instrumented(model: fem.Model) -> dict[str, list[int]]:
    """The nodes a modal survey of this airframe would put sensors on.

    Found from the model rather than written down, because writing them
    down has been wrong twice: the numbering moves whenever the mesh does,
    and a generator holding stale ids produces a test of nowhere.
    """
    def nearest(candidates: Sequence[int], point: ArrayLike) -> int:
        point = np.asarray(point, dtype=float)
        return min(candidates, key=lambda n: float(np.linalg.norm(
            np.asarray(model.position(n))[:len(point)] - point)))

    groups: dict[str, list[int]] = {}
    for node in model.node_ids:
        groups.setdefault(model.group(node), []).append(node)

    found: dict[str, list[int]] = {'motors': [], 'arms': [], 'feet': []}
    for name, angle, radius in ARMS:
        theta = math.radians(angle)
        out = np.array([math.cos(theta), math.sin(theta)])
        nacelle = groups.get(f'nacelle {name}', [])
        if nacelle:
            found['motors'].append(max(nacelle,
                                       key=lambda n: model.position(n)[2]))
        arm = groups.get(f'arm {name}', [])
        if arm:
            found['arms'].append(nearest(arm, out * radius * 0.6))
        leg = groups.get(f'leg {name}', [])
        if leg:
            found['feet'].append(min(leg, key=lambda n: model.position(n)[2]))

    # The shell, not the canopy block: the canopy is the apex cap alone —
    # eight triangles — and the ring at 0.9 of the body radius is the
    # waist below it. Asking one block for both put all four body sensors
    # and the centre one on the same node.
    body = [n for n in model.node_ids if part_of(model, n) == 'body']
    found['body'] = [nearest(body, (BODY_R * 0.9 * math.cos(math.radians(a)),
                                    BODY_R * 0.9 * math.sin(math.radians(a))))
                     for a in (0, 90, 180, 270)]
    found['centre'] = [max(body, key=lambda n: model.position(n)[2])]
    found['payload'] = [
        min(groups.get('battery', [0]), key=lambda n: model.position(n)[2]),
        min(groups.get('camera', [0]), key=lambda n: model.position(n)[0]),
    ]
    return found

describe

describe(
    model: Model,
    maximum_frequency: float = 500.0,
    damping: float = 0.01,
) -> None

What it weighs and where its modes are, by what moves in each.

Source code in src/visualdynamics/demo/drone.py

def describe(model: fem.Model, maximum_frequency: float = 500.0,
             damping: float = 0.01) -> None:
    """What it weighs and where its modes are, by what moves in each."""
    shapes = model.eigensolution(maximum_frequency=maximum_frequency,
                                 damping=damping)
    mass, _ = model.matrices()
    ids = model.node_ids
    order = ['body', 'arm', 'nacelle', 'leg', 'battery', 'camera']

    print(f'{model.num_nodes} nodes, {model.num_dof} DOF, '
          f'{len(model.beams)} members, {len(model.faces)} faces')
    print(f'{model.total_mass * 1000:.0f} g over '
          f'{len(ids)} nodes')
    print()
    print(f'{"f (Hz)":>9}  ' + '  '.join(f'{f:>8s}' for f in order))
    for frequency, shape in zip(shapes.frequency, shapes.shape_matrix):
        if frequency == 0.0:
            continue
        weighted = mass @ shape[:model.num_dof]
        share: dict[str, float] = {}
        for i, node in enumerate(ids):
            key = part_of(model, node)
            share[key] = share.get(key, 0.0) + float(
                np.dot(shape[6 * i:6 * i + 6], weighted[6 * i:6 * i + 6]))
        total = sum(share.values())
        print(f'{frequency:9.2f}  '
              + '  '.join(f'{100 * share.get(k, 0.0) / total:7.1f}%'
                          for k in order)
              + f'   {max(share, key=share.get)}')