visualdynamics.core.tetmesh

tetmesh

Tetrahedral meshing of a geometry's solids, from their recipe.

Mesh Solids meshes an extrusion into bricks (solids.extrusion_mesh). A part that is not one, a hole drilled at an angle, cylinders along three axes, is meshed here into tetrahedra, from the recipe's signed distance (solids.distance) alone: nothing about the shape is assumed beyond where its surface is and which way is out.

The method is the one of Persson and Strang (2004). Points are laid on a body-centered lattice inside the solid, with points fixed along its sharp edges and at its corners; every edge of their Delaunay triangulation pushes its ends apart when shorter than the length aimed at; a point pushed out of the solid is put back onto its surface along the distance's gradient; and the triangulation is made again when the points have moved enough to need it. The relaxation ends when the points have settled, and the tetrahedra whose centroids fall outside the solid are dropped.

Sharp edges are found from the recipe: each block gives six planes, each cylinder a wall and two caps, and an edge is a curve where two of them cross on the solid's surface, a corner a point where three do. Points are placed along the curves a mesh size apart and held there, so the mesh keeps every edge the recipe has and rounds nothing off.

The elements are ten-node tetrahedra (type 118): the corners, then the midside nodes of edges 1-2, 2-3, 1-3, 1-4, 2-4, 3-4, those on a curved wall moved onto it. Four-node ones (type 111) are made on request.

Functions:

Name	Description
tet_count	About how many tetrahedra tet_mesh would make for these parts
tet_mesh	The parts as tetrahedra, ready to join to a geometry
sharp_edges	The points tet_mesh holds fixed: the parts' corners, and points

Classes

Functions:

tet_count

tet_count(geometry: Any, ids, size: float) -> int

About how many tetrahedra tet_mesh would make for these parts at size: their volume, read by sampling, over the cube of the size, times what a relaxed lattice packs into that. Within a quarter or so; the pane reads it to decide whether to preview.

Parameters:

Name	Type	Description	Default
geometry	Geometry	Whose solids.	required
ids	sequence of int	The parts, by id.	required
size	float	The element size, as the geometry holds its coordinates.	required

Returns:

Type	Description
int

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

def tet_count(geometry: Any, ids, size: float) -> int:
    """About how many tetrahedra `tet_mesh` would make for these parts
    at `size`: their volume, read by sampling, over the cube of the
    size, times what a relaxed lattice packs into that. Within a
    quarter or so; the pane reads it to decide whether to preview.

    Parameters
    ----------
    geometry : Geometry
        Whose solids.
    ids : sequence of int
        The parts, by id.
    size : float
        The element size, as the geometry holds its coordinates.

    Returns
    -------
    int
    """
    ids = [int(i) for i in ids]
    low, high = _extent(geometry, ids)
    rail = [np.linspace(a, b, 24) for a, b in zip(low, high)]
    at = np.stack(np.meshgrid(*rail, indexing='ij'), axis=-1).reshape(-1, 3)
    inside = (_distance(geometry, ids, at) < 0.0).mean()
    volume = float(inside * np.prod(high - low))
    return round(TETS_PER_CELL * volume / float(size) ** 3)

tet_mesh

tet_mesh(geometry: Any, ids, size: float, quadratic: bool = True, progress: Any = None, keep_slivers: bool = False) -> tuple[Any, dict]

The parts as tetrahedra, ready to join to a geometry (mesh.join).

Each part gets one element group, named for it (or for its kind and id when it has no name); an element goes to the part whose recipe holds its centroid, the later part where parts overlap. Parts that touch share their nodes there.

Refused: a solid that is part of another, a size that is not positive, and a mesh of more than MESH_LIMIT tetrahedra by tet_count's estimate.

Parameters:

Name	Type	Description	Default
geometry	Geometry	Whose solids.	required
ids	sequence of int	The parts, by id; each a solid nothing else uses.	required
size	float	The element size aimed at, as the geometry holds its coordinates: the edge length the relaxation settles to.	required
quadratic	bool	Ten-node tetrahedra; False for four-node ones.	True
progress	callable	Told (done, total) as the mesh goes: the edges, the points, each step of the relaxation, the elements.	None
keep_slivers	bool	Keep the slivers that lie in the surface rather than remove them. Removed, each leaves a shallow dimple, and the surface pinches along an edge there; kept, the mesh's outer faces are the surface whole, which is what a shell mesh of it wants (shells.shell_mesh).	False

Returns:

