From 4296116f41f5dcf0ef2939a47182cba46b132f4e Mon Sep 17 00:00:00 2001
From: Richard DeVaul <rdevaul@relativityspace.com>
Date: Mon, 6 Oct 2025 16:07:54 -0700
Subject: [PATCH] Addition of toplology analysis functions issolidclosed and
 volumeof

---
 examples/topology_demo.py | 180 +++++++++++++++++++++++
 src/yapcad/geom3d.py      | 191 +++++++++++++++++++++++-
 tests/test_geom3d.py      | 295 +++++++++++++++++++++++++++++++++++++-
 3 files changed, 664 insertions(+), 2 deletions(-)
 create mode 100755 examples/topology_demo.py

diff --git a/examples/topology_demo.py b/examples/topology_demo.py
new file mode 100755
index 0000000..132db7e
--- /dev/null
+++ b/examples/topology_demo.py
@@ -0,0 +1,180 @@
+#!/usr/bin/env python3
+"""
+Demonstration of yapCAD solid topology analysis functions.
+
+This script demonstrates the new issolidclosed() and volumeof() functions
+for analyzing 3D solid geometry.
+"""
+
+from yapcad.geom import *
+from yapcad.geom3d import *
+from yapcad.geom3d_util import *
+
+def main():
+    print("yapCAD Solid Topology Analysis Demo")
+    print("=" * 50)
+    print()
+
+    # Example 1: Simple cube
+    print("Example 1: Unit Cube")
+    print("-" * 30)
+    cube = prism(1, 1, 1)
+    is_closed = issolidclosed(cube)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(cube)
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: 1.0")
+    print()
+
+    # Example 2: 2x3x4 Rectangular prism
+    print("Example 2: Rectangular Prism (2×3×4)")
+    print("-" * 30)
+    box = prism(2, 3, 4)
+    is_closed = issolidclosed(box)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(box)
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: 24.0")
+    print()
+
+    # Example 3: Sphere
+    print("Example 3: Sphere (radius 2.0)")
+    print("-" * 30)
+    import math
+    radius = 2.0
+    sph = sphere(2 * radius, depth=3)
+    is_closed = issolidclosed(sph)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(sph)
+        expected = (4.0/3.0) * math.pi * (radius ** 3)
+        error = abs(vol - expected) / expected * 100
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: {expected:.6f}")
+        print(f"  Error: {error:.2f}%")
+    print()
+
+    # Example 4: Extruded square with hole
+    print("Example 4: Extruded Square with Hole")
+    print("-" * 30)
+    outer = [point(0, 0), point(4, 0), point(4, 4), point(0, 4), point(0, 0)]
+    hole = [point(1, 1), point(3, 1), point(3, 3), point(1, 3), point(1, 1)]
+    surf, _ = poly2surfaceXY(outer, holepolys=[hole])
+    extruded = extrude(surf, 2.0)
+    is_closed = issolidclosed(extruded)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(extruded)
+        # (4*4 - 2*2) * 2 = 24
+        expected = (4 * 4 - 2 * 2) * 2
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: {expected:.1f}")
+    print()
+
+    # Example 5: Cylinder
+    print("Example 5: Cylinder (radius 2, height 5)")
+    print("-" * 30)
+    radius = 2.0
+    height = 5.0
+    cylinder = conic(radius, radius, height)
+    is_closed = issolidclosed(cylinder)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(cylinder)
+        expected = math.pi * (radius ** 2) * height
+        error = abs(vol - expected) / expected * 100
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: {expected:.6f}")
+        print(f"  Error: {error:.2f}%")
+    print()
+
+    # Example 6: Tetrahedron
+    print("Example 6: Regular Tetrahedron")
+    print("-" * 30)
+    # Create tetrahedron vertices
+    tet_verts = [
+        point(1, 1, 1),
+        point(-1, -1, 1),
+        point(-1, 1, -1),
+        point(1, -1, -1)
+    ]
+
+    # Create faces
+    faces = [
+        [0, 2, 1],
+        [0, 1, 3],
+        [1, 2, 3],
+        [2, 0, 3]
+    ]
+
+    # Create surface
+    all_verts = []
+    all_norms = []
+    for f_idx, f in enumerate(faces):
+        tri = [tet_verts[f[0]], tet_verts[f[1]], tet_verts[f[2]]]
+        p0, n = tri2p0n(tri)
+        for v_idx in f:
+            all_verts.append(tet_verts[v_idx])
+            all_norms.append(n)
+
+    new_faces = []
+    for i in range(len(faces)):
+        new_faces.append([i * 3, i * 3 + 1, i * 3 + 2])
+
+    surf_tet = surface(all_verts, all_norms, new_faces)
+    tet = solid([surf_tet])
+
+    is_closed = issolidclosed(tet)
+    print(f"  Is closed: {is_closed}")
+    if is_closed:
+        vol = volumeof(tet)
+        # Calculate edge length and expected volume: V = a³/(6√2)
+        edge = line(tet_verts[0], tet_verts[1])
+        a = length(edge)
+        expected = (a ** 3) / (6 * math.sqrt(2))
+        error = abs(vol - expected) / expected * 100
+        print(f"  Edge length: {a:.6f}")
+        print(f"  Volume: {vol:.6f}")
+        print(f"  Expected: {expected:.6f}")
+        print(f"  Error: {error:.4f}%")
+    print()
+
+    # Example 7: Open box (missing top) - should NOT be closed
+    print("Example 7: Open Box (missing top)")
+    print("-" * 30)
+    bottom = rectangularPlane(2, 2, point(0, 0, 0))
+
+    side1 = rectangularPlane(2, 1, point(0, 0, 0.5))
+    side1 = rotatesurface(side1, 90, axis=point(1, 0, 0))
+    side1 = translatesurface(side1, point(0, -1, 0))
+
+    side2 = rectangularPlane(2, 1, point(0, 0, 0.5))
+    side2 = rotatesurface(side2, 90, axis=point(1, 0, 0))
+    side2 = translatesurface(side2, point(0, 1, 0))
+
+    side3 = rectangularPlane(2, 1, point(0, 0, 0.5))
+    side3 = rotatesurface(side3, 90, axis=point(0, 1, 0))
+    side3 = translatesurface(side3, point(-1, 0, 0))
+
+    side4 = rectangularPlane(2, 1, point(0, 0, 0.5))
+    side4 = rotatesurface(side4, 90, axis=point(0, 1, 0))
+    side4 = translatesurface(side4, point(1, 0, 0))
+
+    open_box = solid([bottom, side1, side2, side3, side4])
+    is_closed = issolidclosed(open_box)
+    print(f"  Is closed: {is_closed}")
+    print(f"  (Expected: False - box has open top)")
+    if is_closed:
+        vol = volumeof(open_box)
+        print(f"  Volume: {vol:.6f}")
+    else:
+        print(f"  Cannot compute volume - solid not closed")
+    print()
+
+    print("Demo complete!")
+    print("=" * 50)
+
+if __name__ == "__main__":
+    main()
diff --git a/src/yapcad/geom3d.py b/src/yapcad/geom3d.py
index aeb348a..38cf5ba 100644
--- a/src/yapcad/geom3d.py
+++ b/src/yapcad/geom3d.py
@@ -125,7 +125,7 @@
 
 ``solid = ['solid', surfaces, material, construction ]``, where:
 
