Bistellar_flip_test.cpp

Go to the documentation of this file.

/*******************************************************************************
 Causal Dynamical Triangulations in C++ using CGAL
 
 Copyright © 2025 Adam Getchell
 ******************************************************************************/
 
 
#include "Ergodic_moves_3.hpp"
 
#include <doctest/doctest.h>
#include <numbers>
 
using namespace std;
using namespace manifolds;
using namespace ergodic_moves;
 
// Constants for testing
static inline std::floating_point auto constexpr SQRT_2 =
    std::numbers::sqrt2_v<double>;
static inline std::floating_point auto constexpr INV_SQRT_2 = 1 / SQRT_2;
 
// Helper function to create a triangulation for testing
auto create_test_triangulation() -> Delaunay
{
    // Create a 6-vertex triangulation that can support bistellar flips
    // This uses the same configuration as the reference implementation
    
    vector vertices{
        Point_t<3>{          0,           0,          0},    // Bottom vertex (t=0)
        Point_t<3>{ INV_SQRT_2,           0, INV_SQRT_2},   // Middle layer (t=1)
        Point_t<3>{          0,  INV_SQRT_2, INV_SQRT_2},   // Middle layer (t=1)
        Point_t<3>{-INV_SQRT_2,           0, INV_SQRT_2},   // Middle layer (t=1)
        Point_t<3>{          0, -INV_SQRT_2, INV_SQRT_2},   // Middle layer (t=1)
        Point_t<3>{          0,           0,          2}     // Top vertex (t=2)
    };
    
    // Create the Delaunay triangulation
    Delaunay triangulation(vertices.begin(), vertices.end());
    
    // Set time values for vertices based on z-coordinate
    for (auto vit = triangulation.finite_vertices_begin(); 
         vit != triangulation.finite_vertices_end(); ++vit)
    {
        auto z = vit->point().z();
        if (z < 0.1) {
            vit->info() = 0;  // Bottom vertex
        }
        else if (z < 1.5) {
            vit->info() = 1;  // Middle layer vertices
        }
        else {
            vit->info() = 2;  // Top vertex
        }
    }
    
    return triangulation;
}
 
// Helper function to verify that a triangulation is valid
auto verify_triangulation_validity(Delaunay const& triangulation) -> bool
{
    return triangulation.is_valid() && 
           triangulation.tds().is_valid();
}
 
// Helper function to verify that all neighbor relationships are bidirectional
auto verify_neighbor_relationships(Delaunay const& triangulation) -> bool
{
    for (auto cit = triangulation.finite_cells_begin(); 
         cit != triangulation.finite_cells_end(); ++cit)
    {
        for (int i = 0; i < 4; ++i)
        {
            auto neighbor = cit->neighbor(i);
            if (triangulation.is_infinite(neighbor))
                continue;
            
            // Check if the neighbor has this cell as its neighbor
            bool has_reciprocal_neighbor = false;
            for (int j = 0; j < 4; ++j)
            {
                Cell_handle cell_handle = cit;
                if (neighbor->neighbor(j) == cell_handle)
                {
                    has_reciprocal_neighbor = true;
                    break;
                }
            }
            
            if (!has_reciprocal_neighbor)
            {
                return false;
            }
        }
    }
    return true;
}
 
SCENARIO("Perform basic bistellar flip on Delaunay triangulation" *
         doctest::test_suite("bistellar_flip"))
{
    GIVEN("A triangulation setup for a bistellar flip")
    {
        auto triangulation = create_test_triangulation();
        
        WHEN("We have a valid triangulation")
        {
            REQUIRE(triangulation.is_valid());
            THEN("We can find a pivot edge for the bistellar flip")
            {
                auto top = triangulation.insert(Point_t<3>{0, 0, 2});
                auto bottom = triangulation.insert(Point_t<3>{0, 0, 0});
                auto edges = foliated_triangulations::collect_edges<3>(triangulation);
                auto pivot_edge = find_pivot_edge(triangulation, edges);
                
                REQUIRE_MESSAGE(pivot_edge, "No pivot edge found.");
                CHECK(triangulation.tds().is_edge(
                    pivot_edge->first, pivot_edge->second, pivot_edge->third));
                
                AND_THEN("We can perform a bistellar flip")
                {
                    auto flipped_triangulation = bistellar_flip(
                        triangulation, pivot_edge.value(), top, bottom);
                        
                    REQUIRE_MESSAGE(flipped_triangulation, "Bistellar flip failed.");
                    
                    // Basic validity check
                    CHECK(flipped_triangulation->is_valid());
                    CHECK(flipped_triangulation->tds().is_valid());
                    
                    // Verify same number of vertices, edges, and cells
                    CHECK_EQ(flipped_triangulation->number_of_vertices(), 
                             triangulation.number_of_vertices());
                    CHECK_EQ(flipped_triangulation->number_of_finite_edges(),
                             triangulation.number_of_finite_edges());
                    CHECK_EQ(flipped_triangulation->number_of_finite_cells(),
                             triangulation.number_of_finite_cells());
                }
            }
        }
    }
}
 