Type	Description
(Geometry, dict)	The tetrahedra, in the geometry's length unit; and what was made: 'tets' and 'nodes', the counts; 'steps', how many the relaxation took; 'quality', the least radius ratio (1 for a regular tetrahedron, 0 for a flat one); 'poor', how many are under 0.3; 'slivers', how many were removed from the surface; and 'shortest', the shortest edge.

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

def tet_mesh(geometry: Any, ids, size: float, quadratic: bool = True,
             progress: Any = None, keep_slivers: bool = False) -> tuple[Any, dict]:
    """The parts as tetrahedra, ready to join to a geometry
    (`mesh.join`).

    Each part gets one element group, named for it (or for its kind and
    id when it has no name); an element goes to the part whose recipe
    holds its centroid, the later part where parts overlap. Parts that
    touch share their nodes there.

    Refused: a solid that is part of another, a size that is not
    positive, and a mesh of more than `MESH_LIMIT` tetrahedra by
    `tet_count`'s estimate.

    Parameters
    ----------
    geometry : Geometry
        Whose solids.
    ids : sequence of int
        The parts, by id; each a solid nothing else uses.
    size : float
        The element size aimed at, as the geometry holds its
        coordinates: the edge length the relaxation settles to.
    quadratic : bool, default True
        Ten-node tetrahedra; False for four-node ones.
    progress : callable, optional
        Told ``(done, total)`` as the mesh goes: the edges, the points,
        each step of the relaxation, the elements.
    keep_slivers : bool, default False
        Keep the slivers that lie in the surface rather than remove
        them. Removed, each leaves a shallow dimple, and the surface
        pinches along an edge there; kept, the mesh's outer faces are
        the surface whole, which is what a shell mesh of it wants
        (`shells.shell_mesh`).

    Returns
    -------
    (Geometry, dict)
        The tetrahedra, in the geometry's length unit; and what was
        made: 'tets' and 'nodes', the counts; 'steps', how many the
        relaxation took; 'quality', the least radius ratio (1 for a
        regular tetrahedron, 0 for a flat one); 'poor', how many are
        under 0.3; 'slivers', how many were removed from the surface;
        and 'shortest', the shortest edge.
    """
    ids = [int(i) for i in ids]
    if not ids:
        raise ValueError('pick a solid to mesh')
    used = used_by(geometry)
    inside_another = [i for i in ids if i in used]
    if inside_another:
        raise ValueError(f'solid {inside_another[0]} is part of solid '
                         f'{used[inside_another[0]]}; mesh that one')
    size = float(size)
    if not size > 0.0:
        raise ValueError('an element size above zero')
    expected = tet_count(geometry, ids, size)
    if expected > MESH_LIMIT:
        raise ValueError(f'about {expected:,} tetrahedra at that size: use a larger one')
    ticker = Ticker(progress)
    ticker.add(MOST_STEPS + 3)

    def fd(points):
        return _distance(geometry, ids, points)

    low, high = _extent(geometry, ids)
    extent = float(np.linalg.norm(high - low))
    fixed = sharp_edges(geometry, ids, size)
    ticker.tick()
    points, nfixed = _seed(fd, fixed, low, high, size)
    ticker.tick()
    points, steps = _relax(fd, points, nfixed, size, extent, ticker)
    points, corners, quality, slivers = _tetrahedra(fd, points, nfixed, size, extent,
                                                    keep_slivers)
    if not len(corners):
        raise ValueError('the parts hold nothing to mesh')
    centroids = points[corners].mean(axis=1)
    groups = _owners(geometry, ids, centroids)
    corners, quality, groups = corners[groups > 0], quality[groups > 0], groups[groups > 0]
    if quadratic:
        points, conn = _quadratic(fd, _surfaces(geometry, ids, size), points, corners, size)
    else:
        conn = corners
    keep = np.unique(conn)
    renumber = np.zeros(len(points), dtype=int)
    renumber[keep] = np.arange(1, len(keep) + 1)
    edges = points[corners[:, [1, 2, 2, 3, 3, 3]]] - points[corners[:, [0, 0, 1, 0, 1, 2]]]
    shortest = float(np.linalg.norm(edges, axis=2).min()) if len(corners) else 0.0
    names, present = [], []
    for slot, part in enumerate(ids, start=1):
        if (groups == slot).any():
            item = solid(geometry, part)
            names.append(item.name.strip() or f'{item.kind.title()} {item.id}')
            present.append(slot)
    code = TET_CODES[conn.shape[1]]
    part = Geometry(node_id=np.arange(1, len(keep) + 1), node_xyz=points[keep],
                    elem_conn=list(renumber[conn]), elem_type=[code] * len(conn),
                    elem_group=[int(g) for g in groups], group_id=present,
                    group_name=names, length_unit=geometry.length_unit)
    ticker.tick()
    return part, {'tets': len(conn), 'nodes': len(keep), 'steps': steps,
                  'quality': float(quality.min()) if len(quality) else 0.0,
                  'poor': int((quality < 0.3).sum()), 'slivers': slivers,
                  'shortest': shortest}