-           ``surfaces`` is a list of surfaces with contiguous boundaries 
+           ``surfaces`` is a list of surfaces with contiguous boundaries
            that completely encloses an interior space,
 
            ``material`` is a list of domain-specific representation of
@@ -136,6 +136,30 @@
            ``construction`` is a list that contains information about
            how the solid was constructed, and may be empty
 
+Topology Analysis
+-----------------
+
+Two key functions are provided for analyzing solid topology:
+
+``issolidclosed(solid)`` -- Verifies that a solid is topologically closed
+by ensuring every edge is shared by exactly two faces across all surfaces.
+This is essential for determining if a solid properly encloses a volume
+without gaps or holes. Returns True if closed, False otherwise.
+
+``volumeof(solid)`` -- Calculates the volume enclosed by a closed solid
+using the divergence theorem. Requires that the solid be topologically
+closed (verified by calling ``issolidclosed()``). Returns the volume
+as a non-negative floating point number.
+
+Example usage::
+
+    from yapcad.geom3d_util import prism
+    from yapcad.geom3d import issolidclosed, volumeof
+
+    cube = prism(2, 2, 2)
+    if issolidclosed(cube):
+        vol = volumeof(cube)  # Returns 8.0
+
 
 Assembly
 --------
@@ -700,6 +724,171 @@ def mirrorsolid(x,plane,preserveNormal=True):
     s2[1] = surfs
     return s2
 
