diff options
Diffstat (limited to 'src/Coxeter_triangulation/example')
5 files changed, 259 insertions, 0 deletions
diff --git a/src/Coxeter_triangulation/example/CMakeLists.txt b/src/Coxeter_triangulation/example/CMakeLists.txt new file mode 100644 index 00000000..7f81c599 --- /dev/null +++ b/src/Coxeter_triangulation/example/CMakeLists.txt @@ -0,0 +1,19 @@ +project(Coxeter_triangulation_example) + +if (NOT EIGEN3_VERSION VERSION_LESS 3.1.0) + # because of random_orthogonal_matrix inclusion + if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable ( Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example manifold_tracing_flat_torus_with_boundary.cpp ) + target_link_libraries(Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example ${CGAL_LIBRARY}) + add_test(NAME Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example + COMMAND $<TARGET_FILE:Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example>) + endif() + + add_executable ( Coxeter_triangulation_manifold_tracing_custom_function_example manifold_tracing_custom_function.cpp ) + add_test(NAME Coxeter_triangulation_manifold_tracing_custom_function_example + COMMAND $<TARGET_FILE:Coxeter_triangulation_manifold_tracing_custom_function_example>) + + add_executable ( Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example cell_complex_from_basic_circle_manifold.cpp ) + add_test(NAME Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example + COMMAND $<TARGET_FILE:Coxeter_triangulation_cell_complex_from_basic_circle_manifold_example>) +endif()
\ No newline at end of file diff --git a/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp new file mode 100644 index 00000000..dfaaffa8 --- /dev/null +++ b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold.cpp @@ -0,0 +1,55 @@ +#include <iostream> + +#include <gudhi/Coxeter_triangulation.h> +#include <gudhi/Implicit_manifold_intersection_oracle.h> // for Gudhi::coxeter_triangulation::make_oracle +#include <gudhi/Manifold_tracing.h> +#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h> +#include <gudhi/Functions/Function_Sm_in_Rd.h> + +using namespace Gudhi::coxeter_triangulation; + +int main(int argc, char** argv) { + // Oracle is a circle of radius 1 + double radius = 1.; + auto oracle = make_oracle(Function_Sm_in_Rd(radius, 1)); + + // Define a Coxeter triangulation. + Coxeter_triangulation<> cox_tr(oracle.amb_d()); + // Theory forbids that a vertex of the triangulation lies exactly on the circle. + // Add some offset to avoid algorithm degeneracies. + cox_tr.change_offset(-Eigen::VectorXd::Random(oracle.amb_d())); + // For a better manifold approximation, one can change the circle radius value or change the linear transformation + // matrix. + // The number of points and edges will increase with a better resolution. + //cox_tr.change_matrix(0.5 * cox_tr.matrix()); + + // Manifold tracing algorithm + using Out_simplex_map = typename Manifold_tracing<Coxeter_triangulation<> >::Out_simplex_map; + + std::vector<Eigen::VectorXd> seed_points(1, oracle.seed()); + Out_simplex_map interior_simplex_map; + manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map); + + // Constructing the cell complex + std::size_t intr_d = oracle.amb_d() - oracle.cod_d(); + Cell_complex<Out_simplex_map> cell_complex(intr_d); + cell_complex.construct_complex(interior_simplex_map); + + // List of Hasse_cell pointers to retrieve vertices values from edges + std::map<Cell_complex<Out_simplex_map>::Hasse_cell*, std::size_t> vi_map; + std::size_t index = 0; + + std::clog << "Vertices:" << std::endl; + for (const auto& cp_pair : cell_complex.cell_point_map()) { + std::clog << index << " : (" << cp_pair.second(0) << ", " << cp_pair.second(1) << ")" << std::endl; + vi_map.emplace(cp_pair.first, index++); + } + + std::clog << "Edges:" << std::endl; + for (const auto& sc_pair : cell_complex.interior_simplex_cell_map(1)) { + Cell_complex<Out_simplex_map>::Hasse_cell* edge_cell = sc_pair.second; + for (const auto& vi_pair : edge_cell->get_boundary()) std::clog << vi_map[vi_pair.first] << " "; + std::clog << std::endl; + } + return 0; +} diff --git a/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt new file mode 100644 index 00000000..b323cca3 --- /dev/null +++ b/src/Coxeter_triangulation/example/cell_complex_from_basic_circle_manifold_for_doc.txt @@ -0,0 +1,26 @@ +Vertices: +0 : (-0.680375, 0.523483) +1 : (0.147642, 0.887879) +2 : (-0.847996, 0.30801) +3 : (-0.881369, 0.0951903) +4 : (0.638494, -0.550215) +5 : (0.415344, 0.843848) +6 : (0.812453, -0.0815816) +7 : (0.319625, -0.7709) +8 : (0.319625, 0.889605) +9 : (0.579487, 0.638553) +10 : (-0.680375, -0.461325) +11 : (-0.364269, -0.760962) +Edges: +3 2 +3 10 +10 11 +11 7 +7 4 +2 0 +0 1 +6 9 +6 4 +1 8 +8 5 +5 9 diff --git a/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp b/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp new file mode 100644 index 00000000..fe2051bb --- /dev/null +++ b/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp @@ -0,0 +1,87 @@ +#include <iostream> + +#include <gudhi/Coxeter_triangulation.h> +#include <gudhi/Functions/Function_Sm_in_Rd.h> +#include <gudhi/Implicit_manifold_intersection_oracle.h> +#include <gudhi/Manifold_tracing.h> +#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h> +#include <gudhi/Functions/Linear_transformation.h> + +#include <gudhi/IO/build_mesh_from_cell_complex.h> +#include <gudhi/IO/output_meshes_to_medit.h> + +using namespace Gudhi::coxeter_triangulation; + +/* A definition of a function that defines a 2d surface embedded in R^4, but that normally + * lives on a complex projective plane. + * In terms of harmonic coordinates [x:y:z] of points on the complex projective plane, + * the equation of the manifold is x^3*y + y^3*z + z^3*x = 0. + * The embedding consists of restricting the manifold to the affine subspace z = 1. + */ +struct Function_surface_on_CP2_in_R4 { + Eigen::VectorXd operator()(const Eigen::VectorXd& p) const { + // The real and imaginary parts of the variables x and y + double xr = p(0), xi = p(1), yr = p(2), yi = p(3); + Eigen::VectorXd result(cod_d()); + + // Squares and cubes of real and imaginary parts used in the computations + double xr2 = xr * xr, xi2 = xi * xi, yr2 = yr * yr, yi2 = yi * yi, xr3 = xr2 * xr, xi3 = xi2 * xi, yr3 = yr2 * yr, + yi3 = yi2 * yi; + + // The first coordinate of the output is Re(x^3*y + y^3 + x) + result(0) = xr3 * yr - 3 * xr * xi2 * yr - 3 * xr2 * xi * yi + xi3 * yi + yr3 - 3 * yr * yi2 + xr; + // The second coordinate of the output is Im(x^3*y + y^3 + x) + result(1) = 3 * xr2 * xi * yr + xr3 * yi - 3 * xr * xi2 * yi - xi3 * yr + 3 * yr2 * yi - yi3 + xi; + return result; + } + + std::size_t amb_d() const { return 4; }; + std::size_t cod_d() const { return 2; }; + + Eigen::VectorXd seed() const { + Eigen::VectorXd result = Eigen::VectorXd::Zero(4); + return result; + } + + Function_surface_on_CP2_in_R4() {} +}; + +int main(int argc, char** argv) { + // The function for the (non-compact) manifold + Function_surface_on_CP2_in_R4 fun; + + // Seed of the function + Eigen::VectorXd seed = fun.seed(); + + // Creating the function that defines the boundary of a compact region on the manifold + double radius = 3.0; + Function_Sm_in_Rd fun_sph(radius, 3, seed); + + // Defining the intersection oracle + auto oracle = make_oracle(fun, fun_sph); + + // Define a Coxeter triangulation scaled by a factor lambda. + // The triangulation is translated by a random vector to avoid violating the genericity hypothesis. + double lambda = 0.2; + Coxeter_triangulation<> cox_tr(oracle.amb_d()); + cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d())); + cox_tr.change_matrix(lambda * cox_tr.matrix()); + + // Manifold tracing algorithm + using MT = Manifold_tracing<Coxeter_triangulation<> >; + using Out_simplex_map = typename MT::Out_simplex_map; + std::vector<Eigen::VectorXd> seed_points(1, seed); + Out_simplex_map interior_simplex_map, boundary_simplex_map; + manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map, boundary_simplex_map); + + // Constructing the cell complex + std::size_t intr_d = oracle.amb_d() - oracle.cod_d(); + Cell_complex<Out_simplex_map> cell_complex(intr_d); + cell_complex.construct_complex(interior_simplex_map, boundary_simplex_map); + + // Output the cell complex to a file readable by medit + output_meshes_to_medit(3, "manifold_on_CP2_with_boundary", + build_mesh_from_cell_complex(cell_complex, Configuration(true, true, true, 1, 5, 3), + Configuration(true, true, true, 2, 13, 14))); + return 0; +} diff --git a/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp b/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp new file mode 100644 index 00000000..59fe2e2b --- /dev/null +++ b/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp @@ -0,0 +1,72 @@ +// workaround for the annoying boost message in boost 1.69 +#define BOOST_PENDING_INTEGER_LOG2_HPP +#include <boost/integer/integer_log2.hpp> +// end workaround + +#include <iostream> + +#include <gudhi/Coxeter_triangulation.h> +#include <gudhi/Functions/Function_affine_plane_in_Rd.h> +#include <gudhi/Functions/Function_Sm_in_Rd.h> +#include <gudhi/Functions/Cartesian_product.h> +#include <gudhi/Functions/Linear_transformation.h> +#include <gudhi/Implicit_manifold_intersection_oracle.h> +#include <gudhi/Manifold_tracing.h> +#include <gudhi/Coxeter_triangulation/Cell_complex/Cell_complex.h> +#include <gudhi/Functions/random_orthogonal_matrix.h> // requires CGAL + +#include <gudhi/IO/build_mesh_from_cell_complex.h> +#include <gudhi/IO/output_meshes_to_medit.h> + +using namespace Gudhi::coxeter_triangulation; + +int main(int argc, char** argv) { + // Creating a circle S1 in R2 of specified radius + double radius = 1.0; + Function_Sm_in_Rd fun_circle(radius, 1); + + // Creating a flat torus S1xS1 in R4 from two circle functions + auto fun_flat_torus = make_product_function(fun_circle, fun_circle); + + // Apply a random rotation in R4 + auto matrix = random_orthogonal_matrix(4); + auto fun_flat_torus_rotated = make_linear_transformation(fun_flat_torus, matrix); + + // Computing the seed of the function fun_flat_torus + Eigen::VectorXd seed = fun_flat_torus_rotated.seed(); + + // Defining a domain function that defines the boundary, which is a hyperplane passing by the origin and orthogonal to + // x. + Eigen::MatrixXd normal_matrix = Eigen::MatrixXd::Zero(4, 1); + for (std::size_t i = 0; i < 4; i++) normal_matrix(i, 0) = -seed(i); + Function_affine_plane_in_Rd fun_bound(normal_matrix, -seed / 2); + + // Defining the intersection oracle + auto oracle = make_oracle(fun_flat_torus_rotated, fun_bound); + + // Define a Coxeter triangulation scaled by a factor lambda. + // The triangulation is translated by a random vector to avoid violating the genericity hypothesis. + double lambda = 0.2; + Coxeter_triangulation<> cox_tr(oracle.amb_d()); + cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d())); + cox_tr.change_matrix(lambda * cox_tr.matrix()); + + // Manifold tracing algorithm + using MT = Manifold_tracing<Coxeter_triangulation<> >; + using Out_simplex_map = typename MT::Out_simplex_map; + std::vector<Eigen::VectorXd> seed_points(1, seed); + Out_simplex_map interior_simplex_map, boundary_simplex_map; + manifold_tracing_algorithm(seed_points, cox_tr, oracle, interior_simplex_map, boundary_simplex_map); + + // Constructing the cell complex + std::size_t intr_d = oracle.amb_d() - oracle.cod_d(); + Cell_complex<Out_simplex_map> cell_complex(intr_d); + cell_complex.construct_complex(interior_simplex_map, boundary_simplex_map); + + // Output the cell complex to a file readable by medit + output_meshes_to_medit(3, "flat_torus_with_boundary", + build_mesh_from_cell_complex(cell_complex, Configuration(true, true, true, 1, 5, 3), + Configuration(true, true, true, 2, 13, 14))); + + return 0; +} |