diff --git a/.gitignore b/.gitignore
index c35eb05..721a58e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,6 +22,8 @@ test_*
 get*
 twot*
 regex*
+run_tests
+run_connectivity_tests
 
 # CodeQL build artifacts
 _codeql_build_dir
diff --git a/src/YgorMathContourConnectivity.cc b/src/YgorMathContourConnectivity.cc
new file mode 100644
index 0000000..36e7648
--- /dev/null
+++ b/src/YgorMathContourConnectivity.cc
@@ -0,0 +1,409 @@
+//YgorMathContourConnectivity.cc
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <list>
+#include <map>
+#include <set>
+#include <stdexcept>
+#include <vector>
+#include <cstdint>
+
+#include "YgorMath.h"
+#include "YgorMathContourConnectivity.h"
+
+// Helper structure for 2D projected vertex with source information.
+template <class T>
+struct ProjectedVertex {
+    vec2<T> pos_2d;          // 2D position in the projection plane
+    vec3<T> pos_3d;          // Original 3D position
+    size_t vertex_index;     // Index in the output mesh vertices array
+    bool is_top;             // True if from top collection, false if from bottom
+    size_t contour_index;    // Index of the contour within the collection
+    size_t point_index;      // Index of the point within the contour
+};
+
+template <class T>
+fv_surface_mesh<T, uint64_t>
+Triangulate_Planar_Contour_Connectivity(
+    const contour_collection<T> &top_cc,
+    const contour_collection<T> &bottom_cc,
+    T distance_eps
+){
+    using IndexType = uint64_t;
+    
+    fv_surface_mesh<T, IndexType> mesh;
+    
+    // Validate inputs
+    if (top_cc.contours.empty()) {
+        throw std::invalid_argument("Top contour collection is empty.");
+    }
+    if (bottom_cc.contours.empty()) {
+        throw std::invalid_argument("Bottom contour collection is empty.");
+    }
+    
+    // Compute average points using existing contour_collection API.
+    // This handles the empty/NaN cases appropriately.
+    const vec3<T> avg_top_point = top_cc.Average_Point();
+    const vec3<T> avg_bottom_point = bottom_cc.Average_Point();
+    
+    // Check if average points are valid (collections may be non-empty but contain no actual vertices)
+    if (!avg_top_point.isfinite() || !avg_bottom_point.isfinite()) {
+        throw std::invalid_argument("Contour collections contain no vertices.");
+    }
+    
+    // Compute the plane separation direction (from bottom to top)
+    vec3<T> separation = avg_top_point - avg_bottom_point;
+    const T separation_dist = separation.length();
+    
+    if (separation_dist < distance_eps) {
+        throw std::runtime_error("Top and bottom planes appear to be coincident.");
+    }
+    
+    // Plane normal is the direction perpendicular to both planes
+    vec3<T> plane_normal = separation.unit();
+    
+    // Create orthonormal basis for 2D projection.
+    // U and V will span the projection plane, N is the normal.
+    vec3<T> U, V;
+    
+    // Find a vector not parallel to plane_normal
+    vec3<T> arbitrary(static_cast<T>(1), static_cast<T>(0), static_cast<T>(0));
+    if (std::abs(plane_normal.Dot(arbitrary)) > static_cast<T>(0.9)) {
+        arbitrary = vec3<T>(static_cast<T>(0), static_cast<T>(1), static_cast<T>(0));
+    }
+    
+    U = plane_normal.Cross(arbitrary).unit();
+    V = plane_normal.Cross(U).unit();
+    
+    // Projection origin (midpoint between average points)
+    vec3<T> proj_origin = (avg_top_point + avg_bottom_point) / static_cast<T>(2);
+    
+    // Lambda to project 3D point to 2D
+    auto project_to_2d = [&](const vec3<T> &pt) -> vec2<T> {
+        vec3<T> rel = pt - proj_origin;
+        return vec2<T>(rel.Dot(U), rel.Dot(V));
+    };
+    
+    // Collect all projected vertices
+    std::vector<ProjectedVertex<T>> projected_verts;
+    
+    // Add all vertices from both collections to the mesh and create projections
+    size_t contour_idx = 0;
+    for (const auto &contour : top_cc.contours) {
+        size_t pt_idx = 0;
+        for (const auto &pt : contour.points) {
+            ProjectedVertex<T> pv;
+            pv.pos_3d = pt;
+            pv.pos_2d = project_to_2d(pt);
+            pv.vertex_index = mesh.vertices.size();
+            pv.is_top = true;
+            pv.contour_index = contour_idx;
+            pv.point_index = pt_idx;
+            
+            mesh.vertices.push_back(pt);
+            projected_verts.push_back(pv);
+            ++pt_idx;
+        }
+        ++contour_idx;
+    }
+    
+    contour_idx = 0;
+    for (const auto &contour : bottom_cc.contours) {
+        size_t pt_idx = 0;
+        for (const auto &pt : contour.points) {
+            ProjectedVertex<T> pv;
+            pv.pos_3d = pt;
+            pv.pos_2d = project_to_2d(pt);
+            pv.vertex_index = mesh.vertices.size();
+            pv.is_top = false;
+            pv.contour_index = contour_idx;
+            pv.point_index = pt_idx;
+            
+            mesh.vertices.push_back(pt);
+            projected_verts.push_back(pv);
+            ++pt_idx;
+        }
+        ++contour_idx;
+    }
+    
+    // Compute total number of top vertices for indexing bottom vertices later.
+    size_t top_vertex_count = 0;
+    for (const auto &contour : top_cc.contours) {
+        top_vertex_count += contour.points.size();
+    }
+    
+    // Now we need to triangulate the gap between top and bottom contours.
+    // We use a greedy approach that connects vertices from both planes:
+    // 1. For each contour edge on one plane, find nearby vertices on the other plane
+    // 2. Create triangles connecting them
+    
+    // Sort projected vertices by angle around centroid for better connectivity
+    vec2<T> centroid_2d(static_cast<T>(0), static_cast<T>(0));
+    for (const auto &pv : projected_verts) {
+        centroid_2d += pv.pos_2d;
+    }
+    centroid_2d /= static_cast<T>(projected_verts.size());
+    
+    // Helper to compute angle from centroid
+    auto angle_from_centroid = [&centroid_2d](const vec2<T> &pos) -> T {
+        vec2<T> delta = pos - centroid_2d;
+        return std::atan2(delta.y, delta.x);
+    };
+    
+    // Create lists of top and bottom vertex indices sorted by angle
+    std::vector<size_t> top_indices, bottom_indices;
+    for (size_t i = 0; i < projected_verts.size(); ++i) {
+        if (projected_verts[i].is_top) {
+            top_indices.push_back(i);
+        } else {
+            bottom_indices.push_back(i);
+        }
+    }
+    
+    // Sort by angle
+    auto sort_by_angle = [&](std::vector<size_t> &indices) {
+        std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) {
+            return angle_from_centroid(projected_verts[a].pos_2d) < 
+                   angle_from_centroid(projected_verts[b].pos_2d);
+        });
+    };
+    
+    sort_by_angle(top_indices);
+    sort_by_angle(bottom_indices);
+    
+    // Use a sweep approach to connect vertices.
+    // For each vertex on the top, find the closest vertex on the bottom and vice versa.
+    // Then create triangles.
+    
+    auto find_closest_on_other_plane = [&](size_t idx, bool from_top) -> size_t {
+        const auto &pos = projected_verts[idx].pos_2d;
+        const auto &other_indices = from_top ? bottom_indices : top_indices;
+        
+        if (other_indices.empty()) {
+            throw std::runtime_error("Cannot find closest vertex: other plane has no vertices.");
+        }
+        
+        size_t closest = other_indices[0];
+        T min_dist = (projected_verts[closest].pos_2d - pos).sq_length();
+        
+        for (size_t i = 1; i < other_indices.size(); ++i) {
+            T dist = (projected_verts[other_indices[i]].pos_2d - pos).sq_length();
+            if (dist < min_dist) {
+                min_dist = dist;
+                closest = other_indices[i];
+            }
+        }
+        return closest;
+    };
+    
+    // Helper function to check if two 2D line segments intersect (excluding shared endpoints).
+    // Uses the cross-product orientation test.
+    auto segments_intersect_2d = [&](const vec2<T> &a1, const vec2<T> &a2, 
+                                     const vec2<T> &b1, const vec2<T> &b2) -> bool {
+        // Cross product of 2D vectors (returns scalar z-component)
+        auto cross_2d = [](const vec2<T> &u, const vec2<T> &v) -> T {
+            return u.x * v.y - u.y * v.x;
+        };
+        
+        const T d1 = cross_2d(b2 - b1, a1 - b1);
+        const T d2 = cross_2d(b2 - b1, a2 - b1);
+        const T d3 = cross_2d(a2 - a1, b1 - a1);
+        const T d4 = cross_2d(a2 - a1, b2 - a1);
+        
+        // Segments straddle each other if they have opposite orientations
+        const T eps = distance_eps;
+        if (((d1 > eps && d2 < -eps) || (d1 < -eps && d2 > eps)) &&
+            ((d3 > eps && d4 < -eps) || (d3 < -eps && d4 > eps))) {
+            return true;
+        }
+        return false;
+    };
+    
+    // Helper to check if triangle is valid:
+    // 1. Non-degenerate (has non-zero 3D area)
+    // 2. Consistent orientation (normal has a meaningful direction)
+    // 3. Edges do not intersect existing face edges in 2D projection
+    
+    // Create a map from vertex_index to projected 2D position for O(1) lookups
+    std::map<IndexType, vec2<T>> vertex_to_2d;
+    for (const auto &pv : projected_verts) {
+        vertex_to_2d[static_cast<IndexType>(pv.vertex_index)] = pv.pos_2d;
+    }
+    
+    auto is_valid_triangle = [&](size_t v0, size_t v1, size_t v2) -> bool {
+        const auto &p0_3d = projected_verts[v0].pos_3d;
+        const auto &p1_3d = projected_verts[v1].pos_3d;
+        const auto &p2_3d = projected_verts[v2].pos_3d;
+        
+        // Calculate 3D triangle area using cross product
+        const vec3<T> edge1 = p1_3d - p0_3d;
+        const vec3<T> edge2 = p2_3d - p0_3d;
+        const vec3<T> tri_normal = edge1.Cross(edge2);
+        const T area = tri_normal.length() / static_cast<T>(2);
+        
+        // Reject degenerate triangles
+        if (area < distance_eps * distance_eps) {
+            return false;
+        }
+        
+        // Orientation check: for bridging triangles connecting top and bottom planes,
+        // we verify the triangle involves vertices from both planes. Triangles entirely
+        // on one plane would be degenerate in 3D (caught by the area check above).
+        // The non-degenerate area check is sufficient for valid bridging triangles.
+        
+        // Check for edge intersections with existing faces in 2D projection
+        const auto &p0_2d = projected_verts[v0].pos_2d;
+        const auto &p1_2d = projected_verts[v1].pos_2d;
+        const auto &p2_2d = projected_verts[v2].pos_2d;
+        
+        // Define edges of the new triangle
+        std::array<std::pair<vec2<T>, vec2<T>>, 3> new_edges = {{
+            {p0_2d, p1_2d},
+            {p1_2d, p2_2d},
+            {p2_2d, p0_2d}
+        }};
+        
+        // Check against all edges of existing faces
+        for (const auto &face : mesh.faces) {
+            if (face.size() < 3) continue;
+            
+            // Get 2D positions for existing face vertices using O(1) map lookup
+            std::vector<vec2<T>> face_verts_2d;
+            for (const auto &idx : face) {
+                auto it = vertex_to_2d.find(idx);
+                if (it != vertex_to_2d.end()) {
+                    face_verts_2d.push_back(it->second);
+                }
+            }
+            
+            if (face_verts_2d.size() < 3) continue;
+            
+            // Check each edge of the existing face
+            for (size_t i = 0; i < face_verts_2d.size(); ++i) {
+                size_t j = (i + 1) % face_verts_2d.size();
+                const auto &e1 = face_verts_2d[i];
+                const auto &e2 = face_verts_2d[j];
+                
+                // Check against each edge of the new triangle
+                for (const auto &new_edge : new_edges) {
+                    // Skip if edges share endpoints (they connect at a vertex)
+                    const T eps_sq = distance_eps * distance_eps;
+                    if ((new_edge.first - e1).sq_length() < eps_sq ||
+                        (new_edge.first - e2).sq_length() < eps_sq ||
+                        (new_edge.second - e1).sq_length() < eps_sq ||
+                        (new_edge.second - e2).sq_length() < eps_sq) {
+                        continue;
+                    }
+                    
+                    if (segments_intersect_2d(new_edge.first, new_edge.second, e1, e2)) {
+                        return false; // Edges intersect, reject triangle
+                    }
+                }
+            }
+        }
+        
+        return true;
+    };
+    
+    // Process each contour edge and try to form triangles with vertices on the other plane
+    std::set<std::tuple<size_t, size_t, size_t>> added_triangles;
+    
+    auto add_triangle_if_valid = [&](size_t v0, size_t v1, size_t v2) {
+        // Normalize triangle vertex order
+        std::array<size_t, 3> tri = {v0, v1, v2};
+        std::sort(tri.begin(), tri.end());
+        
+        auto key = std::make_tuple(tri[0], tri[1], tri[2]);
+        if (added_triangles.count(key) > 0) {
+            return; // Already added
+        }
+        
+        if (is_valid_triangle(v0, v1, v2)) {
+            added_triangles.insert(key);
+            mesh.faces.push_back({static_cast<IndexType>(projected_verts[v0].vertex_index),
+                                  static_cast<IndexType>(projected_verts[v1].vertex_index),
+                                  static_cast<IndexType>(projected_verts[v2].vertex_index)});
+        }
+    };
+    
+    // Triangulate by iterating through contour edges and connecting to nearest vertices on other plane
+    // Process top contour edges
+    size_t global_idx = 0;
+    for (const auto &contour : top_cc.contours) {
+        const size_t n = contour.points.size();
+        if (n < 2) {
+            global_idx += n;
+            continue;
+        }
+        
+        for (size_t i = 0; i < n; ++i) {
+            size_t next = (i + 1) % n;
+            if (!contour.closed && next == 0) continue;
+            
+            size_t ti1 = global_idx + i;
+            size_t ti2 = global_idx + next;
+            
+            // Find nearest bottom vertices to both top vertices
+            size_t bi1 = find_closest_on_other_plane(ti1, true);
+            size_t bi2 = find_closest_on_other_plane(ti2, true);
+            
+            // Create triangles
+            if (bi1 != bi2) {
+                // Two triangles forming a quad
+                add_triangle_if_valid(ti1, ti2, bi1);
+                add_triangle_if_valid(ti2, bi2, bi1);
+            } else {
+                // Single triangle
+                add_triangle_if_valid(ti1, ti2, bi1);
+            }
+        }
+        global_idx += n;
+    }
+    
+    // Process bottom contour edges
+    global_idx = top_vertex_count;
+    for (const auto &contour : bottom_cc.contours) {
+        const size_t n = contour.points.size();
+        if (n < 2) {
+            global_idx += n;
+            continue;
+        }
+        
+        for (size_t i = 0; i < n; ++i) {
+            size_t next = (i + 1) % n;
+            if (!contour.closed && next == 0) continue;
+            
+            size_t bi1 = global_idx + i;
+            size_t bi2 = global_idx + next;
+            
+            // Find nearest top vertices to both bottom vertices
+            size_t ti1 = find_closest_on_other_plane(bi1, false);
+            size_t ti2 = find_closest_on_other_plane(bi2, false);
+            
+            // Create triangles
+            if (ti1 != ti2) {
+                // Two triangles forming a quad
+                add_triangle_if_valid(bi1, bi2, ti1);
+                add_triangle_if_valid(bi2, ti2, ti1);
+            } else {
+                // Single triangle
+                add_triangle_if_valid(bi1, bi2, ti1);
+            }
+        }
+        global_idx += n;
+    }
+    
+    if (mesh.faces.empty()) {
+        throw std::runtime_error("Triangulation produced no valid triangles.");
+    }
+    
+    return mesh;
+}
+
+// Template instantiations
+#ifndef YGORMATH_DISABLE_ALL_SPECIALIZATIONS
+    template fv_surface_mesh<float,  uint64_t> Triangulate_Planar_Contour_Connectivity(const contour_collection<float>  &, const contour_collection<float>  &, float );
+    template fv_surface_mesh<double, uint64_t> Triangulate_Planar_Contour_Connectivity(const contour_collection<double> &, const contour_collection<double> &, double);
+#endif
diff --git a/src/YgorMathContourConnectivity.h b/src/YgorMathContourConnectivity.h
new file mode 100644
index 0000000..e5b6f55
--- /dev/null
+++ b/src/YgorMathContourConnectivity.h
@@ -0,0 +1,48 @@
+//YgorMathContourConnectivity.h
+
+#ifndef YGOR_MATH_CONTOUR_CONNECTIVITY_H_
+#define YGOR_MATH_CONTOUR_CONNECTIVITY_H_
+
+#include "YgorMath.h"
+
+//---------------------------------------------------------------------------------------------------------------------------
+//---------- contour connectivity: triangulate the gap between two planar contour collections on offset planes -------------
+//---------------------------------------------------------------------------------------------------------------------------
+//
+// This function computes the triangulation between two sets of planar contours that lie on parallel, offset planes.
+// The input consists of two contour_collections: 'top' and 'bottom', where all contours within each collection 
+// are coincident on a common plane. The planes are assumed to be parallel but separated by an 
+// arbitrary distance.
+//
+// The function projects all vertices onto a common 2D plane and then computes a triangulation that connects the
+// vertices of the top collection to the vertices of the bottom collection, creating a surface mesh that bridges
+// the gap between the two planes.
+//
+// Key features:
+// - Handles multiple contours within each collection (for example, disjoint polygons)
+// - Treats all input contours uniformly; callers should not rely on orientation to denote holes
+// - Does not assume the polygons are simple or monotone
+// - The output mesh is not watertight, but represents surface patches bridging the planes
+//
+// Parameters:
+// - top_cc: contour_collection with all contours lying on the 'top' plane
+// - bottom_cc: contour_collection with all contours lying on the 'bottom' plane
+// - distance_eps: scale for degeneracy/planarity tolerances (e.g., plane separation and small-triangle rejection) (default: 1E-6)
+//
+// Returns:
+// - fv_surface_mesh containing the triangulated surface bridging the two planes
+//
+// Throws:
+// - std::invalid_argument if either collection is empty
+// - std::runtime_error if triangulation fails
+//
+
+template <class T>
+fv_surface_mesh<T, uint64_t>
+Triangulate_Planar_Contour_Connectivity(
+    const contour_collection<T> &top_cc,
+    const contour_collection<T> &bottom_cc,
+    T distance_eps = static_cast<T>(1E-6)
+);
+
+#endif // YGOR_MATH_CONTOUR_CONNECTIVITY_H_
diff --git a/tests2/YgorMath/contour_connectivity.cc b/tests2/YgorMath/contour_connectivity.cc
new file mode 100644
index 0000000..56c14bd
--- /dev/null
+++ b/tests2/YgorMath/contour_connectivity.cc
@@ -0,0 +1,403 @@
+
+#include <limits>
+#include <list>
+#include <utility>
+#include <iostream>
+#include <cstdint>
+#include <cmath>
+
+#include <YgorMath.h>
+#include <YgorMathContourConnectivity.h>
+
+#include "doctest/doctest.h"
+
+
+TEST_CASE( "Triangulate_Planar_Contour_Connectivity basic tests" ){
+
+    SUBCASE("simple square contours on parallel planes"){
+        // Create a simple square contour on the top plane (z = 1)
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 1.0));
+
+        // Create a matching square contour on the bottom plane (z = 0)
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 0.0));
+
+        // Triangulate the gap
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        // Check that we have vertices (4 from top + 4 from bottom = 8)
+        REQUIRE(mesh.vertices.size() == 8);
+
+        // Check that we have triangles (should have triangles connecting the two squares)
+        REQUIRE(mesh.faces.size() > 0);
+
+        // All faces should be triangles
+        for (const auto &face : mesh.faces) {
+            REQUIRE(face.size() == 3);
+        }
+
+        // Surface area should be positive
+        REQUIRE(mesh.surface_area() > 0.0);
+    }
+
+    SUBCASE("triangular contours on parallel planes"){
+        // Create a simple triangle on the top plane
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 2.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(2.0, 0.0, 2.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.732, 2.0));
+
+        // Create a matching triangle on the bottom plane
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(2.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.732, 0.0));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        // Check vertices (3 + 3 = 6)
+        REQUIRE(mesh.vertices.size() == 6);
+
+        // Check that we have triangles
+        REQUIRE(mesh.faces.size() > 0);
+
+        // Surface area should be positive
+        REQUIRE(mesh.surface_area() > 0.0);
+    }
+
+    SUBCASE("empty top collection throws"){
+        contour_collection<double> top_cc;
+        
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+
+        REQUIRE_THROWS_AS(
+            Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc),
+            std::invalid_argument
+        );
+    }
+
+    SUBCASE("empty bottom collection throws"){
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        
+        contour_collection<double> bottom_cc;
+
+        REQUIRE_THROWS_AS(
+            Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc),
+            std::invalid_argument
+        );
+    }
+
+    SUBCASE("coincident planes throws"){
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+
+        REQUIRE_THROWS_AS(
+            Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc),
+            std::runtime_error
+        );
+    }
+
+    SUBCASE("collections with empty contour_of_points entries throw"){
+        // Create collections that have contours but no actual points
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();  // Empty contour_of_points
+        top_cc.contours.back().closed = true;
+        // No points added
+        
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+
+        REQUIRE_THROWS_AS(
+            Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc),
+            std::invalid_argument
+        );
+    }
+
+    SUBCASE("multiple contours with some empty entries throw"){
+        // Create a collection with one valid contour and one empty
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();  // First contour is empty
+        top_cc.contours.back().closed = true;
+        top_cc.contours.emplace_back();  // Second contour also empty
+        top_cc.contours.back().closed = true;
+        // No points in any contour
+        
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+
+        REQUIRE_THROWS_AS(
+            Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc),
+            std::invalid_argument
+        );
+    }
+}
+
+TEST_CASE( "Triangulate_Planar_Contour_Connectivity with multiple contours" ){
+    
+    SUBCASE("two separate squares on each plane (separate contours)"){
+        // Top plane: two separate squares
+        contour_collection<double> top_cc;
+        
+        // First square
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 1.0));
+        
+        // Second square (offset)
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(2.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(3.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(3.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(2.0, 1.0, 1.0));
+
+        // Bottom plane: matching squares
+        contour_collection<double> bottom_cc;
+        
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 0.0));
+        
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(2.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(3.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(3.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(2.0, 1.0, 0.0));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        // 8 vertices on each plane = 16 total
+        REQUIRE(mesh.vertices.size() == 16);
+        
+        // Should have triangles
+        REQUIRE(mesh.faces.size() > 0);
+        
+        // Surface area should be positive
+        REQUIRE(mesh.surface_area() > 0.0);
+    }
+
+    SUBCASE("polygon with hole pattern (outer and inner contour with opposite orientations)"){
+        // Create an outer square on top (counter-clockwise)
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(4.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(4.0, 4.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 4.0, 1.0));
+        
+        // Create an inner square (clockwise - denotes hole)
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 3.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(3.0, 3.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(3.0, 1.0, 1.0));
+
+        // Bottom plane with matching contours
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(4.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(4.0, 4.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 4.0, 0.0));
+        
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 3.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(3.0, 3.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(3.0, 1.0, 0.0));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        // 8 vertices on each plane = 16 total
+        REQUIRE(mesh.vertices.size() == 16);
+        
+        // Should have triangles
+        REQUIRE(mesh.faces.size() > 0);
+        
+        // Surface area should be positive
+        REQUIRE(mesh.surface_area() > 0.0);
+    }
+}
+
+TEST_CASE( "Triangulate_Planar_Contour_Connectivity with oblique planes" ){
+    
+    SUBCASE("planes not aligned with cardinal axes"){
+        // Create contours on a plane tilted at 45 degrees
+        const double pi = std::acos(-1.0);
+        const double angle = pi / 4.0; // 45 degrees
+        
+        // Top plane points (rotated around Y axis)
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        
+        // Original square at z=1, rotated 45 degrees around Y
+        auto rotate_y = [angle](const vec3<double> &v) {
+            return vec3<double>(
+                v.x * std::cos(angle) + v.z * std::sin(angle),
+                v.y,
+                -v.x * std::sin(angle) + v.z * std::cos(angle)
+            );
+        };
+        
+        vec3<double> offset_top(0.0, 0.0, 1.0);
+        top_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(0.0, 0.0, 0.0) + offset_top));
+        top_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(1.0, 0.0, 0.0) + offset_top));
+        top_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(1.0, 1.0, 0.0) + offset_top));
+        top_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(0.0, 1.0, 0.0) + offset_top));
+
+        // Bottom plane points
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        
+        vec3<double> offset_bottom(0.0, 0.0, 0.0);
+        bottom_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(0.0, 0.0, 0.0) + offset_bottom));
+        bottom_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(1.0, 0.0, 0.0) + offset_bottom));
+        bottom_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(1.0, 1.0, 0.0) + offset_bottom));
+        bottom_cc.contours.back().points.emplace_back(rotate_y(vec3<double>(0.0, 1.0, 0.0) + offset_bottom));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        // Check vertices
+        REQUIRE(mesh.vertices.size() == 8);
+        
+        // Should have triangles
+        REQUIRE(mesh.faces.size() > 0);
+        
+        // Surface area should be positive
+        REQUIRE(mesh.surface_area() > 0.0);
+    }
+}
+
+TEST_CASE( "Triangulate_Planar_Contour_Connectivity with float type" ){
+    
+    SUBCASE("basic square contours with float type"){
+        contour_collection<float> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<float>(0.0f, 0.0f, 1.0f));
+        top_cc.contours.back().points.emplace_back(vec3<float>(1.0f, 0.0f, 1.0f));
+        top_cc.contours.back().points.emplace_back(vec3<float>(1.0f, 1.0f, 1.0f));
+        top_cc.contours.back().points.emplace_back(vec3<float>(0.0f, 1.0f, 1.0f));
+
+        contour_collection<float> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<float>(0.0f, 0.0f, 0.0f));
+        bottom_cc.contours.back().points.emplace_back(vec3<float>(1.0f, 0.0f, 0.0f));
+        bottom_cc.contours.back().points.emplace_back(vec3<float>(1.0f, 1.0f, 0.0f));
+        bottom_cc.contours.back().points.emplace_back(vec3<float>(0.0f, 1.0f, 0.0f));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        REQUIRE(mesh.vertices.size() == 8);
+        REQUIRE(mesh.faces.size() > 0);
+        REQUIRE(mesh.surface_area() > 0.0f);
+    }
+}
+
+TEST_CASE( "Triangulate_Planar_Contour_Connectivity mesh integrity" ){
+    
+    SUBCASE("all face indices are valid"){
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 1.0));
+
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 0.0));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        const size_t n_verts = mesh.vertices.size();
+        for (const auto &face : mesh.faces) {
+            for (const auto &idx : face) {
+                REQUIRE(idx < n_verts);
+            }
+        }
+    }
+
+    SUBCASE("all triangles have non-zero area"){
+        contour_collection<double> top_cc;
+        top_cc.contours.emplace_back();
+        top_cc.contours.back().closed = true;
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 1.0));
+        top_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 1.0));
+
+        contour_collection<double> bottom_cc;
+        bottom_cc.contours.emplace_back();
+        bottom_cc.contours.back().closed = true;
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 0.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(1.0, 1.0, 0.0));
+        bottom_cc.contours.back().points.emplace_back(vec3<double>(0.0, 1.0, 0.0));
+
+        auto mesh = Triangulate_Planar_Contour_Connectivity(top_cc, bottom_cc);
+
+        for (size_t i = 0; i < mesh.faces.size(); ++i) {
+            double area = mesh.surface_area(static_cast<int64_t>(i));
+            REQUIRE(area > 0.0);
+        }
+    }
+}
+