sharp_edges

sharp_edges(geometry: Any, ids, size: float) -> ndarray

The points tet_mesh holds fixed: the parts' corners, and points about size apart along every sharp edge of their surface.

The surfaces are read from the recipe, each block's six planes and each cylinder's wall and caps, unbounded; a sharp edge is where two of them cross on the solid's surface, a corner where three do. Points seeded on one surface are moved onto the other and back until they sit on both, kept where the solid's distance is zero there, and thinned to the size.

Parameters:

Name	Type	Description	Default
geometry	Geometry	Whose solids.	required
ids	sequence of int	The parts, by id.	required
size	float	The mesh size: the spacing along an edge.	required

Returns:

Type	Description
ndarray	(n, 3), the corners first.

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

def sharp_edges(geometry: Any, ids, size: float) -> np.ndarray:
    """The points `tet_mesh` holds fixed: the parts' corners, and points
    about `size` apart along every sharp edge of their surface.

    The surfaces are read from the recipe, each block's six planes and
    each cylinder's wall and caps, unbounded; a sharp edge is where two
    of them cross on the solid's surface, a corner where three do.
    Points seeded on one surface are moved onto the other and back
    until they sit on both, kept where the solid's distance is zero
    there, and thinned to the size.

    Parameters
    ----------
    geometry : Geometry
        Whose solids.
    ids : sequence of int
        The parts, by id.
    size : float
        The mesh size: the spacing along an edge.

    Returns
    -------
    numpy.ndarray
        (n, 3), the corners first.
    """
    ids = [int(i) for i in ids]
    low, high = _extent(geometry, ids)
    extent = float(np.linalg.norm(high - low))
    tol = 1e-9 * extent

    def fd(points):
        return _distance(geometry, ids, points)

    surfaces = _surfaces(geometry, ids, size)
    fine = size / 3.0

    def crossing(which, seeds, rounds=30):
        """`seeds` moved onto every surface in `which` in turn until
        they lie on all, those on the solid's surface and within every
        surface's own extent kept."""
        q = seeds.copy()
        for _ in range(rounds):
            for k in which:
                q = _onto(surfaces[k][0], q, tol)
        on_all = np.all([np.abs(surfaces[k][0](q)) < tol for k in which], axis=0)
        on_all &= np.all([surfaces[k][4](q) for k in which], axis=0)
        on_all &= np.abs(fd(q)) < tol
        return q[on_all]

    curves, corners = [], []
    for i, (gi, seed_i, lo_i, hi_i, _within_i) in enumerate(surfaces):
        for j in range(i + 1, len(surfaces)):
            gj, _seed_j, lo_j, hi_j, _within_j = surfaces[j]
            lo, hi = np.maximum(lo_i, lo_j), np.minimum(hi_i, hi_j)
            if (hi < lo).any():
                continue
            seeds = seed_i(lo, hi, fine)
            if not len(seeds):
                continue
            seeds = seeds[np.abs(gj(seeds)) < 2.0 * size]
            if not len(seeds):
                continue
            found = crossing([i, j], seeds)
            if not len(found):
                continue
            curves.append(found)
            for k, (gk, _seed_k, lo_k, hi_k, _within_k) in enumerate(surfaces):
                if k in (i, j) or (hi_k < lo).any() or (lo_k > hi).any():
                    continue
                near = found[np.abs(gk(found)) < size]
                if len(near):
                    corners.append(crossing([i, j, k], near))
    kept: list[np.ndarray] = []
    if corners:
        for q in np.concatenate(corners):
            if not kept or cKDTree(np.array(kept)).query(q)[0] >= 0.25 * size:
                kept.append(q)
    for curve in curves:
        # along the curve, a point every 0.9 sizes of the way, and clear
        # of the corners and the curves before it by as much: one tree
        # over those, built once a curve
        others = cKDTree(np.array(kept)) if kept else None
        last = None
        for q in _along(curve):
            if last is not None and np.linalg.norm(q - last) < 0.9 * size:
                continue
            if others is not None and others.query(q)[0] < 0.9 * size:
                continue
            kept.append(q)
            last = q
    return np.array(kept).reshape(-1, 3)