diff options
Diffstat (limited to 'src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp')
-rw-r--r-- | src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp | 72 |
1 files changed, 72 insertions, 0 deletions
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; +} |