SCENARIO("Verify neighbor relationships after bistellar flip" *
         doctest::test_suite("bistellar_flip"))
{
    GIVEN("A triangulation setup for a bistellar flip")
    {
        auto triangulation = create_test_triangulation();
        
        WHEN("We perform a bistellar flip")
        {
            auto top = triangulation.insert(Point_t<3>{0, 0, 2});
            auto bottom = triangulation.insert(Point_t<3>{0, 0, 0});
            auto edges = foliated_triangulations::collect_edges<3>(triangulation);
            auto pivot_edge = find_pivot_edge(triangulation, edges);
            
            REQUIRE_MESSAGE(pivot_edge, "No pivot edge found.");
            
            auto flipped_triangulation = bistellar_flip(
                triangulation, pivot_edge.value(), top, bottom);
                
            REQUIRE_MESSAGE(flipped_triangulation, "Bistellar flip failed.");
            
            THEN("All neighbor relationships are bidirectional")
            {
                CHECK(verify_neighbor_relationships(flipped_triangulation.value()));
                
                AND_THEN("Each cell has exactly four neighbors")
                {
                    bool all_cells_have_four_neighbors = true;
                    for (auto cit = flipped_triangulation->finite_cells_begin(); 
                         cit != flipped_triangulation->finite_cells_end(); ++cit)
                    {
                        int neighbor_count = 0;
                        for (int i = 0; i < 4; ++i)
                        {
                            if (!flipped_triangulation->is_infinite(cit->neighbor(i)))
                                neighbor_count++;
                        }
                        
                        if (neighbor_count != 4 && neighbor_count != 3)  // Allow boundary cells with 3 neighbors
                        {
                            all_cells_have_four_neighbors = false;
                            break;
                        }
                    }
                    
                    CHECK(all_cells_have_four_neighbors);
                }
            }
        }
    }
}
 
SCENARIO("Verify cell orientation and vertex ordering after bistellar flip" *
         doctest::test_suite("bistellar_flip"))
{
    GIVEN("A triangulation setup for a bistellar flip")
    {
        auto triangulation = create_test_triangulation();
        
        WHEN("We perform a bistellar flip")
        {
            auto top = triangulation.insert(Point_t<3>{0, 0, 2});
            auto bottom = triangulation.insert(Point_t<3>{0, 0, 0});
            auto edges = foliated_triangulations::collect_edges<3>(triangulation);
            auto pivot_edge = find_pivot_edge(triangulation, edges);
            
            REQUIRE_MESSAGE(pivot_edge, "No pivot edge found.");
            
            auto flipped_triangulation = bistellar_flip(
                triangulation, pivot_edge.value(), top, bottom);
                
            REQUIRE_MESSAGE(flipped_triangulation, "Bistellar flip failed.");
            
            THEN("All cells have correct orientation")
            {
                // CGAL's is_valid() checks for proper orientation
                // FoliatedTriangulation doesn't have is_valid() method
                // Validity will be checked through the Manifold_3 below
                
                AND_THEN("Vertex ordering is consistent")
                {
                    bool consistent_ordering = true;
                    for (auto cit = flipped_triangulation->finite_cells_begin(); 
                         cit != flipped_triangulation->finite_cells_end(); ++cit)
                    {
                        // CGAL ensures that vertices are ordered so that the orientation is positive
                        // This is implicitly checked by is_valid(), but we can verify the determinant is positive
                        auto v0 = cit->vertex(0)->point();
                        auto v1 = cit->vertex(1)->point();
                        auto v2 = cit->vertex(2)->point();
                        auto v3 = cit->vertex(3)->point();
                        
                        // Orientation is checked by CGAL's is_valid, so we don't need to manually calculate it
                        // Just check that each vertex is unique
                        if (v0 == v1 || v0 == v2 || v0 == v3 || 
                            v1 == v2 || v1 == v3 || v2 == v3)
                        {
                            consistent_ordering = false;
                            break;
                        }
                    }
                    
                    CHECK(consistent_ordering);
                }
            }
        }
    }
}
 