+def _point_to_key(p):
+    """
+    Convert a point to a hashable key for edge/vertex identification.
+    Uses rounded coordinates to handle floating point precision issues.
+    """
+    # Round to a reasonable precision to handle floating point comparison
+    return (round(p[0] / epsilon) * epsilon,
+            round(p[1] / epsilon) * epsilon,
+            round(p[2] / epsilon) * epsilon)
+
+def _canonical_edge_key(p1, p2):
+    """
+    Create a canonical edge key from two points.
+    Returns tuple of keys in sorted order so (p1,p2) and (p2,p1) map to same edge.
+    """
+    k1 = _point_to_key(p1)
+    k2 = _point_to_key(p2)
+    return (min(k1, k2), max(k1, k2))
+
+def issolidclosed(x):
+    """
+    Check if solid x is topologically closed.
+
+    A solid is closed if and only if every edge is shared by exactly two
+    faces across all surfaces. This ensures no holes or gaps exist in the
+    solid's boundary.
+
+    The function analyzes face adjacency by:
+    1. Building a global edge map using vertex positions (not indices)
+    2. Counting how many faces share each edge
+    3. Verifying that every edge is shared by exactly 2 faces
+
+    Args:
+        x: A solid data structure
+
+    Returns:
+        True if the solid is topologically closed, False otherwise
+
+    Raises:
+        ValueError: if x is not a valid solid
+
+    Example:
+        >>> from yapcad.geom3d_util import prism, sphere
+        >>> cube = prism(2, 2, 2)
+        >>> issolidclosed(cube)
+        True
+    """
+    # First verify this is a valid solid
+    if not issolid(x, fast=False):
+        raise ValueError('invalid solid passed to issolidclosed')
+
+    surfaces = x[1]
+
+    # Empty solid is trivially closed
+    if not surfaces:
+        return True
+
+    # Build a global edge map across all surfaces
+    # Key: canonical edge tuple (point_key1, point_key2)
+    # Value: count of faces that share this edge
+    global_edge_count = {}
+
+    for surf_idx, surf in enumerate(surfaces):
+        faces = surf[3]
+        vertices = surf[1]
+
+        # Process each face in this surface
+        for face_idx, face in enumerate(faces):
+            if len(face) != 3:
+                raise ValueError(f'non-triangular face in surface {surf_idx}, face {face_idx}')
+
+            # Get the three vertex positions
+            p0 = vertices[face[0]]
+            p1 = vertices[face[1]]
+            p2 = vertices[face[2]]
+
+            # Extract the three edges of this triangular face
+            edges = [
+                _canonical_edge_key(p0, p1),
+                _canonical_edge_key(p1, p2),
+                _canonical_edge_key(p2, p0)
+            ]
+
+            # Count this face's contribution to each edge
+            for edge in edges:
+                if edge not in global_edge_count:
+                    global_edge_count[edge] = 0
+                global_edge_count[edge] += 1
+
+    # Check that every edge is shared by exactly 2 faces
+    for edge, count in global_edge_count.items():
+        if count != 2:
+            return False
+
+    return True
+
+def volumeof(x):
+    """
+    Calculate the volume enclosed by a solid.
+
+    Uses the divergence theorem to compute volume from the surface triangulation.
+    For each triangular face with vertices (p0, p1, p2), the signed volume
+    contribution is: V_i = (1/6) * dot(p0, cross(p1-p0, p2-p0))
+
+    The total volume is the sum of absolute values of all face contributions.
+
+    Args:
+        x: A solid data structure (must be topologically closed)
+
+    Returns:
+        float: The volume of the solid (always non-negative)
+
+    Raises:
+        ValueError: if x is not a valid solid or is not closed
+
+    Example:
+        >>> from yapcad.geom3d_util import prism
+        >>> cube = prism(2, 2, 2)
+        >>> abs(volumeof(cube) - 8.0) < 0.001
+        True
+    """
+    # Verify this is a valid, closed solid
+    if not issolid(x, fast=False):
+        raise ValueError('invalid solid passed to volumeof')
+
+    if not issolidclosed(x):
+        raise ValueError('solid must be topologically closed to compute volume')
+
+    surfaces = x[1]
+
+    # Handle empty solid
+    if not surfaces:
+        return 0.0
+
+    total_volume = 0.0
+
+    # Accumulate signed volume contributions from all faces
+    for surf in surfaces:
+        vertices = surf[1]
+        faces = surf[3]
+
+        for face in faces:
+            if len(face) != 3:
+                raise ValueError('non-triangular face encountered')
+
+            # Get the three vertices of this face
+            p0 = vertices[face[0]]
+            p1 = vertices[face[1]]
+            p2 = vertices[face[2]]
+
+            # Compute vectors from p0 to other vertices
+            v1 = sub(p1, p0)
+            v2 = sub(p2, p0)
+
+            # Signed volume contribution: (1/6) * dot(p0, cross(v1, v2))
+            # This is the volume of the tetrahedron formed by the origin
+            # and the three face vertices
+            cross_product = cross(v1, v2)
+            signed_volume = dot(p0, cross_product) / 6.0
+
+            total_volume += signed_volume
+
+    # Return absolute value (orientation might cause negative result)
+    return abs(total_volume)
+
 def normfunc(tri):
     """
     utility funtion to compute normals for a flat facet triangle
diff --git a/tests/test_geom3d.py b/tests/test_geom3d.py
index c52ad97..ad62d65 100644
--- a/tests/test_geom3d.py
+++ b/tests/test_geom3d.py
@@ -708,5 +708,298 @@ def fcontour(x,ir=10,len=20,fr=5):
         airship2 = mirror(airship2,"yz")
         airship2 = translate(airship2,point(0,0,50))
         dd.draw(airship2,name='airship')
-        
+
         dd.display()
+
+
+class TestSolidTopology:
+    """Test solid topology analysis functions: issolidclosed and volumeof"""
+
+    def test_issolidclosed_empty_solid(self):
+        """Empty solid should be considered closed"""
+        empty = solid([])
+        assert issolidclosed(empty)
+
+    def test_issolidclosed_single_surface(self):
+        """Single open surface is not a closed solid"""
+        # Create a simple triangular surface (not closed)
+        p1 = point(0, 0, 0)
+        p2 = point(1, 0, 0)
+        p3 = point(0.5, 1, 0)
+        surf = surface([p1, p2, p3],
+                      [[0, 0, 1, 0]] * 3,
+                      [[0, 1, 2]])
+
+        s = solid([surf])
+        # Single surface can't be closed - edges not shared by 2 faces
+        assert not issolidclosed(s)
+
+    def test_issolidclosed_tetrahedron(self):
+        """Tetrahedron is a closed solid"""
+        # Create tetrahedron vertices
+        tet_verts = [
+            point(1, 1, 1),
+            point(-1, -1, 1),
+            point(-1, 1, -1),
+            point(1, -1, -1)
+        ]
+
+        # Create four triangular faces with correct winding
+        faces = [
+            [0, 2, 1],  # face pointing outward
+            [0, 1, 3],
+            [1, 2, 3],
+            [2, 0, 3]
+        ]
+
+        # Compute normals for each vertex of each face
+        normals = []
+        for f in faces:
+            tri = [tet_verts[f[0]], tet_verts[f[1]], tet_verts[f[2]]]
+            p0, n = tri2p0n(tri)
+            normals.extend([n, n, n])
+
+        # Create vertices list - each vertex appears in multiple faces
+        all_verts = []
+        all_norms = []
+        for f_idx, f in enumerate(faces):
+            for v_idx in f:
+                all_verts.append(tet_verts[v_idx])
+                all_norms.append(normals[f_idx * 3])
+
+        # Remap face indices to the expanded vertex list
+        new_faces = []
+        for i in range(len(faces)):
+            new_faces.append([i * 3, i * 3 + 1, i * 3 + 2])
+
+        surf = surface(all_verts, all_norms, new_faces)
+        tet = solid([surf])
+
+        assert issolidclosed(tet)
+
+    def test_issolidclosed_cube(self):
+        """Cube created by prism should be closed"""
+        from yapcad.geom3d_util import prism
+        cube = prism(2, 2, 2)
+        assert issolidclosed(cube)
+
+    def test_issolidclosed_sphere(self):
+        """Sphere should be closed"""
+        from yapcad.geom3d_util import sphere
+        sph = sphere(2.0, depth=2)
+        assert issolidclosed(sph)
+
+    def test_issolidclosed_extruded_square(self):
+        """Extruded square should be closed"""
+        from yapcad.geom3d_util import extrude
+        # Create square polygon
+        square = [point(0, 0), point(2, 0), point(2, 2), point(0, 2), point(0, 0)]
+        surf, _ = poly2surfaceXY(square)
+        extruded = extrude(surf, 3.0)
+        assert issolidclosed(extruded)
+
+    def test_issolidclosed_extruded_with_hole(self):
+        """Extruded surface with hole should be closed"""
+        from yapcad.geom3d_util import extrude
+        # Outer square
+        outer = [point(0, 0), point(4, 0), point(4, 4), point(0, 4), point(0, 0)]
+        # Inner square (hole)
+        hole = [point(1, 1), point(3, 1), point(3, 3), point(1, 3), point(1, 1)]
+
+        surf, _ = poly2surfaceXY(outer, holepolys=[hole])
+        extruded = extrude(surf, 2.0)
+        assert issolidclosed(extruded)
+
+    def test_issolidclosed_open_box(self):
+        """Box missing a face should not be closed"""
+        from yapcad.geom3d_util import rectangularPlane, extrude
+
+        # Create bottom surface
+        bottom = rectangularPlane(2, 2, point(0, 0, 0))
+
+        # Create four side surfaces manually (missing top)
+        side1 = rectangularPlane(2, 1, point(0, 0, 0.5))
+        side1 = rotatesurface(side1, 90, axis=point(1, 0, 0))
+        side1 = translatesurface(side1, point(0, -1, 0))
+
+        side2 = rectangularPlane(2, 1, point(0, 0, 0.5))
+        side2 = rotatesurface(side2, 90, axis=point(1, 0, 0))
+        side2 = translatesurface(side2, point(0, 1, 0))
+
+        side3 = rectangularPlane(2, 1, point(0, 0, 0.5))
+        side3 = rotatesurface(side3, 90, axis=point(0, 1, 0))
+        side3 = translatesurface(side3, point(-1, 0, 0))
+
+        side4 = rectangularPlane(2, 1, point(0, 0, 0.5))
+        side4 = rotatesurface(side4, 90, axis=point(0, 1, 0))
+        side4 = translatesurface(side4, point(1, 0, 0))
+
+        # Create open box (no top)
+        open_box = solid([bottom, side1, side2, side3, side4])
+
+        # Should not be closed
+        assert not issolidclosed(open_box)
+
+    def test_volumeof_unit_cube(self):
+        """Volume of unit cube should be 1"""
+        from yapcad.geom3d_util import prism
+        cube = prism(1, 1, 1)
+        vol = volumeof(cube)
+        assert abs(vol - 1.0) < 0.001
+
+    def test_volumeof_cube_2x2x2(self):
+        """Volume of 2x2x2 cube should be 8"""
+        from yapcad.geom3d_util import prism
+        cube = prism(2, 2, 2)
+        vol = volumeof(cube)
+        assert abs(vol - 8.0) < 0.001
+
+    def test_volumeof_rectangular_prism(self):
+        """Volume of 3x4x5 prism should be 60"""
+        from yapcad.geom3d_util import prism
+        box = prism(3, 4, 5)
+        vol = volumeof(box)
+        assert abs(vol - 60.0) < 0.01
+
+    def test_volumeof_sphere(self):
+        """Volume of sphere should match 4πr³/3"""
+        from yapcad.geom3d_util import sphere
+        import math
+
+        radius = 2.0
+        expected = (4.0 / 3.0) * math.pi * (radius ** 3)
+
+        sph = sphere(2 * radius, depth=3)  # depth=3 for better approximation
+        vol = volumeof(sph)
+
+        # Allow 1% error due to triangulation approximation
+        assert abs(vol - expected) / expected < 0.01
+
+    def test_volumeof_extruded_square(self):
+        """Volume of extruded 2x2 square by height 3 should be 12"""
+        from yapcad.geom3d_util import extrude
+
+        square = [point(0, 0), point(2, 0), point(2, 2), point(0, 2), point(0, 0)]
+        surf, _ = poly2surfaceXY(square)
+        extruded = extrude(surf, 3.0)
+
+        vol = volumeof(extruded)
+        assert abs(vol - 12.0) < 0.01
+
+    def test_volumeof_extruded_with_hole(self):
+        """Volume of extruded surface with hole"""
+        from yapcad.geom3d_util import extrude
+
+        # 4x4 outer square with 2x2 inner hole
+        # Expected volume: (4*4 - 2*2) * height = 12 * 2 = 24
+        outer = [point(0, 0), point(4, 0), point(4, 4), point(0, 4), point(0, 0)]
+        hole = [point(1, 1), point(3, 1), point(3, 3), point(1, 3), point(1, 1)]
+
+        surf, _ = poly2surfaceXY(outer, holepolys=[hole])
+        extruded = extrude(surf, 2.0)
+
+        vol = volumeof(extruded)
+        expected = (4 * 4 - 2 * 2) * 2
+        assert abs(vol - expected) < 0.1
+
+    def test_volumeof_truncated_cone(self):
+        """Volume of truncated cone (frustum) should match formula"""
+        from yapcad.geom3d_util import conic
+        import math
+
+        base_radius = 3.0
+        top_radius = 1.0
+        height = 4.0
+        # Volume of frustum: (π*h/3)(r1² + r1*r2 + r2²)
+        expected = (math.pi * height / 3.0) * (
+            base_radius**2 + base_radius*top_radius + top_radius**2
+        )
+
+        frustum = conic(base_radius, top_radius, height, angr=10)
+        vol = volumeof(frustum)
+
+        # Allow 2% error due to triangulation
+        assert abs(vol - expected) / expected < 0.02
+
+    def test_volumeof_cylinder(self):
+        """Volume of cylinder should match πr²h"""
+        from yapcad.geom3d_util import conic
+        import math
+
+        radius = 2.0
+        height = 5.0
+        expected = math.pi * (radius ** 2) * height
+
+        cylinder = conic(radius, radius, height, angr=5)
+        vol = volumeof(cylinder)
+
+        # Allow 1% error
+        assert abs(vol - expected) / expected < 0.01
+
+    def test_volumeof_requires_closed_solid(self):
+        """volumeof should raise error for non-closed solid"""
+        # Create a single triangular surface (not closed)
+        p1 = point(0, 0, 0)
+        p2 = point(1, 0, 0)
+        p3 = point(0.5, 1, 0)
+        surf = surface([p1, p2, p3],
+                      [[0, 0, 1, 0]] * 3,
+                      [[0, 1, 2]])
+
+        s = solid([surf])
+
+        with pytest.raises(ValueError, match="must be topologically closed"):
+            volumeof(s)
+
+    def test_volumeof_tetrahedron(self):
+        """Volume of tetrahedron should match a³/(6√2)"""
+        import math
+
+        # Create a tetrahedron with vertices at specific positions
+        # This forms a tetrahedron with all edges of equal length
+        tet_verts = [
+            point(1, 1, 1),
+            point(-1, -1, 1),
+            point(-1, 1, -1),
+            point(1, -1, -1)
+        ]
+
+        # Create faces
+        faces = [
+            [0, 2, 1],
+            [0, 1, 3],
+            [1, 2, 3],
+            [2, 0, 3]
+        ]
+
+        # Create surface
+        all_verts = []
+        all_norms = []
+        for f_idx, f in enumerate(faces):
+            tri = [tet_verts[f[0]], tet_verts[f[1]], tet_verts[f[2]]]
+            p0, n = tri2p0n(tri)
+            for v_idx in f:
+                all_verts.append(tet_verts[v_idx])
+                all_norms.append(n)
+
+        new_faces = []
+        for i in range(len(faces)):
+            new_faces.append([i * 3, i * 3 + 1, i * 3 + 2])
+
+        surf = surface(all_verts, all_norms, new_faces)
+        tet = solid([surf])
+
+        # Verify it's closed
+        assert issolidclosed(tet)
+
+        # Calculate actual edge length using the length() function
+        edge = line(tet_verts[0], tet_verts[1])
+        a = length(edge)
+
+        # Volume of regular tetrahedron: V = a³/(6√2)
+        expected = (a ** 3) / (6 * math.sqrt(2))
+
+        vol = volumeof(tet)
+
+        # Check volume matches expected within small tolerance
+        assert abs(vol - expected) < 0.001