Added the files from the coxeter branch

3 files changed, 182 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..345949ae
--- /dev/null
+++ b/src/Coxeter_triangulation/example/CMakeLists.txt
@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 2.6)
+project(Coxeter_triangulation_example)
+
+add_executable ( Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example manifold_tracing_flat_torus_with_boundary.cpp )
+add_executable ( Coxeter_triangulation_manifold_tracing_custom_function_example manifold_tracing_custom_function.cpp )
+
+if (TBB_FOUND)
+   target_link_libraries(Coxeter_triangulation_manifold_tracing_flat_torus_with_boundary_example ${TBB_LIBRARIES})
+   target_link_libraries(Coxeter_triangulation_manifold_tracing_custom_function_example ${TBB_LIBRARIES})
+endif()
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..7a89a32f
--- /dev/null
+++ b/src/Coxeter_triangulation/example/manifold_tracing_custom_function.cpp
@@ -0,0 +1,97 @@
+#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/Cell_complex.h>
+#include <gudhi/Functions/random_orthogonal_matrix.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 : public Function {
+
+  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..c83fdd5d
--- /dev/null
+++ b/src/Coxeter_triangulation/example/manifold_tracing_flat_torus_with_boundary.cpp
@@ -0,0 +1,75 @@
+// 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/Cell_complex.h>
+#include <gudhi/Functions/random_orthogonal_matrix.h>
+
+#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;
+}