SCENARIO("Test edge cases and error conditions for bistellar flip" *
         doctest::test_suite("bistellar_flip"))
{
    GIVEN("A triangulation setup for a bistellar flip")
    {
        auto triangulation = create_test_triangulation();
        auto top = triangulation.insert(Point_t<3>{0, 0, 2});
        auto bottom = triangulation.insert(Point_t<3>{0, 0, 0});
        auto edges = foliated_triangulations::collect_edges<3>(triangulation);
        auto pivot_edge = find_pivot_edge(triangulation, edges);
        
        REQUIRE_MESSAGE(pivot_edge, "No pivot edge found.");
        
        WHEN("We provide an invalid edge")
        {
            // Create an invalid edge
            Edge_handle invalid_edge;
            invalid_edge.first = nullptr;
            invalid_edge.second = 0;
            invalid_edge.third = 1;
            
            auto result = bistellar_flip(triangulation, invalid_edge, top, bottom);
            
            THEN("The bistellar flip should fail")
            {
                CHECK_FALSE(result.has_value());
            }
        }
        
        WHEN("We provide null vertex handles")
        {
            auto result = bistellar_flip(triangulation, pivot_edge.value(), nullptr, bottom);
            
            THEN("The bistellar flip should fail")
            {
                CHECK_FALSE(result.has_value());
            }
            
            auto result2 = bistellar_flip(triangulation, pivot_edge.value(), top, nullptr);
            
            THEN("The bistellar flip should fail")
            {
                CHECK_FALSE(result2.has_value());
            }
        }
        
        WHEN("We provide infinite vertices")
        {
            auto infinite_vertex = triangulation.infinite_vertex();
            auto result = bistellar_flip(triangulation, pivot_edge.value(), infinite_vertex, bottom);
            
            THEN("The bistellar flip should fail")
            {
                CHECK_FALSE(result.has_value());
            }
            
            auto result2 = bistellar_flip(triangulation, pivot_edge.value(), top, infinite_vertex);
            
            THEN("The bistellar flip should fail")
            {
                CHECK_FALSE(result2.has_value());
            }
        }
    }
}
 
SCENARIO("Verify that a flipped triangulation can be used in a Manifold" *
         doctest::test_suite("bistellar_flip"))
{
    GIVEN("A valid triangulation that has been flipped")
    {
        auto triangulation = create_test_triangulation();
        auto top = triangulation.insert(Point_t<3>{0, 0, 2});
        auto bottom = triangulation.insert(Point_t<3>{0, 0, 0});
        auto edges = foliated_triangulations::collect_edges<3>(triangulation);
        auto pivot_edge = find_pivot_edge(triangulation, edges);
        
        REQUIRE_MESSAGE(pivot_edge, "No pivot edge found.");
        
        auto flipped_triangulation = bistellar_flip(
            triangulation, pivot_edge.value(), top, bottom);
            
        REQUIRE_MESSAGE(flipped_triangulation, "Bistellar flip failed.");
        
        WHEN("We create a foliated triangulation from the flipped triangulation")
        {
            // Create a foliated triangulation from the flipped Delaunay triangulation
            auto foliated_triangulation = foliated_triangulations::FoliatedTriangulation<3>(
                flipped_triangulation.value(), 0, 1);
                
            THEN("The foliated triangulation is valid")
            {
                // The FoliatedTriangulation class doesn't have is_valid() method
                // Use is_correct() method which verifies the class invariants
                CHECK(foliated_triangulation.is_correct());
                
                AND_THEN("We can create a manifold from the foliated triangulation")
                {
                    Manifold_3 manifold(foliated_triangulation);
                    
                    CHECK(manifold.is_correct());
                    CHECK(manifold.is_valid());
                }
            }
        }
    }
}