Merge Coxeter Siargey branch

58 files changed, 6835 insertions, 0 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8b82036b..06902168 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -27,6 +27,7 @@ add_gudhi_module(Bitmap_cubical_complex)
 add_gudhi_module(Bottleneck_distance)
 add_gudhi_module(Collapse)
 add_gudhi_module(Contraction)
+add_gudhi_module(Coxeter_triangulation)
 add_gudhi_module(Cech_complex)
 add_gudhi_module(Hasse_complex)
 add_gudhi_module(Persistence_representations)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 9e4d78ac..079e3230 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -27,6 +27,7 @@ add_gudhi_module(Bottleneck_distance)
 add_gudhi_module(Cech_complex)
 add_gudhi_module(Contraction)
 add_gudhi_module(Collapse)
+add_gudhi_module(Coxeter_triangulation)
 add_gudhi_module(Hasse_complex)
 add_gudhi_module(Persistence_representations)
 add_gudhi_module(Persistent_cohomology)
diff --git a/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h b/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h
new file mode 100644
index 00000000..2339817e
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/FunctionForImplicitManifold.h
@@ -0,0 +1,46 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_FUNCTION_FOR_IMPLICIT_MANIFOLD_H_
+#define CONCEPT_COXETER_TRIANGULATION_FUNCTION_FOR_IMPLICIT_MANIFOLD_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept FunctionForImplicitManifold describes the requirements 
+ * for a type to implement an implicit function class used for example in Manifold_tracing.
+ */
+struct FunctionForImplicitManifold {
+
+  /** \brief Value of the function at a specified point 'p'. 
+   *  @param[in] p The input point given by its Cartesian coordinates.
+   *  Its size needs to be equal to amb_d().
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const;
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const;
+  
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const;
+
+  /** \brief Returns a point on the zero-set of the function. */
+  Eigen::VectorXd seed() const;
+};
+
+}  // namespace coxeter_triangulation
+
+}  // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/IntersectionOracle.h b/src/Coxeter_triangulation/concept/IntersectionOracle.h
new file mode 100644
index 00000000..642317f2
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/IntersectionOracle.h
@@ -0,0 +1,109 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_INTERSECTION_ORACLE_H_
+#define CONCEPT_COXETER_TRIANGULATION_INTERSECTION_ORACLE_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept IntersectionOracle describes the requirements 
+ * for a type to implement an intersection oracle class used for example in Manifold_tracing.
+ *
+ */
+struct IntersectionOracle {
+
+  /** \brief Returns the domain (ambient) dimension of the underlying manifold. */
+  std::size_t amb_d() const;
+  
+  /** \brief Returns the codomain dimension of the underlying manifold. */
+  std::size_t cod_d() const;
+
+  /** \brief Intersection query with the relative interior of the manifold.
+   *
+   *  \details The returned structure Query_result contains the boolean value
+   *   that is true only if the intersection point of the query simplex and
+   *   the relative interior of the manifold exists, the intersection point
+   *   and the face of the query simplex that contains 
+   *   the intersection point.
+   *   
+   *  \tparam Simplex_handle The class of the query simplex.
+   *   Needs to be a model of the concept SimplexInCoxeterTriangulation.
+   *  \tparam Triangulation The class of the triangulation.
+   *   Needs to be a model of the concept TriangulationForManifoldTracing.
+   *  
+   *  @param[in] simplex The query simplex. The dimension of the simplex
+   *   should be the same as the codimension of the manifold 
+   *   (the codomain dimension of the function).
+   *  @param[in] triangulation The ambient triangulation. The dimension of 
+   *   the triangulation should be the same as the ambient dimension of the manifold 
+   *   (the domain dimension of the function).
+   */
+  template <class Simplex_handle,
+	    class Triangulation>
+  Query_result<Simplex_handle> intersects(const Simplex_handle& simplex,
+					  const Triangulation& triangulation) const;
+
+  /** \brief Intersection query with the boundary of the manifold.
+   *  
+   *  \details The returned structure Query_result contains the boolean value
+   *   that is true only if the intersection point of the query simplex and
+   *   the boundary of the manifold exists, the intersection point
+   *   and the face of the query simplex that contains 
+   *   the intersection point.
+   *   
+   *  \tparam Simplex_handle The class of the query simplex.
+   *   Needs to be a model of the concept SimplexInCoxeterTriangulation.
+   *  \tparam Triangulation The class of the triangulation.
+   *   Needs to be a model of the concept TriangulationForManifoldTracing.
+   *
+   *  @param[in] simplex The query simplex. The dimension of the simplex
+   *   should be the same as the codimension of the boundary of the manifold 
+   *   (the codomain dimension of the function + 1).
+   *  @param[in] triangulation The ambient triangulation. The dimension of 
+   *   the triangulation should be the same as the ambient dimension of the manifold 
+   *   (the domain dimension of the function).
+   */
+  template <class Simplex_handle,
+	    class Triangulation>
+  Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex,
+						   const Triangulation& triangulation) const;
+
+  /** \brief Returns true if the input point lies inside the piecewise-linear
+   *   domain induced by the given ambient triangulation that defines the relative
+   *   interior of the piecewise-linear approximation of the manifold.
+   *
+   * @param p The input point. Needs to have the same dimension as the ambient
+   *  dimension of the manifold (the domain dimension of the function).
+   * @param triangulation The ambient triangulation. Needs to have the same
+   *  dimension as the ambient dimension of the manifold 
+   *  (the domain dimension of the function).
+   */
+  template <class Triangulation>
+  bool lies_in_domain(const Eigen::VectorXd& p,
+		      const Triangulation& triangulation) const {
+    Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p);
+    return pl_p(0) < 0;
+  }
+
+  /** \brief Returns the function that defines the interior of the manifold */
+  const Function_& function() const;
+
+};
+
+}  // namespace coxeter_triangulation
+
+}  // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h b/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h
new file mode 100644
index 00000000..c1562066
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/SimplexInCoxeterTriangulation.h
@@ -0,0 +1,82 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_SIMPLEX_IN_COXETER_TRIANGULATION_H_
+#define CONCEPT_COXETER_TRIANGULATION_SIMPLEX_IN_COXETER_TRIANGULATION_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept SimplexInCoxeterTriangulation describes the requirements 
+ * for a type to implement a representation of simplices in Freudenthal_triangulation
+ * or in Coxeter_triangulation.
+ */
+struct SimplexInCoxeterTriangulation {
+
+  /** \brief Type of the vertex. */
+  typedef Vertex_ Vertex;
+  
+  /** \brief Type of the ordered partition. */
+  typedef Ordered_set_partition_ OrderedSetPartition;
+
+  /** \brief Dimension of the simplex. */
+  unsigned dimension() const;
+
+  /** \brief Type of a range of vertices, each of type Vertex. */
+  typedef Vertex_range;
+  
+  /** \brief Returns a range of vertices of the simplex.
+   */
+  Vertex_range vertex_range() const;
+
+  /** \brief Type of a range of faces, each of type that 
+   *  is a model of the concept SimplexInCoxeterTriangulation. 
+   */
+  typedef Face_range;
+
+  /** \brief Returns a range of permutahedral representations of k-dimensional faces
+   *   of the simplex for some given integer parameter 'k'.
+   */
+  Face_range face_range(std::size_t k) const;
+
+  /** \brief Returns a range of permutahedral representations of facets of the simplex.
+   * The dimension of the simplex must be strictly positive.
+   */
+  Face_range facet_range() const;
+
+  /** \brief Type of a range of cofaces, each of type that 
+   *  is a model of the concept SimplexInCoxeterTriangulation. 
+   */
+  typedef Coface_range;
+
+  /** \brief Returns a range of permutahedral representations of k-dimensional cofaces
+   *   of the simplex for some given integer parameter 'k'.
+   */
+  Coface_range coface_range(std::size_t k) const;
+
+  /** \brief Returns a range of permutahedral representations of cofacets of the simplex.
+   * The dimension of the simplex must be strictly different from the ambient dimension.
+   */
+  Coface_range cofacet_range() const;
+
+  /** \brief Returns true, if the simplex is a face of other simplex. */
+  bool is_face_of(const Permutahedral_representation& other) const;
+  
+};
+
+}  // namespace coxeter_triangulation
+
+}  // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h b/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h
new file mode 100644
index 00000000..e64d3b97
--- /dev/null
+++ b/src/Coxeter_triangulation/concept/TriangulationForManifoldTracing.h
@@ -0,0 +1,59 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CONCEPT_COXETER_TRIANGULATION_TRIANGULATION_FOR_MANIFOLD_TRACING_H_
+#define CONCEPT_COXETER_TRIANGULATION_TRIANGULATION_FOR_MANIFOLD_TRACING_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief The concept TriangulationForManifoldTracing describes the requirements 
+ * for a type to implement a triangulation class used for example in Manifold_tracing.
+ */
+struct TriangulationForManifoldTracing {
+
+  /** \brief Type of the simplices in the triangulation. 
+   *  Needs to be a model of the concept SimplexInCoxeterTriangulation. */
+  typedef Simplex_handle;
+
+  /** \brief Type of the vertices in the triangulation. 
+   *  Needs to be a random-access range of integer values. */
+  typedef Vertex_handle;
+  
+  /** \brief Returns the permutahedral representation of the simplex in the
+   *  triangulation that contains a given query point 'p'.
+   * \tparam Point_d A class that represents a point in d-dimensional Euclidean space.
+   * The coordinates should be random-accessible. Needs to provide the method size().
+   * @param[in] point The query point.
+   */
+  template <class Point_d>
+  Simplex_handle locate_point(const Point_d& point) const;
+
+  /** \brief Returns the Cartesian coordinates of the given vertex 'v'. 
+   * @param[in] v The input vertex.
+   */
+  Eigen::VectorXd cartesian_coordinates(const Vertex_handle& v) const;
+
+  /** \brief Returns the Cartesian coordinates of the barycenter of a given simplex 's'. 
+   * @param[in] s The input simplex given by permutahedral representation.
+   */
+  Eigen::VectorXd barycenter(const Simplex_handle& s) const;
+  
+};
+
+}  // namespace coxeter_triangulation
+
+}  // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/doc/custom_function.png b/src/Coxeter_triangulation/doc/custom_function.png
new file mode 100644
index 00000000..8bb8ba9a
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/custom_function.png
diff --git a/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png
new file mode 100644
index 00000000..338b39fe
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/flat_torus_with_boundary.png
diff --git a/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h b/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h
new file mode 100644
index 00000000..80ca5b37
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/intro_coxeter_triangulation.h
@@ -0,0 +1,146 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef DOC_TANGENTIAL_COMPLEX_INTRO_COXETER_TRIANGULATION_H_
+#define DOC_TANGENTIAL_COMPLEX_INTRO_COXETER_TRIANGULATION_H_
+
+// needs namespaces for Doxygen to link on classes
+namespace Gudhi {
+namespace tangential_complex {
+
+/**  \defgroup coxeter_triangulation Coxeter triangulation
+
+\author    Siargey Kachanovich
+
+@{
+
+\section overview Module overview
+
+Coxeter triangulation module is designed to provide tools for constructing a piecewise-linear approximation of an \f$m\f$-dimensional smooth manifold embedded in \f$ \mathbb{R}^d \f$ using an ambient triangulation.
+For a more detailed description of the module see \cite KachanovichThesis.
+
+\section manifoldtracing Manifold tracing algorithm
+The central piece of the module is the manifold tracing algorithm represented by the class \ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing".
+The manifold tracing algorithm takes as input a manifold of some dimension \f$m\f$ embedded in \f$\mathbb{R}^d\f$ represented by an intersection oracle (see Section \ref intersectionoracle "Intersection oracle"), a point on the manifold and an ambient triangulation (see Section \ref ambienttriangulations "Ambient triangulations").
+The output consists of one map (or two maps in the case of manifolds with boundary) from the \f$(d-m)\f$-dimensional (and \f$(d-m+1)\f$-dimensional in the case of manifolds with boundary) simplices in the ambient triangulation that intersect the manifold to their intersection points.
+From this output, it is possible to construct the cell complex of the piecewise-linear approximation of the input manifold.
+
+There are two methods that execute the manifold tracing algorithm: the method \ref Gudhi::coxeter_triangulation::Manifold_tracing::manifold_tracing_algorithm() "Manifold_tracing::manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map)" for manifolds without boundary and \ref Gudhi::coxeter_triangulation::Manifold_tracing::manifold_tracing_algorithm() "Manifold_tracing::manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map)" for manifolds with boundary.
+The algorithm functions as follows.
+It starts at the specified seed points and inserts a \f$(d-m)\f$-dimensional simplices nearby each seed point that intersect the manifold into the output.
+Starting from this simplex, the algorithm propagates the search for other \f$(d-m)\f$-dimensional simplices that intersect the manifold by marching from a simplex to neighbouring simplices via their common cofaces.
+
+This class \ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing" has one template parameter Triangulation_ which specifies the ambient triangulation which is used by the algorithm.
+The template type Triangulation_ has to be a model of the concept \ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing".
+
+The module also provides two static methods: \ref Gudhi::coxeter_triangulation::manifold_tracing_algorithm() "manifold_tracing_algorithm(seed_points, triangulation, oracle, out_simplex_map)" for manifolds without boundary and \ref manifold_tracing_algorithm() "manifold_tracing_algorithm(seed_points, triangulation, oracle, interior_simplex_map, boundary_simplex_map)" for manifolds with boundary.
+For these static methods it is not necessary to specify any template arguments.
+
+\section ambienttriangulations Ambient triangulations
+
+The ambient triangulations supported by the manifold tracing algorithm have to be models of the concept \ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing". 
+This module offers two such models: the class \ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" and the derived class \ref Gudhi::coxeter_triangulation::Coxeter_triangulation "Coxeter_triangulation". 
+
+Both these classes encode affine transformations of the so-called Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$.
+The Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$ is defined as the simplicial subdivision of the unit cubic partition of \f$\mathbb{R}^d\f$.
+Each simplex is encoded using the permutahedral representation, which consists of an integer-valued vector \f$y\f$ that positions the simplex in a specific cube in the cubical partition and an ordered partition \f$\omega\f$ of the set \f$\{1,\ldots,d+1\}\f$, which positions the simplex in the simplicial subdivision of the cube.
+The default constructor \ref Gudhi::coxeter_triangulation::Freudenthal_triangulation::Freudenthal_triangulation(std::size_t) "Freudenthal_triangulation(d)" the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$.
+The class \ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" can also encode any affine transformation of the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$ using an invertible matrix \f$\Lambda\f$ and an offset vector \f$b\f$ that can be specified in the constructor and which can be changed using the methods change_matrix and change_offset.
+The class \ref Gudhi::coxeter_triangulation::Coxeter_triangulation "Coxeter_triangulation" is derived from \ref Gudhi::coxeter_triangulation::Freudenthal_triangulation "Freudenthal_triangulation" and its default constructor \ref Gudhi::coxeter_triangulation::Coxeter_triangulation::Coxeter_triangulation(std::size_t) "Coxeter_triangulation(d)" builds a Coxeter triangulation of type \f$\tilde{A}_d\f$, which has the best simplex quality of all linear transformations of the Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$.
+
+\image html two_triangulations.png "On the left: Coxeter triangulation of type \f$\tilde{A}_2\f$. On the right: Freudenthal-Kuhn triangulation of \f$\mathbb{R}^d\f$."
+
+
+\section intersectionoracle Intersection oracle
+
+The input \f$m\f$-dimensional manifold in \f$\mathbb{R}^d\f$ needs to be given via the intersection oracle that answers the following query: given a \f$(d-m)\f$-dimensional simplex, does it intersect the manifold?
+The concept \ref Gudhi::coxeter_triangulation::IntersectionOracle "IntersectionOracle" describes all requirements for an intersection oracle class to be compatible with the class \ref Gudhi::coxeter_triangulation::Manifold_tracing "Manifold_tracing".
+This module offers one model of the concept \ref Gudhi::coxeter_triangulation::IntersectionOracle "IntersectionOracle", which is the class \ref Gudhi::coxeter_triangulation::Implicit_manifold_intersection_oracle "Implicit_manifold_intersection_oracle".
+This class represents a manifold given as the zero-set of a specified function \f$F: \mathbb{R}^d \rightarrow \mathbb{R}^{d-m}\f$.
+The function \f$F\f$ is given by a class which is a model of the concept \ref Gudhi::coxeter_triangulation::FunctionForImplicitManifold "FunctionForImplicitManifold".
+There are multiple function classes that are already implemented in this module.
+
+\li \ref Gudhi::coxeter_triangulation::Constant_function(std::size_t, std::size_t, Eigen::VectorXd) "Constant_function(d,k,v)" defines a constant function \f$F\f$ such that for all \f$x \in \mathbb{R}^d\f$, we have \f$F(x) = v \in \mathbb{R}^k\f$.
+  The class Constant_function does not define an implicit manifold, but is useful as the domain function when defining boundaryless implicit manifolds.
+\li \ref Gudhi::coxeter_triangulation::Function_affine_plane_in_Rd(N,b) "Function_affine_plane_in_Rd(N,b)" defines an \f$m\f$-dimensional implicit affine plane in the \f$d\f$-dimensional Euclidean space given by a normal matrix \f$N\f$ and an offset vector \f$b\f$.
+\li \ref Gudhi::coxeter_triangulation::Function_Sm_in_Rd(r,m,d,center) "Function_Sm_in_Rd(r,m,d,center)" defines an \f$m\f$-dimensional implicit sphere embedded in the \f$d\f$-dimensional Euclidean space of radius \f$r\f$ centered at the point 'center'.
+\li \ref Gudhi::coxeter_triangulation::Function_moment_curve_in_Rd(r,d) "Function_moment_curve(r,d)" defines the moment curve in the \f$d\f$-dimensional Euclidean space of radius \f$r\f$ given as the parameterized curve (but implemented as an implicit curve):
+  \f[ (r, rt, \ldots, rt^{d-1}) \in \mathbb{R}^d,\text{ for $t \in \mathbb{R}$.} \f]
+\li \ref Gudhi::coxeter_triangulation::Function_torus_in_R3(R, r) "Function_torus_in_R3(R, r)" defines a torus in \f$\mathbb{R}^3\f$ with the outer radius \f$R\f$ and the inner radius, given by the equation:
+  \f[ z^2 + (\sqrt{x^2 + y^2} - r)^2 - R^2 = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_chair_in_R3(a, b, k) "Function_chair_in_R3(a, b, k)" defines the ``Chair'' surface in \f$\mathbb{R}^3\f$ defined by the equation:
+  \f[ (x^2 + y^2 + z^2 - ak^2)^2 - b((z-k)^2 - 2x^2)((z+k)^2 - 2y^2) = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_iron_in_R3() "Function_iron_in_R3()" defines the ``Iron'' surface in \f$\mathbb{R}^3\f$ defined by the equation:
+  \f[ \frac{-x^6-y^6-z^6}{300} + \frac{xy^2z}{2.1} + y^2 + (z-2)^2 = 1. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_lemniscate_revolution_in_R3(a) "Function_lemniscate_revolution_in_R3(a)" defines a revolution surface in \f$\mathbb{R}^3\f$ obtained from the lemniscate of Bernoulli defined by the equation:
+  \f[ (x^2 + y^2 + z^2)^2 - 2a^2(x^2 - y^2 - z^2) = 0. \f]
+\li \ref Gudhi::coxeter_triangulation::Function_whitney_umbrella_in_R3() "Function_whitney_umbrella_in_R3()" defines the Whitney umbrella surface in \f$\mathbb{R}^3\f$ defined by the equation:
+  \f[ x^2 - y^2z = 0. \f]
+
+The base function classes above can be composed or modified into new functions using the following classes and methods:
+
+\li \ref Gudhi::coxeter_triangulation::Cartesian_product "Cartesian_product(functions...)" expresses the Cartesian product \f$F_1^{-1}(0) \times \ldots \times F_k^{-1}(0)\f$ of multiple implicit manifolds as an implicit manifold.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::make_product_function() "make_product_function(functions...)" is provided that takes a pack of function-typed objects as the argument.
+\li \ref Gudhi::coxeter_triangulation::Embed_in_Rd "Embed_in_Rd(F, d)" expresses an implicit manifold given as the zero-set of a function \f$F\f$ embedded in a higher-dimensional Euclidean space \f$\mathbb{R}^d\f$.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::make_embedding() "make_embedding(F, d)" is provided.
+\li \ref Gudhi::coxeter_triangulation::Linear_transformation "Linear_transformation(F, M)" applies a linear transformation given by a matrix \f$M\f$ on an implicit manifold given as the zero-set of the function \f$F\f$.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::linear_transformation() "linear_transformation(F, M)" is provided.
+\li \ref Gudhi::coxeter_triangulation::Translate "Translate(F, v)" translates an implicit manifold given as the zero-set of ththe function \f$F\f$ by a vector \f$v\f$.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::translate(F, v) "translate(F, v)" is provided.
+\li \ref Gudhi::coxeter_triangulation::Negation() "Negation(F)" defines the negative of the given function \f$F\f$.
+  This class is useful to define the complementary of a given domain, when defining a manifold with boundary.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::negation() "negation(F)" is provided.
+\li \ref Gudhi::coxeter_triangulation::PL_approximation "PL_approximation(F, T)" defines a piecewise-linear approximation of a given function \f$F\f$ induced by an ambient triangulation \f$T\f$.
+  The purpose of this class is to define a piecewise-linear function that is compatible with the requirements for the domain function \f$D\f$ when defining a manifold with boundary.
+  For convenience, a static function \ref Gudhi::coxeter_triangulation::make_pl_approximation() "make_pl_approximation(F, T)" is provided.
+  The type of \f$T\f$ is required to be a model of the concept \ref Gudhi::coxeter_triangulation::TriangulationForManifoldTracing "TriangulationForManifoldTracing".
+
+\section cellcomplex Cell complex construction
+
+The output of the manifold tracing algorithm can be transformed into the Hasse diagram of a cell complex that approximates the input manifold using the class \ref Gudhi::coxeter_triangulation::Cell_complex "Cell_complex".
+The type of the cells in the Hasse diagram is \ref Gudhi::Hasse_diagram::Hasse_diagram_cell "Hasse_cell<int, double, bool>" provided by the module Hasse diagram.
+The cells in the cell complex given by an object of the class \ref Gudhi::coxeter_triangulation::Cell_complex "Cell_complex" are accessed through several maps that are accessed through the following methods.
+
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::interior_simplex_cell_maps() "interior_simplex_cell_maps()" returns a vector of maps from the cells of various dimensions in the interior of the cell complex to the permutahedral representations of the corresponding simplices in the ambient triangulation.
+Each individual map for cells of a specific dimension \f$l\f$ can be accessed using the method \ref Gudhi::coxeter_triangulation::Cell_complex::interior_simplex_cell_map() "interior_simplex_cell_map(l)".
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::boundary_simplex_cell_maps() "boundary_simplex_cell_maps()" returns a vector of maps from the cells of various dimensions on the boundary of the cell complex to the permutahedral representations of the corresponding simplices in the ambient triangulation.
+Each individual map for cells of a specific dimension \f$l\f$ can be accessed using the method \ref Gudhi::coxeter_triangulation::Cell_complex::boundary_simplex_cell_map() "boundary_simplex_cell_map(l)".
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::cell_simplex_map() "cell_simplex_map()" returns a map from the cells in the cell complex to the permutahedral representations of the corresponding simplices in the ambient triangulation.
+\li The method \ref Gudhi::coxeter_triangulation::Cell_complex::cell_point_map() "cell_point_map()" returns a map from the vertex cells in the cell complex to their Cartesian coordinates.
+
+\section example Examples
+
+Here is an example of constructing a piecewise-linear approximation of a flat torus embedded in \f$\mathbb{R}^4\f$, rotated by a random rotation in \f$\mathbb{R}^4\f$ and cut by a hyperplane.
+
+\include Coxeter_triangulation/manifold_tracing_flat_torus_with_boundary.cpp
+
+The output in medit is:
+
+\image html "flat_torus_with_boundary.png" "Output from the example of a flat torus with boundary"
+
+In the following more complex example, we define a custom function for the implicit manifold.
+
+\include Coxeter_triangulation/manifold_tracing_custom_function.cpp
+
+The output in medit looks as follows:
+
+\image html "custom_function.png" "Output from the example with a custom function"
+
+
+ */
+/** @} */  // end defgroup coxeter_triangulation
+
+}  // namespace coxeter_triangulation
+
+}  // namespace Gudhi
+
+#endif  // DOC_TANGENTIAL_COMPLEX_INTRO_TANGENTIAL_COMPLEX_H_
diff --git a/src/Coxeter_triangulation/doc/two_triangulations.png b/src/Coxeter_triangulation/doc/two_triangulations.png
new file mode 100644
index 00000000..055d93e7
--- /dev/null
+++ b/src/Coxeter_triangulation/doc/two_triangulations.png
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;
+}
diff --git a/src/Coxeter_triangulation/include/gudhi/Cell_complex.h b/src/Coxeter_triangulation/include/gudhi/Cell_complex.h
new file mode 100644
index 00000000..51109bcb
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Cell_complex.h
@@ -0,0 +1,376 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef CELL_COMPLEX_H_
+#define CELL_COMPLEX_H_
+
+#include <Eigen/Dense> 
+
+#include <algorithm>
+#include <gudhi/Permutahedral_representation/Simplex_comparator.h>
+#include <fstream> // for Hasse_diagram_persistence.h
+
+#include <gudhi/Cell_complex/Hasse_diagram_cell.h> // for Hasse_cell
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ *  \ingroup coxeter_triangulation
+ */
+
+/** \class Cell_complex 
+ *  \brief A class that constructs the cell complex from the output provided by the class 
+ *   Gudhi::Manifold_tracing<Triangulation_>.
+ *
+ *  \tparam Out_simplex_map_ The type of a map from a simplex type that is a
+ *   model of SimplexInCoxeterTriangulation to Eigen::VectorXd.
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template <class Out_simplex_map_>
+class Cell_complex {
+public:
+  
+  /** \brief Type of a simplex in the ambient triangulation.
+   *  Is a model of the concept SimplexInCoxeterTriangulation.
+   */
+  typedef typename Out_simplex_map_::key_type Simplex_handle;
+  /** \brief Type of a cell in the cell complex. 
+   *  Always is Gudhi::Hasse_cell from the Hasse diagram module.
+   *  The additional information is the boolean that is true if and only if the cell lies
+   *  on the boundary.
+   */
+  typedef Gudhi::Hasse_diagram::Hasse_diagram_cell<int, double, bool> Hasse_cell;
+  /** \brief Type of a map from permutahedral representations of simplices in the
+   *  ambient triangulation to the corresponding cells in the cell complex of some
+   *  specific dimension.
+   */
+  typedef std::map<Simplex_handle,
+		   Hasse_cell*,
+		   Simplex_comparator<Simplex_handle> > Simplex_cell_map;
+  /** \brief Type of a vector of maps from permutahedral representations of simplices in the
+   *  ambient triangulation to the corresponding cells in the cell complex of various dimensions.
+   */
+  typedef std::vector<Simplex_cell_map> Simplex_cell_maps;
+  
+  /** \brief Type of a map from cells in the cell complex to the permutahedral representations
+   *  of the corresponding simplices in the ambient triangulation.
+   */
+  typedef std::map<Hasse_cell*, Simplex_handle> Cell_simplex_map;
+
+  /** \brief Type of a map from vertex cells in the cell complex to the permutahedral representations
+   *  of their Cartesian coordinates.
+   */
+  typedef std::map<Hasse_cell*, Eigen::VectorXd> Cell_point_map;
+
+private:
+    
+  Hasse_cell* insert_cell(const Simplex_handle& simplex,
+			  std::size_t cell_d,
+			  bool is_boundary) {
+    Simplex_cell_maps& simplex_cell_maps
+      = (is_boundary? boundary_simplex_cell_maps_ : interior_simplex_cell_maps_); 
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    CC_detail_list& cc_detail_list =
+      (is_boundary? cc_boundary_detail_lists[cell_d] : cc_interior_detail_lists[cell_d]);
+    cc_detail_list.emplace_back(CC_detail_info(simplex));
+#endif
+    Simplex_cell_map& simplex_cell_map = simplex_cell_maps[cell_d];
+    auto map_it = simplex_cell_map.find(simplex);
+    if (map_it == simplex_cell_map.end()) {
+      hasse_cells_.push_back(new Hasse_cell(is_boundary, cell_d));
+      Hasse_cell* new_cell = hasse_cells_.back();    
+      simplex_cell_map.emplace(std::make_pair(simplex, new_cell));
+      cell_simplex_map_.emplace(std::make_pair(new_cell, simplex));
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      cc_detail_list.back().status_ = CC_detail_info::Result_type::inserted;
+#endif
+      return new_cell;
+    }
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    CC_detail_info& cc_info = cc_detail_list.back();
+    cc_info.trigger_ = to_string(map_it->first);
+    cc_info.status_ = CC_detail_info::Result_type::self;
+#endif
+    return map_it->second;
+  }
+
+  void expand_level(std::size_t cell_d) {
+    bool is_manifold_with_boundary = boundary_simplex_cell_maps_.size() > 0;
+    for (auto& sc_pair: interior_simplex_cell_maps_[cell_d - 1]) {
+      const Simplex_handle& simplex = sc_pair.first;
+      Hasse_cell* cell = sc_pair.second;
+      for (Simplex_handle coface: simplex.coface_range(cod_d_ + cell_d)) {
+	Hasse_cell* new_cell = insert_cell(coface, cell_d, false);
+	new_cell->get_boundary().emplace_back(std::make_pair(cell, 1));
+      }
+    }
+    
+    if (is_manifold_with_boundary) {
+      for (auto& sc_pair: boundary_simplex_cell_maps_[cell_d - 1]) {
+	const Simplex_handle& simplex = sc_pair.first;
+	Hasse_cell* cell = sc_pair.second;
+	if (cell_d != intr_d_)
+	  for (Simplex_handle coface: simplex.coface_range(cod_d_ + cell_d + 1)) {
+	    Hasse_cell* new_cell = insert_cell(coface, cell_d, true);
+	    new_cell->get_boundary().emplace_back(std::make_pair(cell, 1));
+	  }
+	auto map_it = interior_simplex_cell_maps_[cell_d].find(simplex);
+	if (map_it == interior_simplex_cell_maps_[cell_d].end())
+	  std::cerr << "Cell_complex::expand_level error: A boundary cell does not have an interior counterpart.\n";
+	else {
+	  Hasse_cell* i_cell = map_it->second;
+	  i_cell->get_boundary().emplace_back(std::make_pair(cell, 1));
+	}
+      }
+    }
+  }
+  
+  void construct_complex_(const Out_simplex_map_& out_simplex_map) {
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size());
+    cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size());
+    cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size());
+#endif    
+    for (auto& os_pair: out_simplex_map) {
+      const Simplex_handle& simplex = os_pair.first;
+      const Eigen::VectorXd& point = os_pair.second;
+      Hasse_cell* new_cell = insert_cell(simplex, 0, false);
+      cell_point_map_.emplace(std::make_pair(new_cell, point));
+    }
+    for (std::size_t cell_d = 1;
+	 cell_d < interior_simplex_cell_maps_.size() &&
+	   !interior_simplex_cell_maps_[cell_d - 1].empty();
+	 ++cell_d) {
+      expand_level(cell_d);
+    }
+  }
+  
+  void construct_complex_(const Out_simplex_map_& interior_simplex_map,
+			  const Out_simplex_map_& boundary_simplex_map) {
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    cc_interior_summary_lists.resize(interior_simplex_cell_maps_.size());
+    cc_interior_prejoin_lists.resize(interior_simplex_cell_maps_.size());
+    cc_interior_detail_lists.resize(interior_simplex_cell_maps_.size());
+    cc_boundary_summary_lists.resize(boundary_simplex_cell_maps_.size());
+    cc_boundary_prejoin_lists.resize(boundary_simplex_cell_maps_.size());
+    cc_boundary_detail_lists.resize(boundary_simplex_cell_maps_.size());
+#endif    
+    for (auto& os_pair: boundary_simplex_map) {
+      const Simplex_handle& simplex = os_pair.first;
+      const Eigen::VectorXd& point = os_pair.second;
+      Hasse_cell* new_cell = insert_cell(simplex, 0, true);
+      cell_point_map_.emplace(std::make_pair(new_cell, point));
+    }
+    for (auto& os_pair: interior_simplex_map) {
+      const Simplex_handle& simplex = os_pair.first;
+      const Eigen::VectorXd& point = os_pair.second;
+      Hasse_cell* new_cell = insert_cell(simplex, 0, false);
+      cell_point_map_.emplace(std::make_pair(new_cell, point));
+    }
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    for (const auto& sc_pair: interior_simplex_cell_maps_[0])
+      cc_interior_summary_lists[0].push_back(CC_summary_info(sc_pair));
+    for (const auto& sc_pair: boundary_simplex_cell_maps_[0])
+      cc_boundary_summary_lists[0].push_back(CC_summary_info(sc_pair));
+#endif        
+
+    for (std::size_t cell_d = 1;
+	 cell_d < interior_simplex_cell_maps_.size() &&
+	   !interior_simplex_cell_maps_[cell_d - 1].empty();
+	 ++cell_d) {
+      expand_level(cell_d);
+      
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      for (const auto& sc_pair: interior_simplex_cell_maps_[cell_d])
+	cc_interior_summary_lists[cell_d].push_back(CC_summary_info(sc_pair));
+      if (cell_d < boundary_simplex_cell_maps_.size())
+	for (const auto& sc_pair: boundary_simplex_cell_maps_[cell_d])
+	  cc_boundary_summary_lists[cell_d].push_back(CC_summary_info(sc_pair));
+#endif
+    }
+  }
+  
+public:  
+  
+  /**
+   * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold
+   *  without boundary embedded in the \f$ d \f$-dimensional Euclidean space
+   *  from the output of the class Gudhi::Manifold_tracing.
+   *
+   * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$ 
+   * in the ambient triangulation that intersect the relative interior of the manifold 
+   * to the intersection points.
+   */
+  void construct_complex(const Out_simplex_map_& out_simplex_map) {
+    interior_simplex_cell_maps_.resize(intr_d_ + 1);
+    if (!out_simplex_map.empty())
+      cod_d_ = out_simplex_map.begin()->first.dimension();
+    construct_complex_(out_simplex_map);
+  }
+  
+  /**
+   * \brief Constructs the skeleton of the cell complex that approximates
+   *  an \f$m\f$-dimensional manifold without boundary embedded 
+   *  in the \f$d\f$-dimensional Euclidean space
+   *  up to a limit dimension from the output of the class Gudhi::Manifold_tracing.
+   *
+   * \param[in] out_simplex_map A map from simplices of dimension \f$(d-m)\f$ 
+   * in the ambient triangulation that intersect the relative interior of the manifold 
+   * to the intersection points.
+   * \param[in] limit_dimension The dimension of the constructed skeleton.
+   */
+  void construct_complex(const Out_simplex_map_& out_simplex_map,
+			 std::size_t limit_dimension) {
+    interior_simplex_cell_maps_.resize(limit_dimension + 1);
+    if (!out_simplex_map.empty())
+      cod_d_ = out_simplex_map.begin()->first.dimension();
+    construct_complex_(out_simplex_map);
+  }
+
+  /**
+   * \brief Constructs the the cell complex that approximates an \f$m\f$-dimensional manifold
+   *  with boundary embedded in the \f$ d \f$-dimensional Euclidean space
+   *  from the output of the class Gudhi::Manifold_tracing.
+   *
+   * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$ 
+   * in the ambient triangulation that intersect the relative interior of the manifold 
+   * to the intersection points.
+   * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$ 
+   * in the ambient triangulation that intersect the boundary of the manifold 
+   * to the intersection points.
+   */
+  void construct_complex(const Out_simplex_map_& interior_simplex_map,
+			 const Out_simplex_map_& boundary_simplex_map) {
+    interior_simplex_cell_maps_.resize(intr_d_ + 1);
+    boundary_simplex_cell_maps_.resize(intr_d_);
+    if (!interior_simplex_map.empty())
+      cod_d_ = interior_simplex_map.begin()->first.dimension();
+    construct_complex_(interior_simplex_map, boundary_simplex_map);
+  }
+
+  /**
+   * \brief Constructs the skeleton of the cell complex that approximates
+   *  an \f$m\f$-dimensional manifold with boundary embedded 
+   *  in the \f$d\f$-dimensional Euclidean space
+   *  up to a limit dimension from the output of the class Gudhi::Manifold_tracing.
+   *
+   * \param[in] interior_simplex_map A map from simplices of dimension \f$(d-m)\f$ 
+   * in the ambient triangulation that intersect the relative interior of the manifold 
+   * to the intersection points.
+   * \param[in] boundary_simplex_map A map from simplices of dimension \f$(d-m+1)\f$ 
+   * in the ambient triangulation that intersect the boundary of the manifold 
+   * to the intersection points.
+   * \param[in] limit_dimension The dimension of the constructed skeleton.
+   */
+  void construct_complex(const Out_simplex_map_& interior_simplex_map,
+			 const Out_simplex_map_& boundary_simplex_map,
+			 std::size_t limit_dimension) {
+    interior_simplex_cell_maps_.resize(limit_dimension + 1);
+    boundary_simplex_cell_maps_.resize(limit_dimension);
+    if (!interior_simplex_map.empty())
+      cod_d_ = interior_simplex_map.begin()->first.dimension();
+    construct_complex_(interior_simplex_map, boundary_simplex_map);
+  }
+
+  /**
+   * \brief Returns the dimension of the cell complex.
+   */
+  std::size_t intrinsic_dimension() const {
+    return intr_d_;
+  }
+
+  /**
+   * \brief Returns a vector of maps from the cells of various dimensions in the interior
+   *  of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations 
+   *  of the corresponding simplices in the ambient triangulation.
+   */
+  const Simplex_cell_maps& interior_simplex_cell_maps() const {
+    return interior_simplex_cell_maps_;
+  }
+
+  /**
+   * \brief Returns a vector of maps from the cells of various dimensions on the boundary
+   *  of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations 
+   *  of the corresponding simplices in the ambient triangulation.
+   */
+  const Simplex_cell_maps& boundary_simplex_cell_maps() const {
+    return boundary_simplex_cell_maps_;
+  }
+  
+  /**
+   * \brief Returns a map from the cells of a given dimension in the interior 
+   *  of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations 
+   *  of the corresponding simplices in the ambient triangulation.
+   *
+   * \param[in] cell_d The dimension of the cells.
+   */
+  const Simplex_cell_map& interior_simplex_cell_map(std::size_t cell_d) const {
+    return interior_simplex_cell_maps_[cell_d];
+  }
+
+  /**
+   * \brief Returns a map from the cells of a given dimension on the boundary 
+   *  of the cell complex of type Gudhi::Hasse_cell to the permutahedral representations 
+   *  of the corresponding simplices in the ambient triangulation.
+   *
+   * \param[in] cell_d The dimension of the cells.
+   */
+  const Simplex_cell_map& boundary_simplex_cell_map(std::size_t cell_d) const {
+    return boundary_simplex_cell_maps_[cell_d];
+  }
+
+  /**
+   * \brief Returns a map from the cells in the cell complex of type Gudhi::Hasse_cell
+   *  to the permutahedral representations of the corresponding simplices in the 
+   *  ambient triangulation.
+   */
+  const Cell_simplex_map& cell_simplex_map() const {
+    return cell_simplex_map_;
+  }
+
+  /**
+   * \brief Returns a map from the vertex cells in the cell complex of type Gudhi::Hasse_cell
+   *  to their Cartesian coordinates.
+   */
+  const Cell_point_map& cell_point_map() const {
+    return cell_point_map_;
+  }
+
+  /**
+   * \brief Conxtructor for the class Cell_complex.
+   *
+   * \param[in] intrinsic_dimension The dimension of the cell complex.
+   */
+  Cell_complex(std::size_t intrinsic_dimension)
+    : intr_d_(intrinsic_dimension) {}
+
+  ~Cell_complex() {
+    for (Hasse_cell* hs_ptr: hasse_cells_)
+      delete hs_ptr;
+  }
+  
+private:
+  std::size_t intr_d_, cod_d_;
+  Simplex_cell_maps interior_simplex_cell_maps_, boundary_simplex_cell_maps_;
+  Cell_simplex_map cell_simplex_map_;
+  Cell_point_map cell_point_map_;
+  std::vector<Hasse_cell*> hasse_cells_;
+};
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h b/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h
new file mode 100644
index 00000000..bd730368
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Cell_complex/Hasse_diagram_cell.h
@@ -0,0 +1,310 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Pawel Dlotko
+ *
+ *    Copyright (C) 2017 Swansea University UK
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <iostream>
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <sstream>
+#include <type_traits>
+#include <cstdlib>
+#include <cstdio>
+
+
+
+
+#ifndef HASSE_DIAGRAM_CELL_H
+#define HASSE_DIAGRAM_CELL_H
+
+
+namespace Gudhi {
+namespace Hasse_diagram {
+
+
+template < typename Cell_type > class Hasse_diagram;
+
+
+
+/**
+ * \class Hasse_diagram_cell
+ * \brief Data structure to store a cell in a Hasse diagram.
+ *
+ * \ingroup Hasse_diagram
+ *
+ * \details
+ * This is a data structure to store a cell in a general Hasse diagram data structure. 	It store the following information about the cell:
+ * References to boundary and coBoundary elements, dimension of a cell and its filtration. It also allow to store any additional
+ * information of a type Additional_information which is a template parameter of the class (set by default to void). 
+ *
+ * Please refer to \ref Hasse_diagram for examples.
+ *
+ * The complex is a template class requiring the following parameters:
+ * Incidence_type_ - determine the type of incidence coefficients. Use integers in most general case. 
+ * Filtration_type_ - type of filtration of cells.
+ * Additional_information_ (set by default to void) - allows to store any 
+ * additional information in the cells of Hasse diagrams.
+ *
+ */
+template <typename Incidence_type_, typename Filtration_type_ , typename Additional_information_ = void>
+class Hasse_diagram_cell
+{
+public:
+	typedef Incidence_type_ Incidence_type;
+	typedef Filtration_type_ Filtration_type;
+	typedef Additional_information_ Additional_information;
+	using Cell_range = std::vector< std::pair<Hasse_diagram_cell*,Incidence_type> >;
+
+    /**
+     * Default constructor.
+    **/
+	Hasse_diagram_cell():dimension(0),position(0),deleted_(false){}
+
+	/**
+     * Constructor of a cell of dimension dim.
+    **/
+	Hasse_diagram_cell( int dim ):dimension(dim),position(0),deleted_(false){}
+	
+	/**
+     * Constructor of a cell of dimension dim.
+    **/
+	Hasse_diagram_cell( int dim , Filtration_type filt_ ):dimension(dim),position(0),deleted_(false),filtration(filt_){}
+
+	/**
+     * Constructor of a cell of dimension dim with a given boundary.
+    **/	
+	Hasse_diagram_cell( const Cell_range& boundary_ , int dim ):
+	dimension(dim),boundary(boundary_),position(0),deleted_(false){}
+
+	/**
+     * Constructor of a cell of dimension dim with a given boundary and coboundary.
+    **/
+	Hasse_diagram_cell( const Cell_range&  boundary_ , const Cell_range& coboundary_,
+		 int dim ):dimension(dim),boundary(boundary_),coBoundary(coboundary_),
+		 position(0),deleted_(false){}
+
+	/**
+     * Constructor of a cell of dimension dim with a given boundary, coboundary and
+     * additional information.
+    **/
+	Hasse_diagram_cell( const Cell_range&  boundary_ , const Cell_range&  coboundary_,
+	const Additional_information& ai, int dim ):
+	dimension(dim),boundary(boundary_),coBoundary(coboundary_),additional_info(ai),
+	position(0),deleted_(false){}
+
+	/**
+     * Construcor of a cell of dimension dim having given additional information.
+    **/
+	Hasse_diagram_cell(Additional_information ai, int dim ):
+	dimension(dim),additional_info(ai),position(0),deleted_(false){}
+
+	/**
+     * Procedure to get the boundary of a fiven cell. The output format
+     * is a vector of pairs of pointers to boundary elements and incidence
+     * coefficients.
+    **/
+	inline Cell_range& get_boundary(){return this->boundary;}
+
+	/**
+     * Procedure to get the coboundary of a fiven cell. The output format
+     * is a vector of pairs of pointers to coboundary elements and incidence
+     * coefficients.
+    **/
+	inline Cell_range& get_coBoundary(){return this->coBoundary;}
+
+	/**
+     * Procedure to get the dimension of a cell.
+    **/
+	inline int& get_dimension(){return this->dimension;}
+
+	/**
+     * Procedure to get additional information about the cell.s
+    **/
+	inline Additional_information& get_additional_information(){return this->additional_info;}
+
+	/**
+	 * Procedure to retrive position of the cell in the structure. It is used in
+	 * the implementation of Hasse diagram and set by it. Note that removal of
+	 * cell and subsequent call of clean_up_the_structure will change those
+	 * positions.
+	**/
+	inline unsigned& get_position(){return this->position;}
+	
+	/**
+	 * Accessing the filtration of the cell.
+	**/
+	inline Filtration_type& get_filtration()
+	{
+		//std::cout << "Accessing the filtration of a cell : " << *this << std::endl;
+		return this->filtration;
+	}		
+
+	/**
+	 * A procedure used to check if the cell is deleted. It is used by the
+	 * subsequent implementation of Hasse diagram that is absed on lazy
+	 * delete.
+	**/
+	inline bool deleted(){return this->deleted_;}
+
+	template < typename Cell_type >
+	friend class Hasse_diagram;
+	
+	template < typename Cell_type >
+	friend class is_before_in_filtration;
+	
+	
+	template <typename Complex_type , typename Cell_type >  
+	friend std::vector<Cell_type*> convert_to_vector_of_Cell_type( Complex_type& cmplx );
+
+	/**
+	 * Procedure to remove deleted boundary and coboundary elements from the
+	 * vectors of boundary and coboundary elements of this cell.
+	**/
+	void remove_deleted_elements_from_boundary_and_coboundary()
+	{
+		Cell_range new_boundary;
+		new_boundary.reserve( this->boundary.size() );
+		for ( size_t bd = 0 ; bd != this->boundary.size() ; ++bd )
+		{
+			if ( !this->boundary[bd].first->deleted() )
+			{
+				new_boundary.push_back( this->boundary[bd] );
+			}
+		}
+		this->boundary.swap( new_boundary );
+
+		Cell_range new_coBoundary;
+		new_coBoundary.reserve( this->coBoundary.size() );
+		for ( size_t cbd = 0 ; cbd != this->coBoundary.size() ; ++cbd )
+		{
+			if ( !this->coBoundary[cbd].first->deleted() )
+			{
+				new_coBoundary.push_back( this->coBoundary[cbd] );
+			}
+		}
+		this->coBoundary.swap( new_coBoundary );
+	}
+
+	/**
+	 * Writing to a stream operator.
+	**/
+	friend std::ostream& operator<<( std::ostream& out, const Hasse_diagram_cell<Incidence_type,Filtration_type,Additional_information>& c )
+	{
+		 //cout << "position : " << c.position << ", dimension : " << c.dimension << ", filtration: " << c.filtration << ", size of boudary : " <<  c.boundary.size() << "\n";
+		 out << c.position << " " << c.dimension << " " << c.filtration << std::endl;
+		 for ( size_t bd = 0 ; bd != c.boundary.size() ; ++bd )
+	     {
+			 //do not write out the cells that has been deleted
+			 if ( c.boundary[bd].first->deleted() )continue;
+			 out <<  c.boundary[bd].first->position << " " << c.boundary[bd].second << " ";
+		 }		
+		 out << std::endl;
+		return out;
+	}
+	
+		
+	/**
+	 * Procedure that return vector of pointers to boundary elements of a given cell.
+	**/  	
+	inline std::vector< Hasse_diagram_cell* > get_list_of_boundary_elements()
+	{
+		std::vector< Hasse_diagram_cell* > result;	
+		size_t size_of_boundary = this->boundary.size();
+		result.reserve( size_of_boundary );	
+		for ( size_t bd = 0 ; bd != size_of_boundary ; ++bd )
+		{
+			result.push_back( this->boundary[bd].first );
+		}		
+		return result;
+	}
+	
+	/**
+	 * Procedure that return vector of positios of boundary elements of a given cell.
+	**/  	
+	inline std::vector< unsigned > get_list_of_positions_of_boundary_elements()
+	{
+		std::vector< unsigned > result;	
+		size_t size_of_boundary = this->boundary.size();
+		result.reserve( size_of_boundary );			
+		for ( size_t bd = 0 ; bd != size_of_boundary ; ++bd )
+		{
+			result.push_back( this->boundary[bd].first->position );
+		}		
+		return result;
+	}
+	
+	/**
+	 * Function that display a string being a signature of a structure. 
+	 * Used mainly for debugging purposes. 
+	**/ 
+	std::string full_signature_of_the_structure()
+	{
+		std::string result;
+		result += "dimension: ";
+		result += std::to_string(this->dimension);
+		result += " filtration: ";
+		result += std::to_string(this->filtration);
+		result += " position: ";
+		result += std::to_string(this->position);
+		result += " deleted_: ";
+		result += std::to_string(this->deleted_);
+		
+		//if the Additional_information is not void, add them to
+		//the signature as well.
+		if ( std::is_same<Additional_information, void>::value )
+		{			
+			result += " Additional_information: ";
+			result += std::to_string(this->additional_info);
+		}
+		result += " boundary ";
+		for ( size_t bd = 0 ; bd != this->boundary.size() ; ++bd )
+		{
+			result += "( " + std::to_string(this->boundary[bd].first->position);
+			result += " " + std::to_string(this->boundary[bd].second);
+			result += ") ";			
+		}
+		
+		result += " coBoundary ";
+		for ( size_t cbd = 0 ; cbd != this->coBoundary.size() ; ++cbd )
+		{
+			result += "( " + std::to_string(this->coBoundary[cbd].first->position);
+			result += " " + std::to_string(this->coBoundary[cbd].second);
+			result += ") ";			
+		}
+		
+		return result;
+	}
+		
+	
+protected:
+	Cell_range boundary;
+	Cell_range coBoundary;
+	int dimension;
+	Additional_information additional_info;
+	unsigned position;
+	bool deleted_;
+	Filtration_type filtration;
+
+	/**
+	 * A procedure to delete a cell. It is a private function of the Hasse_diagram_cell
+	 * class, since in the Hasse_diagram class I want to have a control
+	 * of removal of cells. Therefore, to remove cell please use
+	 * remove_cell in the Hasse_diagram structure.
+	**/
+	void delete_cell(){ this->deleted_ = true; }
+};//Hasse_diagram_cell
+
+
+
+}//namespace Hasse_diagram
+}//namespace Gudhi
+
+#endif //CELL_H
diff --git a/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h
new file mode 100644
index 00000000..bc59ac12
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Coxeter_triangulation.h
@@ -0,0 +1,94 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef COXETER_TRIANGULATION_H_
+#define COXETER_TRIANGULATION_H_
+
+#include <stack>
+#include <map>
+#include <numeric> //iota
+
+#include <boost/range/iterator_range.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+
+#include <Eigen/Eigenvalues>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+
+#include <gudhi/Freudenthal_triangulation.h>
+#include <gudhi/Permutahedral_representation.h>
+
+namespace Gudhi {
+  
+namespace coxeter_triangulation {
+
+/** 
+ * \class Coxeter_triangulation
+ * \brief A class that stores Coxeter triangulation of type \f$\tilde{A}_d\f$.
+ * This triangulation has the greatest simplex quality out of all linear transformations
+ * of the Freudenthal-Kuhn triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation.
+ *  Needs to be a model of SimplexInCoxeterTriangulation.
+ */
+template <class Permutahedral_representation_
+	  = Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
+class Coxeter_triangulation : public Freudenthal_triangulation<Permutahedral_representation_> {
+
+  typedef Eigen::MatrixXd Matrix;
+
+  Matrix root_matrix(unsigned d) {
+    Matrix cartan(d,d);
+    for (unsigned i = 0; i < d; i++) {
+      cartan(i,i) = 1.0;
+    }
+    for (unsigned i = 1; i < d; i++) {
+      cartan(i-1,i) = -0.5;
+      cartan(i,i-1) = -0.5;
+    }
+    for (unsigned i = 0; i < d; i++)
+      for (unsigned j = 0; j < d; j++)
+	if (j+1 < i || j > i+1) 
+	  cartan(i,j) = 0;    
+    Eigen::SelfAdjointEigenSolver<Matrix> saes(cartan);
+    Eigen::VectorXd sqrt_diag(d);
+    for (unsigned i = 0; i < d; ++i)
+      sqrt_diag(i) = std::sqrt(saes.eigenvalues()[i]);
+
+    Matrix lower(d,d);
+    for (unsigned i = 0; i < d; i++)
+      for (unsigned j = 0; j < d; j++)
+	if (i < j)
+	  lower(i,j) = 0;
+	else
+	  lower(i,j) = 1;    
+    Matrix result = (lower * saes.eigenvectors()*sqrt_diag.asDiagonal()).inverse();
+    return result;
+  }
+  
+public:
+
+  /** \brief Constructor of Coxeter triangulation of a given dimension. 
+   * @param[in] dimension The dimension of the triangulation.
+   */
+  Coxeter_triangulation(std::size_t dimension)
+    : Freudenthal_triangulation<Permutahedral_representation_>(dimension, root_matrix(dimension)) {}
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h
new file mode 100644
index 00000000..be5d275c
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Freudenthal_triangulation.h
@@ -0,0 +1,239 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FREUDENTHAL_TRIANGULATION_H_
+#define FREUDENTHAL_TRIANGULATION_H_
+
+#include <stack>
+#include <map>
+#include <numeric> //iota
+
+#include <Eigen/Eigenvalues>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+
+#include <gudhi/Permutahedral_representation.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Freudenthal_triangulation
+ * \brief A class that stores any affine transformation of the Freudenthal-Kuhn
+ * triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \details The data structure is a record that consists of a matrix
+ * that represents the linear transformation of the Freudenthal-Kuhn triangulation
+ * and a vector that represents the offset.
+ *
+ * \tparam Permutahedral_representation_ Type of a simplex given by a permutahedral representation.
+ * Needs to be a model of SimplexInCoxeterTriangulation.
+ */
+template <class Permutahedral_representation_
+	  = Permutahedral_representation<std::vector<int>, std::vector<std::vector<std::size_t> > > >
+class Freudenthal_triangulation {
+
+  typedef Eigen::MatrixXd Matrix;
+  typedef Eigen::VectorXd Vector;
+  typedef Eigen::SparseMatrix<double> SparseMatrix;
+  typedef Eigen::Triplet<double> Triplet;
+  
+public:
+
+  /** \brief Type of the simplices in the triangulation. */
+  typedef Permutahedral_representation_ Simplex_handle;
+
+  /** \brief Type of the vertices in the triangulation. */
+  typedef typename Permutahedral_representation_::Vertex Vertex_handle;
+  
+  
+  /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension. 
+   * @param[in] dimension The dimension of the triangulation.
+   */
+  Freudenthal_triangulation(std::size_t dimension)
+    : dimension_(dimension),
+      matrix_(Matrix::Identity(dimension, dimension)),
+      offset_(Vector::Zero(dimension)),
+      colpivhouseholderqr_(matrix_.colPivHouseholderQr()),
+      is_freudenthal(true)  {
+  }
+
+  /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under
+   * a linear transformation by a given matrix.
+   * @param[in] dimension The dimension of the triangulation.
+   * @param[in] matrix The matrix that defines the linear transformation.
+   * Needs to be invertible.
+   */
+  Freudenthal_triangulation(std::size_t dimension, const Matrix& matrix)
+    : dimension_(dimension),
+      matrix_(matrix),
+      offset_(Vector::Zero(dimension)),
+      colpivhouseholderqr_(matrix_.colPivHouseholderQr()),
+      is_freudenthal(false) {
+  }
+
+  /** \brief Constructor of the Freudenthal-Kuhn triangulation of a given dimension under
+   * an affine transformation by a given matrix and a translation vector.
+   * @param[in] dimension The dimension of the triangulation.
+   * @param[in] matrix The matrix that defines the linear transformation.
+   * Needs to be invertible.
+   * @param[in] offset The offset vector.
+   */
+  Freudenthal_triangulation(unsigned dimension, const Matrix& matrix, const Vector& offset)
+    : dimension_(dimension),
+      matrix_(matrix),
+      offset_(offset),
+      colpivhouseholderqr_(matrix_.colPivHouseholderQr()),
+      is_freudenthal(false) {
+  }
+  
+
+  /** \brief Dimension of the triangulation. */
+  unsigned dimension() const {
+    return dimension_;
+  }
+
+  /** \brief Matrix that defines the linear transformation of the triangulation. */
+  const Matrix& matrix() const {
+    return matrix_;
+  }
+
+  /** \brief Vector that defines the offset of the triangulation. */
+  const Vector& offset() const {
+    return offset_;
+  }
+
+  /** \brief Change the linear transformation matrix to a given value. 
+   *  @param[in] matrix New value of the linear transformation matrix.
+   */
+  void change_matrix(const Eigen::MatrixXd& matrix) {
+    matrix_ = matrix;
+    colpivhouseholderqr_ = matrix.colPivHouseholderQr();
+    is_freudenthal = false;
+  }  
+
+  /** \brief Change the offset vector to a given value. 
+   *  @param[in] offset New value of the offset vector.
+   */
+  void change_offset(const Eigen::VectorXd& offset) {
+    offset_ = offset;
+    is_freudenthal = false;
+  }
+  
+  /** \brief Returns the permutahedral representation of the simplex in the
+   *  triangulation that contains a given query point.
+   * \details Using the additional parameter scale, the search can be done in a 
+   * triangulation that shares the origin, but is scaled by a given factor.
+   * This parameter can be useful to simulate the point location in a subdivided
+   * triangulation.
+   * The returned simplex is always minimal by inclusion.
+   * 
+   * \tparam Point_d A class that represents a point in d-dimensional Euclidean space.
+   * The coordinates should be random-accessible. Needs to provide the method size().
+   *
+   * @param[in] point The query point.
+   * @param[in] scale The scale of the triangulation. 
+   */
+  template <class Point_d>
+  Simplex_handle locate_point(const Point_d& point, double scale=1) const {
+    using Ordered_set_partition = typename Simplex_handle::OrderedSetPartition;      
+    using Part = typename Ordered_set_partition::value_type;      
+    unsigned d = point.size();
+    assert(d == dimension_);
+    double error = 1e-9;
+    Simplex_handle output;
+    std::vector<double> z;
+    if (is_freudenthal) {
+      for (std::size_t i = 0; i < d; i++) {
+	double x_i = scale * point[i];
+	int y_i = std::floor(x_i);
+	output.vertex().push_back(y_i);
+	z.push_back(x_i - y_i);
+      }
+    }
+    else {
+      Eigen::VectorXd p_vect(d);
+      for (std::size_t i = 0; i < d; i++)
+	p_vect(i) = point[i];
+      Eigen::VectorXd x_vect = colpivhouseholderqr_.solve(p_vect - offset_);
+      for (std::size_t i = 0; i < d; i++) {
+	double x_i = scale * x_vect(i);
+	int y_i = std::floor(x_i);
+	output.vertex().push_back(y_i);
+	z.push_back(x_i - y_i);
+      }
+    }
+    z.push_back(0);
+    Part indices(d+1);
+    std::iota(indices.begin(), indices.end(), 0);
+    std::sort(indices.begin(),
+	      indices.end(),
+	      [&z](std::size_t i1, std::size_t i2) {return z[i1] > z[i2];});
+
+    output.partition().push_back(Part(1, indices[0]));
+    for (std::size_t i = 1; i <= d; ++i)
+      if (z[indices[i-1]] > z[indices[i]] + error)
+	output.partition().push_back(Part(1, indices[i]));
+      else
+	output.partition().back().push_back(indices[i]);
+    return output;
+  }
+
+  /** \brief Returns the Cartesian coordinates of the given vertex.
+   * \details Using the additional parameter scale, the search can be done in a 
+   * triangulation that shares the origin, but is scaled by a given factor.
+   * This parameter can be useful to simulate the computation of Cartesian coordinates 
+   * of a vertex in a subdivided triangulation.
+   * @param[in] vertex The query vertex.
+   * @param[in] scale The scale of the triangulation. 
+   */
+  Eigen::VectorXd cartesian_coordinates(const Vertex_handle& vertex, double scale = 1) const {
+    Eigen::VectorXd v_vect(dimension_);
+    for (std::size_t j = 0; j < dimension_; j++)
+      v_vect(j) = vertex[j] / scale;
+    return matrix_ * v_vect + offset_;
+  }
+
+  /** \brief Returns the Cartesian coordinates of the barycenter of a given simplex.
+   * \details Using the additional parameter scale, the search can be done in a 
+   * triangulation that shares the origin, but is scaled by a given factor.
+   * This parameter can be useful to simulate the computation of Cartesian coordinates 
+   * of the barycenter of a simplex in a subdivided triangulation.
+   * @param[in] simplex The query simplex.
+   * @param[in] scale The scale of the triangulation. 
+   */
+  Eigen::VectorXd barycenter(const Simplex_handle& simplex, double scale = 1) const {
+    Eigen::VectorXd res_vector(dimension_);
+    for (size_t i = 0; i < dimension_; ++i)
+      res_vector(i) = 0;
+    for (auto v: simplex.vertex_range()) {
+      res_vector += cartesian_coordinates(v, scale);
+    }
+    return (1./(simplex.dimension()+1)) * res_vector;
+  }
+  
+protected:
+  unsigned dimension_;
+  Matrix matrix_;
+  Vector offset_;
+  Eigen::ColPivHouseholderQR<Matrix> colpivhouseholderqr_;
+  bool is_freudenthal;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h
new file mode 100644
index 00000000..43198b89
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Cartesian_product.h
@@ -0,0 +1,171 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_CARTESIAN_PRODUCT_H_
+#define FUNCTIONS_CARTESIAN_PRODUCT_H_
+
+#include <cstdlib>
+#include <tuple>
+#include <utility>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* Get the domain dimension of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof... (T), std::size_t>::type
+get_amb_d (const std::tuple<T...>& tuple) {
+  return 0;
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof... (T), std::size_t>::type
+get_amb_d (const std::tuple<T...>& tuple) {
+  return std::get<I>(tuple).amb_d() + get_amb_d<I+1, T...>(tuple);
+}
+
+/* Get the codomain dimension of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof... (T), std::size_t>::type
+get_cod_d (const std::tuple<T...>& tuple) {
+  return 0;
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof... (T), std::size_t>::type
+get_cod_d (const std::tuple<T...>& tuple) {
+  return std::get<I>(tuple).cod_d() + get_cod_d<I+1, T...>(tuple);
+}
+
+/* Get the seed of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof... (T), void>::type
+get_seed (const std::tuple<T...>& tuple, Eigen::VectorXd& point, std::size_t i = 0) {
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof... (T), void>::type
+get_seed (const std::tuple<T...>& tuple, Eigen::VectorXd& point, std::size_t i = 0) {
+  const auto& f = std::get<I>(tuple);
+  std::size_t n = f.amb_d();
+  Eigen::VectorXd seed = f.seed();
+  for (std::size_t j = 0; j < n; ++j)
+    point(i+j) = seed(j);
+  get_seed<I+1, T...>(tuple, point, i+n);  
+}
+
+/* Get the seed of the tuple of functions.
+ */
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I == sizeof... (T), void>::type
+get_value (const std::tuple<T...>& tuple, const Eigen::VectorXd& x, Eigen::VectorXd& point, std::size_t i = 0, std::size_t j = 0) {
+}
+
+template <std::size_t I = 0, typename... T>
+inline typename std::enable_if<I != sizeof... (T), void>::type
+get_value (const std::tuple<T...>& tuple, const Eigen::VectorXd& x, Eigen::VectorXd& point, std::size_t i = 0, std::size_t j = 0) {
+  const auto& f = std::get<I>(tuple);
+  std::size_t n = f.amb_d();
+  std::size_t k = f.cod_d();
+  Eigen::VectorXd x_i(n);
+  for (std::size_t l = 0; l < n; ++l)
+    x_i(l) = x(i+l);
+  Eigen::VectorXd res = f(x_i);
+  for (std::size_t l = 0; l < k; ++l)
+    point(j+l) = res(l);
+  get_value<I+1, T...>(tuple, x, point, i+n, j+k);  
+}
+
+
+/** 
+ * \class Cartesian_product
+ * \brief Constructs the function the zero-set of which is the Cartesian product
+ * of the zero-sets of some given functions.
+ *
+ * \tparam Functions A pack template parameter for functions. All functions should be models of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class... Functions>
+struct Cartesian_product : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result(cod_d_);
+    get_value(function_tuple_, p, result, 0, 0);
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return amb_d_;}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return cod_d_;}
+
+  /** \brief Returns a point on the zero-set. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result(amb_d_);
+    get_seed(function_tuple_, result, 0);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the Cartesian product function.
+   *
+   * @param[in] functions The functions the zero-sets of which are factors in the
+   * Cartesian product of the resulting function.
+   */
+  Cartesian_product(const Functions&... functions)
+    : function_tuple_(std::make_tuple(functions...)) {
+    amb_d_ = get_amb_d(function_tuple_);
+    cod_d_ = get_cod_d(function_tuple_);
+  }
+
+private:
+  std::tuple<Functions...> function_tuple_;
+  std::size_t amb_d_, cod_d_;
+};
+
+
+/** 
+ * \brief Static constructor of a Cartesian product function.
+ *
+ * @param[in] functions The functions the zero-sets of which are factors in the
+ * Cartesian product of the resulting function.
+ *
+ * \tparam Functions A pack template parameter for functions. All functions should be models of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <typename... Functions>
+Cartesian_product<Functions...> make_product_function(const Functions&... functions) {
+  return Cartesian_product<Functions...>(functions...);
+}
+
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h
new file mode 100644
index 00000000..ea354fac
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Constant_function.h
@@ -0,0 +1,71 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_CONSTANT_FUNCTION_H_
+#define FUNCTIONS_CONSTANT_FUNCTION_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Constant_function 
+ * \brief A class that encodes a constant function from R^d to R^k.
+ * This class does not have any implicit manifold in correspondence.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Constant_function : public Function {
+
+  /** \brief Value of the function at a specified point. The value is constant.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result = value_;
+    return result;
+  }
+  
+  /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+  std::size_t amb_d() const {return d_;};
+
+  /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+  std::size_t cod_d() const {return k_;};
+
+  /** \brief No seed point is available. Throws an exception on evocation. */
+  Eigen::VectorXd seed() const {
+    throw "Seed invoked on a constant function.\n";
+  }
+
+  Constant_function() {}
+  
+  /** 
+   * \brief Constructor of a constant function from R^d to R^m.
+   *
+   * @param[in] d The domain dimension.
+   * @param[in] k The codomain dimension.
+   * @param[in] value The constant value of the function.
+   */
+  Constant_function(std::size_t d, std::size_t k, const Eigen::VectorXd& value)
+    : d_(d), k_(k), value_(value) {}
+  
+  std::size_t d_, k_;
+  Eigen::VectorXd value_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h
new file mode 100644
index 00000000..9476944b
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Embed_in_Rd.h
@@ -0,0 +1,106 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_EMBED_IN_RD_H_
+#define FUNCTIONS_EMBED_IN_RD_H_
+
+#include <cstdlib>
+#include <random>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ * \class Embed_in_Rd
+ * \brief Embedding of an implicit manifold in a higher dimension.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+struct Embed_in_Rd : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd x = p;
+    Eigen::VectorXd x_k(fun_.amb_d()), x_rest(d_ - fun_.amb_d());
+    for (std::size_t i = 0; i < fun_.amb_d(); ++i)
+      x_k(i) = x(i);
+    for (std::size_t i = fun_.amb_d(); i < d_; ++i)
+      x_rest(i - fun_.amb_d()) = x(i);
+    Eigen::VectorXd result = fun_(x_k);
+    result.conservativeResize(this->cod_d());
+    for (std::size_t i = fun_.cod_d(); i < this->cod_d(); ++i)
+      result(i) = x_rest(i - fun_.cod_d());
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return d_;}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return d_-(fun_.amb_d() - fun_.cod_d());}
+
+  /** \brief Returns a point on the zero-set of the embedded function. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = fun_.seed();
+    result.conservativeResize(d_);
+    for (std::size_t l = fun_.amb_d(); l < d_; ++l)
+      result(l) = 0;
+    return result;
+  }
+
+  /** 
+   * \brief Constructor of the embedding function.
+   *
+   * @param[in] function The function to be embedded in higher dimension.
+   * @param[in] d Embedding dimension.
+   */
+  Embed_in_Rd(const Function_& function, std::size_t d) :
+    fun_(function), d_(d) {
+  }
+  Function_ fun_;
+  std::size_t d_;
+};
+
+
+/** 
+ * \brief Static constructor of an embedding function.
+ *
+ * @param[in] function The function to be embedded in higher dimension.
+ * @param[in] d Embedding dimension.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Embed_in_Rd<Function_> make_embedding(const Function_& function, std::size_t d) {
+  return Embed_in_Rd<Function_>(function, d);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function.h
new file mode 100644
index 00000000..dbcb0142
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function.h
@@ -0,0 +1,54 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_H_
+#define FUNCTIONS_FUNCTION_H_
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function
+ * \brief The parent class for all functions implemented in the module.
+ *  Contains virtual methods needed to be a model of the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function {
+  
+  /** \brief Virtual method for the value of the function at a specified point.
+   * @param[in] p The input point.
+   */
+  virtual Eigen::VectorXd evaluate(const Eigen::VectorXd& p) const {return Eigen::VectorXd(0);}
+  
+  /** \brief Virtual method for the domain dimension. */
+  virtual std::size_t amb_d() const {return 0;};
+
+  /** \brief Virtual method for the codomain dimension. */
+  virtual std::size_t cod_d() const {return 0;};
+
+  /** \brief Virtual method for the seed point. */
+  virtual Eigen::VectorXd seed() const {return Eigen::VectorXd(0);}
+
+  /** \brief Virtual destructor. */
+  virtual ~Function() {}
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h
new file mode 100644
index 00000000..a7d4a965
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_Sm_in_Rd.h
@@ -0,0 +1,131 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_SM_IN_RD_H_
+#define FUNCTIONS_FUNCTION_SM_IN_RD_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_Sm_in_Rd 
+ * \brief A class for the function that defines an m-dimensional implicit sphere embedded
+ * in the d-dimensional Euclidean space.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_Sm_in_Rd: public Function {
+
+/** \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd x = p;
+    for (std::size_t i = 0; i < d_; ++i)
+      x(i) -= center_[i];
+    Eigen::VectorXd result = Eigen::VectorXd::Zero(k_); 
+    for (std::size_t i = 0; i < m_+1; ++i)
+      result(0) += x(i)*x(i);
+    result(0) -= r_*r_;
+    for (std::size_t j = 1; j < k_; ++j)
+      result(j) = x(m_+j);
+    return result;
+  }
+  
+  /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+  std::size_t amb_d() const {return d_;};
+
+  /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+  std::size_t cod_d() const {return k_;};
+
+  /** \brief Returns a point on the sphere. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = Eigen::VectorXd::Zero(d_);
+    result(0) += r_;
+    for (std::size_t i = 0; i < d_; ++i)
+      result(i) += center_[i];
+    return result;
+  }
+    
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+   * in the d-dimensional Euclidean space.
+   *
+   * @param[in] r The radius of the sphere.
+   * @param[in] m The dimension of the sphere.
+   * @param[in] d The ambient dimension of the sphere.
+   * @param[in] center The center of the sphere.
+   */
+  Function_Sm_in_Rd(double r,
+		    std::size_t m,
+		    std::size_t d,
+		    Eigen::VectorXd center)
+    : m_(m), k_(d-m), d_(d), r_(r), center_(center) {}
+
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+   * in the d-dimensional Euclidean space centered at the origin.
+   *
+   * @param[in] r The radius of the sphere.
+   * @param[in] m The dimension of the sphere.
+   * @param[in] d The ambient dimension of the sphere.
+   */
+  Function_Sm_in_Rd(double r,
+		    std::size_t m,
+		    std::size_t d)
+    : m_(m), k_(d-m), d_(d), r_(r), center_(Eigen::VectorXd::Zero(d_)) {}
+
+  
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+   * in the (m+1)-dimensional Euclidean space.
+   *
+   * @param[in] r The radius of the sphere.
+   * @param[in] m The dimension of the sphere.
+   * @param[in] center The center of the sphere.
+   */
+  Function_Sm_in_Rd(double r,
+		    std::size_t m,
+		    Eigen::VectorXd center)
+    : m_(m), k_(1), d_(m_+1), r_(r), center_(center) {}
+
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit sphere embedded
+   * in the (m+1)-dimensional Euclidean space centered at the origin.
+   *
+   * @param[in] r The radius of the sphere.
+   * @param[in] m The dimension of the sphere.
+   */
+  Function_Sm_in_Rd(double r,
+		    std::size_t m)
+    : m_(m), k_(1), d_(m_+1), r_(r), center_(Eigen::VectorXd::Zero(d_)) {}
+
+  Function_Sm_in_Rd(const Function_Sm_in_Rd& rhs)
+    : Function_Sm_in_Rd(rhs.r_, rhs.m_, rhs.d_, rhs.center_) {}
+
+  
+protected:
+  std::size_t m_, k_, d_;
+  double r_;
+  Eigen::VectorXd center_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h
new file mode 100644
index 00000000..cb3af848
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_affine_plane_in_Rd.h
@@ -0,0 +1,101 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_
+#define FUNCTIONS_FUNCTION_AFFINE_PLANE_IN_RD_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_affine_plane_in_Rd 
+ * \brief A class for the function that defines an m-dimensional implicit affine plane 
+ * embedded in d-dimensional Euclidean space.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_affine_plane_in_Rd : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result = normal_matrix_.transpose() * (p - off_);
+    return result;
+  }
+  
+  /** \brief Returns the domain dimension. Same as the ambient dimension of the sphere. */
+  std::size_t amb_d() const {return d_;};
+
+  /** \brief Returns the codomain dimension. Same as the codimension of the sphere. */
+  std::size_t cod_d() const {return k_;};
+
+  /** \brief Returns a point on the affine plane. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = off_;
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit affine
+   * plane in the d-dimensional Euclidean space.
+   *
+   * @param[in] normal_matrix A normal matrix of the affine plane. The number of rows should
+   * correspond to the ambient dimension, the number of columns should corespond to 
+   * the size of the normal basis (codimension). 
+   * @param[in] offset The offset vector of the affine plane.
+   * The dimension of the vector should be the ambient dimension of the manifold.
+   */
+  Function_affine_plane_in_Rd(const Eigen::MatrixXd& normal_matrix,
+			      const Eigen::VectorXd& offset)
+    : normal_matrix_(normal_matrix),
+      d_(normal_matrix.rows()),
+      k_(normal_matrix.cols()),
+      m_(d_ - k_),
+      off_(offset) {
+    normal_matrix_.colwise().normalize();
+  }
+
+  /** 
+   * \brief Constructor of the function that defines an m-dimensional implicit affine
+   * plane in the d-dimensional Euclidean space that passes through origin.
+   *
+   * @param[in] normal_matrix A normal matrix of the affine plane. The number of rows should
+   * correspond to the ambient dimension, the number of columns should corespond to 
+   * the size of the normal basis (codimension). 
+   */
+  Function_affine_plane_in_Rd(const Eigen::MatrixXd& normal_matrix)
+    : normal_matrix_(normal_matrix),
+      d_(normal_matrix.rows()),
+      k_(normal_matrix.cols()),
+      m_(d_ - k_),
+      off_(Eigen::VectorXd::Zero(d_)) {
+    normal_matrix_.colwise().normalize();
+  }
+
+protected:
+  Eigen::MatrixXd normal_matrix_;
+  std::size_t d_, k_, m_;
+  Eigen::VectorXd off_;  
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h
new file mode 100644
index 00000000..2a76e631
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_chair_in_R3.h
@@ -0,0 +1,87 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_
+#define FUNCTIONS_FUNCTION_CHAIR_IN_R3_H_
+
+#include <cstdlib>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_chair_in_R3 
+ * \brief A class that encodes the function, the zero-set of which is a so-called
+ * "chair" surface embedded in R^3.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_chair_in_R3 : public Function {
+
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    double x = p(0)-off_[0], y = p(1)-off_[1], z = p(2)-off_[2];
+    Eigen::VectorXd result(cod_d());
+    result(0) = std::pow(x*x + y*y + z*z - a_*k_*k_, 2) - b_*((z-k_)*(z-k_) - 2*x*x)*((z+k_)*(z+k_) - 2*y*y);
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return 3;}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return 1;}
+
+  /** \brief Returns a point on the surface. */
+  Eigen::VectorXd seed() const {
+    double t1 = a_-b_;
+    double discr = t1*t1 - (1.0 - b_)*(a_*a_ - b_);
+    double z0 = k_*std::sqrt((t1+std::sqrt(discr))/(1-b_));
+    Eigen::Vector3d result(off_[0], off_[1], z0+off_[2]);
+    return result;
+  }
+
+  /** 
+   * \brief Constructor of the function that defines the 'chair' surface
+   * embedded in R^3.
+   *
+   * @param[in] a A numerical parameter.
+   * @param[in] b A numerical parameter.
+   * @param[in] k A numerical parameter.
+   * @param[in] off Offset vector.
+   */
+  Function_chair_in_R3(double a = 0.8,
+		       double b = 0.4,
+		       double k = 1.0,
+		       Eigen::Vector3d off = Eigen::Vector3d::Zero()) :
+    a_(a), b_(b), k_(k), off_(off) {}
+  
+protected:
+  double a_, b_, k_;
+  Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
+
+// (x^2 + y^2 + z^2 - a*k^2)^2 - b*((z-k)^2 - 2*x^2)*((z+k)^2 - 2*y^2)
+// sqrt(k/(1-b))*sqrt(a-b + sqrt((a-b)^2 - (1-b)*(a^2 - b)*k^2))
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h
new file mode 100644
index 00000000..af323c94
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_iron_in_R3.h
@@ -0,0 +1,73 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_IRON_IN_R3_H_
+#define FUNCTIONS_FUNCTION_IRON_IN_R3_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_iron_in_R3 
+ * \brief A class that encodes the function, the zero-set of which is a surface
+ * embedded in R^3 that ressembles an iron.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_iron_in_R3 : public Function {
+
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    double x = p(0), y = p(1), z = p(2);
+    Eigen::VectorXd result(cod_d());
+    result(0) = - std::pow(x,6)/300. - std::pow(y,6)/300. - std::pow(z,6)/300. + x*y*y*z/2.1 + y*y + std::pow(z-2, 4) - 1;
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return 3;};
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return 1;};
+
+  /** \brief Returns a point on the surface. */
+  Eigen::VectorXd seed() const {
+    Eigen::Vector3d result(std::pow(4500, 1./6), 0, 0);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines a surface embedded in R^3
+   * that ressembles an iron.
+   *
+   * @param[in] off Offset vector.
+   */
+  Function_iron_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero()) :
+    off_(off) {}
+  
+protected:
+  Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h
new file mode 100644
index 00000000..cd03a0a5
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_lemniscate_revolution_in_R3.h
@@ -0,0 +1,90 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_
+#define FUNCTIONS_FUNCTION_LEMNISCATE_REVOLUTION_IN_R3_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_lemniscate_revolution_in_R3 
+ * \brief A class that encodes the function, the zero-set of which is a surface of revolution
+ * around the x axis based on the lemniscate of Bernoulli embedded in R^3.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_lemniscate_revolution_in_R3 : public Function {
+
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    double x = p(0)-off_[0], y = p(1)-off_[1], z = p(2)-off_[2];
+    Eigen::VectorXd result(cod_d());
+    double x2 = x*x, y2 = y*y, z2 = z*z, a2 = a_*a_;
+    double t1 = x2 + y2 + z2;
+    result(0) = t1*t1 - 2*a2*(x2 - y2 - z2);
+    return result;
+  }
+
+  /** \brief Returns the (ambient) domain dimension.*/
+  std::size_t amb_d() const {return 3;};
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return 1;};
+
+  /** \brief Returns a point on the surface. This seed point is only one of 
+   * two necessary seed points for the manifold tracing algorithm.
+   * See the method seed2() for the other point. 
+   */
+  Eigen::VectorXd seed() const {
+    Eigen::Vector3d result(sqrt(2*a_)+off_[0], off_[1], off_[2]);
+    return result;
+  }
+
+  /** \brief Returns a point on the surface. This seed point is only one of 
+   * two necessary seed points for the manifold tracing algorithm. 
+   * See the method seed() for the other point. 
+   */
+  Eigen::VectorXd seed2() const {
+    Eigen::Vector3d result(-sqrt(2*a_)+off_[0], off_[1], off_[2]);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines a surface of revolution
+   * around the x axis based on the lemniscate of Bernoulli embedded in R^3.
+   *
+   * @param[in] a A numerical parameter.
+   * @param[in] off Offset vector.
+   */
+  Function_lemniscate_revolution_in_R3(double a = 1,
+				       Eigen::Vector3d off = Eigen::Vector3d::Zero())
+    : a_(a), off_(off) {}
+
+protected:
+  double a_;
+  Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h
new file mode 100644
index 00000000..91eb016b
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_moment_curve_in_Rd.h
@@ -0,0 +1,89 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_
+#define FUNCTIONS_FUNCTION_MOMENT_CURVE_IN_RD_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_moment_curve_in_Rd 
+ * \brief A class for the function that defines an implicit moment curve
+ * in the d-dimensional Euclidean space.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_moment_curve_in_Rd : public Function {
+
+/** \brief Value of the function at a specified point.
+ * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+ */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result(k_);
+    for (std::size_t i = 1; i < d_; ++i)
+      result(i-1) = p(i) - p(0) * p(i-1);
+    return result;
+  }
+  
+  /** \brief Returns the domain (ambient) dimension.. */
+  std::size_t amb_d() const {return d_;};
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return k_;};
+
+  /** \brief Returns a point on the moment curve. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = Eigen::VectorXd::Zero(d_);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines an implicit moment curve
+   * in the d-dimensional Euclidean space.
+   *
+   * @param[in] r Numerical parameter.
+   * @param[in] d The ambient dimension.
+   */
+  Function_moment_curve_in_Rd(double r,
+			      std::size_t d)
+    : m_(1), k_(d-1), d_(d), r_(r) {}
+
+  /** 
+   * \brief Constructor of the function that defines an implicit moment curve
+   * in the d-dimensional Euclidean space.
+   *
+   * @param[in] r Numerical parameter.
+   * @param[in] d The ambient dimension.
+   * @param[in] offset The offset of the moment curve.
+   */
+  Function_moment_curve_in_Rd(double r,
+			      std::size_t d,
+			      Eigen::VectorXd& offset)
+    : m_(1), k_(d-1), d_(d), r_(r), off_(offset) {}
+  
+protected:
+  std::size_t m_, k_, d_;
+  double r_;
+  Eigen::VectorXd off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h
new file mode 100644
index 00000000..3cda2518
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_torus_in_R3.h
@@ -0,0 +1,77 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_TORUS_IN_R3_H_
+#define FUNCTIONS_FUNCTION_TORUS_IN_R3_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_torus_in_R3 
+ * \brief A class that encodes the function, the zero-set of which is a torus
+ * surface embedded in R^3.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_torus_in_R3 : public Function {
+
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    double x = p(0)-off_[0], y = p(1)-off_[1], z = p(2)-off_[2];
+    Eigen::VectorXd result(cod_d());
+    result(0) = (z*z + (std::sqrt(x*x + y*y) - r_)*(std::sqrt(x*x + y*y) - r_) - R_*R_);
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return 3;};
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return 1;};
+
+  /** \brief Returns a point on the surface. */
+  Eigen::VectorXd seed() const {
+    Eigen::Vector3d result(R_ + r_ +off_[0], off_[1], off_[2]);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines a torus embedded in R^3.
+   *
+   * @param[in] R The outer radius of the torus.
+   * @param[in] r The inner radius of the torus.
+   * @param[in] off Offset vector.
+   */
+  Function_torus_in_R3(double R = 1,
+		       double r = 0.5,
+		       Eigen::Vector3d off = Eigen::Vector3d::Zero()) :
+    R_(R), r_(r), off_(off) {}
+  
+protected:
+  double R_, r_;
+  Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h
new file mode 100644
index 00000000..3d825f8f
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Function_whitney_umbrella_in_R3.h
@@ -0,0 +1,82 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_
+#define FUNCTIONS_FUNCTION_WHITNEY_UMBRELLA_IN_R3_H_
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Function_whitney_umbrella_in_R3 
+ * \brief A class that encodes the function, the zero-set of which is the Whitney umbrella
+ * surface embedded in R^3.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Function_whitney_umbrella_in_R3 : public Function {
+
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    double x = p(0)-off_[0], y = p(1)-off_[1], z = p(2)-off_[2];
+    Eigen::VectorXd result(cod_d());
+    result(0) = x*x - y*y*z;
+    return result;
+  }
+
+  /** \brief Returns the (ambient) domain dimension.*/
+  std::size_t amb_d() const {return 3;};
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return 1;};
+
+  /** \brief Returns a point on the surface. This seed point is only one of 
+   * two necessary seed points for the manifold tracing algorithm.
+   * See the method seed2() for the other point. 
+   */
+  Eigen::VectorXd seed() const {
+    Eigen::Vector3d result(1+off_[0], 1+off_[1], 1+off_[2]);
+    return result;
+  }
+
+  /** \brief Returns a point on the surface. This seed point is only one of 
+   * two necessary seed points for the manifold tracing algorithm. 
+   * See the method seed() for the other point. 
+   */
+  Eigen::VectorXd seed2() const {
+    Eigen::Vector3d result(-1+off_[0], -1+off_[1], 1+off_[2]);
+    return result;
+  }
+  
+  /** 
+   * \brief Constructor of the function that defines the Whitney umbrella in R^3.
+   *
+   * @param[in] off Offset vector.
+   */
+  Function_whitney_umbrella_in_R3(Eigen::Vector3d off = Eigen::Vector3d::Zero())
+    : off_(off) {}
+  Eigen::Vector3d off_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h
new file mode 100644
index 00000000..4cb5ca84
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Linear_transformation.h
@@ -0,0 +1,96 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_LINEAR_TRANSFORMATION_H_
+#define FUNCTIONS_LINEAR_TRANSFORMATION_H_
+
+#include <cstdlib>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Linear_transformation
+ * \brief Transforms the zero-set of the function by a given linear transformation.
+ * The underlying function corresponds to f(M*x), where M is the transformation matrix.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+struct Linear_transformation : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result = fun_(matrix_.householderQr().solve(p));
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return fun_.amb_d();}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return fun_.cod_d();}
+
+  /** \brief Returns a point on the zero-set. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = fun_.seed();
+    result = matrix_ * result;
+    return result;
+  }
+
+  /** 
+   * \brief Constructor of a linearly transformed function.
+   *
+   * @param[in] function The function to be linearly transformed.
+   * @param[in] matrix The transformation matrix. Its dimension should be d*d,
+   * where d is the domain (ambient) dimension of 'function'.
+   */
+  Linear_transformation(const Function_& function, const Eigen::MatrixXd& matrix) :
+    fun_(function), matrix_(matrix) {
+  }
+
+private:
+  Function_ fun_;
+  Eigen::MatrixXd matrix_;
+};
+
+
+/** 
+ * \brief Static constructor of a linearly transformed function.
+ *
+ * @param[in] function The function to be linearly transformed.
+ * @param[in] matrix The transformation matrix. Its dimension should be d*d,
+ * where d is the domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ */
+template <class Function_>
+Linear_transformation<Function_> make_linear_transformation(const Function_& function,
+						       const Eigen::MatrixXd& matrix) {
+  return Linear_transformation<Function_>(function, matrix); 
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h
new file mode 100644
index 00000000..6b7feff5
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Negation.h
@@ -0,0 +1,92 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_NEGATION_H_
+#define FUNCTIONS_NEGATION_H_
+
+#include <cstdlib>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ *\class Negation
+ * \brief Constructs the "minus" function. The zero-set is the same, but
+ * the values at other points are the negative of their original value.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+struct Negation : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result = -fun_(p);
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return fun_.amb_d();}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return fun_.cod_d();}
+
+  /** \brief Returns a point on the zero-set. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = fun_.seed();
+    return result;
+  }
+
+  /** 
+   * \brief Constructor of the negative function.
+   *
+   * @param[in] function The function to be negated.
+   */
+  Negation(const Function_& function) :
+    fun_(function) {
+  }
+  Function_ fun_;
+};
+
+
+/** 
+ * \brief Static constructor of the negative function.
+ *
+ * @param[in] function The function to be translated.
+ * @param[in] off The offset vector. The dimension should correspond to the 
+ * domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Negation<Function_>  negation(const Function_& function) {
+  return Negation<Function_>(function);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h
new file mode 100644
index 00000000..25c8a139
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/PL_approximation.h
@@ -0,0 +1,125 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_PL_APPROXIMATION_H_
+#define FUNCTIONS_PL_APPROXIMATION_H_
+
+#include <cstdlib>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class PL_approximation
+ * \brief Constructs a piecewise-linear approximation of a function induced by
+ * an ambient triangulation.
+ *
+ * \tparam Function The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ * \tparam Triangulation The triangulation template parameter. Should be a model of
+ * the concept TriangulationForManifoldTracing.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_,
+	  class Triangulation_>
+struct PL_approximation : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    std::size_t cod_d = this->cod_d();
+    std::size_t amb_d = this->amb_d();
+    auto s = tr_.locate_point(p);
+    Eigen::MatrixXd matrix(cod_d, s.dimension() + 1);
+    Eigen::MatrixXd vertex_matrix(amb_d + 1, s.dimension() + 1);
+    for (std::size_t i = 0; i < s.dimension() + 1; ++i)
+      vertex_matrix(0, i) = 1;
+    std::size_t j = 0;
+    for (auto v: s.vertex_range()) {
+      Eigen::VectorXd pt_v = tr_.cartesian_coordinates(v);
+      Eigen::VectorXd fun_v = fun_(pt_v);
+      for (std::size_t i = 1; i < amb_d + 1; ++i)
+	vertex_matrix(i, j) = pt_v(i-1); 
+      for (std::size_t i = 0; i < cod_d; ++i)
+	matrix(i, j) = fun_v(i);
+      j++;
+    }
+    assert(j == s.dimension()+1);
+    Eigen::VectorXd z(amb_d + 1);
+    z(0) = 1;
+    for (std::size_t i = 1; i < amb_d + 1; ++i)
+      z(i) = p(i-1);
+    Eigen::VectorXd lambda = vertex_matrix.colPivHouseholderQr().solve(z);    
+    Eigen::VectorXd result = matrix * lambda;
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return fun_.amb_d();}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return fun_.cod_d();}
+
+  /** \brief Returns a point on the zero-set. */
+  Eigen::VectorXd seed() const {
+    // TODO: not finished. Should use an oracle.
+    return Eigen::VectorXd(amb_d());
+  }
+
+  /** 
+   * \brief Constructor of the piecewise-linear approximation of a function 
+   * induced by an ambient triangulation.
+   *
+   * @param[in] function The function.
+   * @param[in] triangulation The ambient triangulation.
+   */
+  PL_approximation(const Function_& function, const Triangulation_& triangulation)
+    : fun_(function), tr_(triangulation) {}
+
+private:  
+  Function_ fun_;
+  Triangulation_ tr_;
+};
+
+
+/** 
+ * \brief Static constructor of the piecewise-linear approximation of a function 
+ * induced by an ambient triangulation.
+ *
+ * @param[in] function The function.
+ * @param[in] triangulation The ambient triangulation.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_,
+	  class Triangulation_>
+PL_approximation<Function_, Triangulation_>
+make_pl_approximation(const Function_& function,
+		      const Triangulation_& triangulation) {
+  return PL_approximation<Function_, Triangulation_>(function, triangulation);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h
new file mode 100644
index 00000000..6edac28c
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/Translate.h
@@ -0,0 +1,96 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_TRANSLATE_H_
+#define FUNCTIONS_TRANSLATE_H_
+
+#include <cstdlib>
+
+#include <gudhi/Functions/Function.h>
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Translate
+ * \brief Translates the zero-set of the function by a vector.
+ * The underlying function corresponds to f(x-off), where off is the offset vector.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+struct Translate : public Function {
+  
+  /** 
+   * \brief Value of the function at a specified point.
+   * @param[in] p The input point. The dimension needs to coincide with the ambient dimension.
+   */
+  Eigen::VectorXd operator()(const Eigen::VectorXd& p) const {
+    Eigen::VectorXd result = fun_(p - off_);
+    return result;
+  }
+
+  /** \brief Returns the domain (ambient) dimension. */
+  std::size_t amb_d() const {return fun_.amb_d();}
+
+  /** \brief Returns the codomain dimension. */
+  std::size_t cod_d() const {return fun_.cod_d();}
+
+  /** \brief Returns a point on the zero-set. */
+  Eigen::VectorXd seed() const {
+    Eigen::VectorXd result = fun_.seed();
+    result += off_;
+    return result;
+  }
+
+  /** 
+   * \brief Constructor of the translated function.
+   *
+   * @param[in] function The function to be translated.
+   * @param[in] off The offset vector. The dimension should correspond to the 
+   * domain (ambient) dimension of 'function'.
+   */
+  Translate(const Function_& function, const Eigen::VectorXd& off) :
+    fun_(function), off_(off) {
+  }
+  Function_ fun_;
+  Eigen::VectorXd off_;
+};
+
+
+/** 
+ * \brief Static constructor of a translated function.
+ *
+ * @param[in] function The function to be translated.
+ * @param[in] off The offset vector. The dimension should correspond to the 
+ * domain (ambient) dimension of 'function'.
+ *
+ * \tparam Function_ The function template parameter. Should be a model of 
+ * the concept FunctionForImplicitManifold.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Function_>
+Translate<Function_> translate(const Function_& function, Eigen::VectorXd off) {
+  return Translate<Function_>(function, off);
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h
new file mode 100644
index 00000000..cf0b00ac
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Functions/random_orthogonal_matrix.h
@@ -0,0 +1,70 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
+#define FUNCTIONS_RANDOM_ORTHOGONAL_MATRIX_H_
+
+#include <cstdlib>
+#include <random>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SVD>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/point_generators_d.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief Generates a uniform random orthogonal matrix using the "subgroup algorithm" by 
+ * Diaconis & Shashahani. 
+ * \details Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Uniform_random_rotation_matrices.
+ * The idea: take a random rotation matrix of dimension d-1, embed it 
+ * as a d*d matrix M with the last column (0,...,0,1).
+ * Pick a random vector v on a sphere S^d. rotate the matrix M so that its last column is v.
+ * The determinant of the matrix can be either 1 or -1
+ */
+// Note: the householderQR operation at the end seems to take a lot of time at compilation.
+// The CGAL headers are another source of long compilation time.
+Eigen::MatrixXd random_orthogonal_matrix(std::size_t d) {
+  typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> Kernel;
+  typedef typename Kernel::Point_d Point_d;
+  if (d == 1)
+    return Eigen::VectorXd::Constant(1, 1.0);
+  if (d == 2) {
+    double X = 2 * 3.14159265358;
+    double alpha = static_cast <float> (rand()) / (static_cast <float> (RAND_MAX/X));
+    Eigen::Matrix2d rot;
+    rot << cos(alpha), -sin(alpha), sin(alpha), cos(alpha);
+    return rot;
+  }
+  Eigen::MatrixXd low_dim_rot = random_orthogonal_matrix(d-1);
+  Eigen::MatrixXd rot(d,d);
+  Point_d v = *CGAL::Random_points_on_sphere_d<Point_d>(d, 1);
+  for (std::size_t i = 0; i < d; ++i)
+    rot(i,0) = v[i];
+  for (std::size_t i = 0; i < d-1; ++i)
+    for (std::size_t j = 1; j < d-1; ++j)
+      rot(i,j) = low_dim_rot(i,j-1);
+  for (std::size_t j = 1; j < d; ++j)
+    rot(d-1,j) = 0;
+  rot = rot.householderQr().householderQ(); // a way to do Gram-Schmidt, see https://forum.kde.org/viewtopic.php?f=74&t=118568#p297246
+  return rot;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
new file mode 100644
index 00000000..b1cc0fe5
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/Mesh_medit.h
@@ -0,0 +1,57 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_MESH_MEDIT_H_
+#define IO_MESH_MEDIT_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* \class Mesh_medit
+ * \brief Structure to store a mesh that can be output in Medit .mesh file format
+ *  using the output_meshes_to_medit method.
+ *
+ * \ingroup coxeter_triangulation
+ */
+struct Mesh_medit {
+  /** \brief Type of a range of vertices. */
+  typedef std::vector<Eigen::VectorXd> Vertex_points;
+  /** \brief Type of a mesh element.
+   *  A pair consisting of a vector of vertex indices of type std::size_t
+   *  and of an integer that represents the common reference number for
+   *  the mesh elements of this type. */
+  typedef std::pair<std::vector<std::size_t>, std::size_t> Mesh_element;
+  /** \brief Type of a range of mesh elements. */
+  typedef std::vector<Mesh_element> Mesh_elements;
+  /** \brief Type of a range of scalar field . */
+  typedef std::vector<double> Scalar_field_range;
+
+  /** \brief Range of vertices of type Eigen::VectorXd to output. */
+  Vertex_points vertex_points;
+  /** \brief Range of edges. */
+  Mesh_elements edges;
+  /** \brief Range of triangles. */
+  Mesh_elements triangles;
+  /** \brief Range of tetrahedra. */
+  Mesh_elements tetrahedra;
+  /** \brief Range of scalar values over triangles. */
+  Scalar_field_range triangles_scalar_range;
+  /** \brief Range of scalar values over tetrahedra. */
+  Scalar_field_range tetrahedra_scalar_range;
+};
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
new file mode 100644
index 00000000..71070c05
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/build_mesh_from_cell_complex.h
@@ -0,0 +1,176 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+#define IO_BUILD_MESH_FROM_CELL_COMPLEX_H_
+
+#include <gudhi/IO/Mesh_medit.h>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+struct Configuration {
+  bool toggle_edges = true,
+    toggle_triangles = true,
+    toggle_tetrahedra = true;
+  std::size_t ref_edges = 1,
+    ref_triangles = 1,
+    ref_tetrahedra = 1;
+
+  Configuration(bool t_edges, bool t_triangles, bool t_tetrahedra,
+		std::size_t r_edges, std::size_t r_triangles, std::size_t r_tetrahedra)
+    : toggle_edges(t_edges), toggle_triangles(t_triangles), toggle_tetrahedra(t_tetrahedra),
+      ref_edges(r_edges), ref_triangles(r_triangles), ref_tetrahedra(r_tetrahedra) {}
+
+  Configuration() {}
+};
+
+
+template <class Hasse_cell,
+	  class Simplex_cell_map>
+void populate_mesh(Mesh_medit& output,
+		   Simplex_cell_map& sc_map,
+		   Configuration configuration,
+		   std::size_t amb_d,
+		   std::map<Hasse_cell*, std::size_t> vi_map) {
+  using Mesh_element_vertices = Mesh_medit::Mesh_elements::value_type::first_type;
+  std::map<Hasse_cell*, std::size_t> ci_map;
+  std::size_t index = vi_map.size() + 1; // current size of output.vertex_points
+  if (sc_map.size() >= 3) 
+    for (const auto& sc_pair: sc_map[2]) {
+      Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+      std::set<std::size_t> vertex_indices;
+      Hasse_cell* cell = sc_pair.second;
+      for (const auto& ei_pair: cell->get_boundary())
+	for (const auto& vi_pair: ei_pair.first->get_boundary())
+	  vertex_indices.emplace(vi_map[vi_pair.first]);
+      for (const std::size_t& v: vertex_indices)
+	barycenter += output.vertex_points[v-1];
+      ci_map.emplace(std::make_pair(cell, index++));
+      output.vertex_points.emplace_back((1./vertex_indices.size()) * barycenter);    
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist = " (" + std::to_string(index-1) + ")";
+      for (const std::size_t& v: vertex_indices)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(cell), vlist));
+#endif
+    }
+
+  if (configuration.toggle_edges && sc_map.size() >= 2)
+    for (const auto& sc_map: sc_map[1]) {
+      Hasse_cell* edge_cell = sc_map.second;
+      Mesh_element_vertices edge;
+      for (const auto& vi_pair: edge_cell->get_boundary())
+	edge.push_back(vi_map[vi_pair.first]);
+      output.edges.emplace_back(std::make_pair(edge, configuration.ref_edges));
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist;
+      for (const std::size_t& v: edge)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(edge_cell), vlist));
+#endif
+    }
+  
+  if (configuration.toggle_triangles && sc_map.size() >= 3)
+    for (const auto& sc_pair: sc_map[2]) {
+      for (const auto& ei_pair: sc_pair.second->get_boundary()) {
+	Mesh_element_vertices triangle(1, ci_map[sc_pair.second]);
+	for (const auto& vi_pair: ei_pair.first->get_boundary())
+	  triangle.push_back(vi_map[vi_pair.first]);
+	output.triangles.emplace_back(std::make_pair(triangle, configuration.ref_triangles));
+      }
+    }
+  
+  if (configuration.toggle_tetrahedra && sc_map.size() >= 4)
+    for (const auto& sc_pair: sc_map[3]) {
+      Eigen::VectorXd barycenter = Eigen::VectorXd::Zero(amb_d);
+      std::set<std::size_t> vertex_indices;
+      Hasse_cell* cell = sc_pair.second;
+      for (const auto& ci_pair: cell->get_boundary())
+	for (const auto& ei_pair: ci_pair.first->get_boundary())
+	  for (const auto& vi_pair: ei_pair.first->get_boundary())
+	    vertex_indices.emplace(vi_map[vi_pair.first]);
+      for (const std::size_t& v: vertex_indices)
+	barycenter += output.vertex_points[v-1];
+      output.vertex_points.emplace_back((1./vertex_indices.size()) * barycenter);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+      std::string vlist = " (" + std::to_string(index) + ")";
+      for (const std::size_t& v: vertex_indices)
+	vlist += " " + std::to_string(v);
+      cell_vlist_map.emplace(std::make_pair(to_string(cell), vlist));
+#endif
+
+      for (const auto& ci_pair: cell->get_boundary())
+	for (const auto& ei_pair: ci_pair.first->get_boundary()) {
+	  Mesh_element_vertices tetrahedron = {index, ci_map[sc_pair.second]};
+	  for (const auto& vi_pair: ei_pair.first->get_boundary())
+	    tetrahedron.push_back(vi_map[vi_pair.first]);
+	  output.tetrahedra.emplace_back(std::make_pair(tetrahedron, configuration.ref_tetrahedra));
+	}
+      index++;
+    }  
+}
+
+template <class Cell_complex>
+Mesh_medit build_mesh_from_cell_complex(const Cell_complex& cell_complex,
+					Configuration i_configuration = Configuration(),
+					Configuration b_configuration = Configuration()) {
+  using Hasse_cell = typename Cell_complex::Hasse_cell;
+  Mesh_medit output;
+  std::map<Hasse_cell*, std::size_t> vi_map; // one for vertices, other for 2d-cells
+  std::size_t index = 1; // current size of output.vertex_points
+
+  if (cell_complex.cell_point_map().empty())
+    return output;
+  std::size_t amb_d = std::min((int) cell_complex.cell_point_map().begin()->second.size(), 3);
+  
+  for (const auto& cp_pair: cell_complex.cell_point_map()) {
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+    std::string vlist;
+    vlist += " " + std::to_string(index);
+    cell_vlist_map.emplace(std::make_pair(to_string(cp_pair.first), vlist));
+#endif
+    vi_map.emplace(std::make_pair(cp_pair.first, index++));
+    output.vertex_points.push_back(cp_pair.second);
+    output.vertex_points.back().conservativeResize(amb_d);
+  }
+  
+
+  populate_mesh(output, cell_complex.interior_simplex_cell_maps(), i_configuration, amb_d, vi_map);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+  for (const auto& sc_map: cell_complex.interior_simplex_cell_maps())    
+    for (const auto& sc_pair: sc_map) {
+      std::string simplex = "I" + to_string(sc_pair.first);
+      std::string cell = to_string(sc_pair.second);
+      std::string vlist = cell_vlist_map.at(cell).substr(1);
+      simplex_vlist_map.emplace(std::make_pair(simplex, vlist));
+    }
+#endif  
+  populate_mesh(output, cell_complex.boundary_simplex_cell_maps(), b_configuration, amb_d, vi_map);  
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+  for (const auto& sc_map: cell_complex.boundary_simplex_cell_maps())    
+    for (const auto& sc_pair: sc_map) {
+      std::string simplex = "B" + to_string(sc_pair.first);
+      std::string cell = to_string(sc_pair.second);
+      std::string vlist = cell_vlist_map.at(cell).substr(1);
+      simplex_vlist_map.emplace(std::make_pair(simplex, vlist));
+    }
+#endif  
+  return output;
+}
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
new file mode 100644
index 00000000..43cca967
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_debug_traces_to_html.h
@@ -0,0 +1,589 @@
+#ifdef GUDHI_DEBUG
+#define GUDHI_COX_OUTPUT_TO_HTML
+
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <list>
+#include <string>
+#include <regex>
+
+#include <Eigen/Dense>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+template <class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T>& vector) {
+  os << "(";  
+  if (vector.empty()) {
+    os << ")";
+    return os;
+  }
+  auto v_it = vector.begin();
+  os << *v_it++;
+  for (; v_it != vector.end(); ++v_it)
+    os << ", " << *v_it;
+  os << ")";
+  return os;
+}
+
+/* A class to make the vector horizontal instead of vertical */
+struct Straighten {
+  Straighten (const Eigen::VectorXd& vector) : vector_(vector) {}
+  const Eigen::VectorXd& vector_;
+};
+
+std::ostream& operator<<(std::ostream& os, const Straighten& str) {
+  std::size_t size = str.vector_.size();
+  os << "(" << str.vector_(0);
+  if (size == 0) {
+    os << ")";
+    return os;
+  }
+  for (std::size_t i = 1; i < size; ++i)
+    os << ", " << str.vector_(i);
+  os << ")";
+  return os;
+}
+
+std::string id_from_simplex(const std::string& simplex) {
+  std::regex r("\\s+"), r2("\\(|\\)|\\{|\\}"), r3(","), r4("\\["), r5("\\]");
+  std::string output = std::regex_replace(simplex, r, "");
+  output = std::regex_replace(output, r2, ":");
+  output = std::regex_replace(output, r3, ".");
+  output = std::regex_replace(output, r4, "_");
+  output = std::regex_replace(output, r5, "");
+  return output;
+}
+
+template <typename T>
+std::string to_string(const T& t) {
+  std::ostringstream oss;
+  oss << t;
+  return oss.str();
+}
+
+struct MT_inserted_info {
+  std::string qr_face_, init_face_, qr_intersection_;
+  bool qr_success_, is_boundary_;
+  template <class Query_result,
+	    class Simplex_handle>
+  MT_inserted_info(const Query_result& qr, const Simplex_handle& face, bool is_boundary)
+    : qr_face_(to_string(face)), init_face_(to_string(face)),
+      qr_intersection_(to_string(qr.intersection)),
+      qr_success_(qr.success), is_boundary_(is_boundary) {}
+};
+std::list<MT_inserted_info> mt_seed_inserted_list, mt_inserted_list;
+
+struct CC_summary_info {
+  std::string face_, cell_;
+  template <class SC_pair>
+  CC_summary_info(const SC_pair& sc_pair)
+    : face_(to_string(sc_pair.first)), cell_(to_string(sc_pair.second)) {}
+};
+using CC_summary_list = std::list<CC_summary_info>;
+std::vector<CC_summary_list> cc_interior_summary_lists, cc_boundary_summary_lists;
+
+struct CC_detail_info {
+  enum class Result_type {self, face, coface, inserted, join_single, join_is_face};
+  std::string simplex_, trigger_, init_simplex_;
+  Result_type status_;
+  bool join_trigger_ = false;
+  std::list<std::string> faces_, post_faces_, cofaces_;
+  template <class Simplex_handle>
+  CC_detail_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)) {}
+};
+using CC_detail_list = std::list<CC_detail_info>;
+std::vector<CC_detail_list> cc_interior_detail_lists, cc_boundary_detail_lists;
+std::vector<CC_detail_list> cc_interior_insert_detail_lists, cc_boundary_insert_detail_lists;
+
+struct CC_prejoin_info {
+  enum class Result_type {join_single, join_is_face, join_different, join_same};
+  std::string simplex_, join_;
+  std::vector<std::string> faces_;
+  std::size_t dimension_;
+  Result_type status_;
+  template <class Simplex_handle>
+  CC_prejoin_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)), dimension_(simplex.dimension()) {}
+};
+using CC_prejoin_list = std::list<CC_prejoin_info>;
+std::vector<CC_prejoin_list> cc_interior_prejoin_lists, cc_boundary_prejoin_lists;
+
+
+struct CC_join_info {
+  enum class Result_type {self, face, coface, inserted, join_single, join_is_face};
+  std::string simplex_, join_, trigger_;
+  Result_type status_;
+  std::list<std::string> boundary_faces_;
+  std::list<std::string> faces_, post_faces_, cofaces_;
+  template <class Simplex_handle>
+  CC_join_info(const Simplex_handle& simplex)
+    : simplex_(to_string(simplex)) {}
+};
+bool join_switch = false;
+std::vector<CC_detail_list> cc_interior_join_detail_lists, cc_boundary_join_detail_lists;
+
+std::map<std::string, std::string> cell_vlist_map;
+std::map<std::string, std::string> simplex_vlist_map;
+
+std::ostringstream mt_ostream, vis_ostream;
+std::vector<std::ostringstream> cc_summary_ostream, cc_traces_ostream;
+
+std::string simplex_format(const std::string& simplex, bool is_boundary) {
+  std::string b_simplex = (is_boundary? "B": "I") + simplex;
+  std::string tooltiptext;
+  auto it = simplex_vlist_map.find(b_simplex);
+  if (it == simplex_vlist_map.end())
+    tooltiptext = "deleted";
+  else
+    tooltiptext = simplex_vlist_map.at(b_simplex);
+  return (std::string)"<a class=\"" + (is_boundary? "boundary": "interior")
+    + "\" href=\"#" + id_from_simplex(b_simplex) + "\">" + b_simplex
+    + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+std::string simplex_format(const std::string& b_simplex) {
+  bool is_boundary = b_simplex[0] == 'B';
+  std::string tooltiptext;
+  auto it = simplex_vlist_map.find(b_simplex);
+  if (it == simplex_vlist_map.end())
+    tooltiptext = "deleted";
+  else
+    tooltiptext = simplex_vlist_map.at(b_simplex);
+  return (std::string)"<a class=\"" + (is_boundary? "boundary": "interior")
+    + "\" href=\"#" + id_from_simplex(b_simplex) + "\">" + b_simplex
+    + "<span class=\"tooltiptext\">" + tooltiptext + "</span></a>";
+}
+
+
+void write_head(std::ofstream& ofs) {
+  ofs << "  <head>\n"
+      << "    <title>Cell complex debug trace</title>\n"
+      << "    <style>\n"
+      << "      a.boundary {\n"
+      << "        position: relative;\n"
+      << "        display: inline-block;\n"
+      << "	  color: darkred;\n"
+      << "	  background-color: lightgreen\n"
+      << "      }\n"
+      << "      a.interior {\n"
+      << "        position: relative;\n"
+      << "        display: inline-block;\n"
+      << "	  color: navy;\n"
+      << "	  background-color: yellow\n"
+      << "      }\n"
+      << "      .tooltiptext {\n"
+      << "        visibility: hidden;\n"
+      << "        width: 120px;\n"
+      << "        background-color: #555;\n"
+      << "        color: #fff;\n"
+      << "        text-align: center;\n"
+      << "        padding: 5px 0;\n"
+      << "        border-radius: 6px;\n"
+      << "        position: absolute;\n"
+      << "        z-index: 1;\n"
+      << "        bottom: 125%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -60px;\n"
+      << "        opacity: 0;\n"
+      << "        transition: opacity 0.3s;\n"
+      << "      }\n"
+      << "      .boundary .tooltiptext::after {\n"
+      << "        content: \"\";\n"
+      << "        position: absolute;\n"
+      << "        top: 100%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -5px;\n"
+      << "        border-width: 5px;\n"
+      << "        border-style: solid;\n"
+      << "        border-color: #555 transparent transparent transparent;\n"
+      << "      }\n"
+      << "      .interior .tooltiptext::after {\n"
+      << "        content: \"\";\n"
+      << "        position: absolute;\n"
+      << "        top: 100%;\n"
+      << "        left: 50%;\n"
+      << "        margin-left: -5px;\n"
+      << "        border-width: 5px;\n"
+      << "        border-style: solid;\n"
+      << "        border-color: #555 transparent transparent transparent;\n"
+      << "      }\n"
+      << "      .boundary:hover .tooltiptext {\n"
+      << "        visibility: visible;\n"
+      << "        opacity: 1;\n"
+      << "      }\n"    
+      << "      .interior:hover .tooltiptext {\n"
+      << "        visibility: visible;\n"
+      << "        opacity: 1;\n"
+      << "      }\n"    
+      << "      ul.nav {\n"
+      << "	  list-style-type: none;\n"
+      << "	  margin: 0;\n"
+      << "	  padding: 0;\n"
+      << "	  overflow: auto;\n"
+      << "	  background-color: #333;\n"
+      << "	  position: fixed;\n"
+      << "	  height: 100%;\n"
+      << "	  width: 15%;\n"
+      << "      }\n"
+      << "      ul.nav li a {\n"
+      << "	  display: block;\n"
+      << "	  color: white;\n"
+      << "	  text-align: left;\n"
+      << "	  padding: 14px 16px;\n"
+      << "	  text-decoration: none;\n"
+      << "      }\n"
+      << "      .active {\n"
+      << "	  background-color: #4CAF50;\n"
+      << "      }\n"
+      << "      div {\n"
+      << "        margin-left: 15%;\n"
+      << "        padding: 1px 16px\n"
+      << "      }\n"
+      << "      div.navi {\n"
+      << "        margin-left: 0%;\n"
+      << "        padding: 0px 0px\n"
+      << "      }\n"
+      << "      h1 {\n"
+      << "        margin-left: 15%;\n"
+      << "        padding: 1px 16px\n"
+      << "      }\n"
+      << "    </style>\n"
+      << "  </head>\n";
+}
+
+void write_nav(std::ofstream& ofs) {
+  ofs << "    <div class=\"navi\" style=\"margin-top:30px;background-color:#1abc9c;\">\n"
+      << "    <ul class=\"nav\">\n"
+      << "      <li><a href=\"#mant\">Manifold tracing</a></li>\n"
+      << "      <li><a href=\"#cell\">Cell complex</a>\n"
+      << "        <ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "          <li><a href=\"#dim" << i << "\">Dimension " << i << "</a>\n";
+    ofs << "            <ul>\n";
+    ofs << "              <li><a href=\"#dim" << i << "i\">Interior</a></li>\n";
+    if (i < cc_boundary_summary_lists.size()) {
+      ofs << "              <li><a href=\"#dim" << i << "b\">Boundary</a></li>\n";
+    }
+    ofs << "            </ul>\n";
+    ofs << "          </li>\n";
+  }
+  ofs << "        </ul>\n"
+      << "      </li>\n"
+      << "      <li><a href=\"#visu\">Visualization details</a></li>\n"
+      << "    </ul>\n"
+      << "    </div>\n";
+}
+
+void write_mt(std::ofstream& ofs) {
+  ofs << "    <div id=\"mant\">\n";
+  ofs << "      <h2> Manifold debug trace </h2>\n";
+  ofs << "      <h3> Simplices inserted during the seed phase </h3>\n";
+  ofs << "      <ul>\n";
+  for (const MT_inserted_info& mt_info: mt_seed_inserted_list) {
+    if (mt_info.qr_success_) {
+      ofs << "        <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+      if (mt_info.qr_face_ != mt_info.init_face_)
+	ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+      ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+    }
+    else
+      ofs << "        <li>Failed to insert " << mt_info.init_face_
+	  << "</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "      <h3> Simplices inserted during the while loop phase </h3>\n";
+  ofs << "      <ul>\n";
+  for (const MT_inserted_info& mt_info: mt_inserted_list) {
+    if (mt_info.qr_success_) {
+      ofs << "        <li>Inserted " << simplex_format(mt_info.qr_face_, mt_info.is_boundary_);
+      if (mt_info.qr_face_ != mt_info.init_face_)
+	ofs << " (initially " << simplex_format(mt_info.init_face_, mt_info.is_boundary_) << ")";
+      ofs << " intersection point is " << mt_info.qr_intersection_ << "</li>\n";
+    }
+    else
+      ofs << "        <li>Failed to insert " << mt_info.init_face_
+	  << ")</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "    </div>\n";
+}
+
+void write_cc(std::ofstream& ofs) {
+  ofs << "    <div id=\"cell\">\n"
+      << "      <h2> Cell complex debug trace </h2>\n"
+      << "      <p>Go to:</p>\n"
+      << "      <ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "	<li><a href=\"#dim" << i << "\">Dimension " << i << "</a></li>\n";
+  }
+  ofs << "      </ul>\n";
+  for (std::size_t i = 0; i < cc_interior_summary_lists.size(); ++i) {
+    ofs << "      <h3 id=\"dim" << i << "\"> Dimension " << i << "</h3>\n";
+    ofs << "      <h4 id=\"dim" << i << "i\"> Summary for interior simplices</h4>\n";
+    if (i < cc_boundary_summary_lists.size())
+      ofs << "      <p><a href=\"#dim" << i << "b\">Go to boundary</a></p>\n";
+    ofs << "        <ul>\n";
+    for (const CC_summary_info& cc_info: cc_interior_summary_lists[i])
+      ofs << "          <li id = \"" << id_from_simplex("I" + cc_info.face_) << "\">"
+	  << simplex_format(cc_info.face_, false)
+	  << " cell =" << cc_info.cell_ << "</li>\n";
+    ofs << "        </ul>\n";
+    ofs << "      <h4> Prejoin state of the interior cells of dimension " << i << "</h4>\n";
+    auto prejoin_it = cc_interior_prejoin_lists[i].begin();
+    while (prejoin_it != cc_interior_prejoin_lists[i].end()) {
+      std::size_t j = prejoin_it->dimension_;
+      ofs << "        <h5>" << j << "-dimensional ambient simplices</h5>\n";
+      ofs << "          <ul>\n";
+      for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+	ofs << "          <li>" << simplex_format(prejoin_it->simplex_, false)
+	    << " join = " << simplex_format(prejoin_it->join_, false)
+	    << " boundary:\n"
+	    << "            <ul>\n";
+	for (const auto& face: prejoin_it->faces_)
+	  ofs << "              <li>" << simplex_format(face) << "</li>";
+	ofs << "            </ul>\n";
+	switch (prejoin_it->status_) {
+	case (CC_prejoin_info::Result_type::join_single):
+	  ofs << "            <p style=\"color: red\">Deleted "
+	      << simplex_format(prejoin_it->simplex_, false)
+	      << " as it has a single face.</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_is_face):
+	  ofs << "            <p style=\"color: red\">Deleted "
+	      << simplex_format(prejoin_it->simplex_, false)
+	      << " as its join " << simplex_format(prejoin_it->join_, false)
+	      << " is one of the faces.</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_different):
+	  ofs << "            <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, false)
+	      << " and replaced by its join " << simplex_format(prejoin_it->join_, false)
+	      << ".</p>";
+	  break;
+	case (CC_prejoin_info::Result_type::join_same):
+	  ofs << "            <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, false)
+	      << ".</p>";
+	}
+	ofs << "          </li>";
+      }
+      ofs << "          </ul>\n";
+    }
+    ofs << "      <h4> Details for interior simplices</h4>\n";
+    ofs << "        <ul>\n";
+    for (const CC_detail_info& cc_info: cc_interior_detail_lists[i]) {
+      if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	ofs << "          <li style=\"color:magenta\" id = \""
+	    << id_from_simplex("I" + cc_info.simplex_) << "\"> Simplex "
+	    << simplex_format(cc_info.simplex_, false) << " has only one face ("
+	    << simplex_format(cc_info.trigger_, false) << ") and is deleted.";
+	continue;
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	ofs << "          <li style=\"color:darkmagenta\" id = \""
+	    << id_from_simplex("I" + cc_info.simplex_) << "\"> The join of the simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is one of its faces ("
+	    << simplex_format(cc_info.trigger_, false) << "), hence it is is deleted.";
+	continue;
+      }
+      ofs << "          <li> Insert_cell called for " << simplex_format(cc_info.simplex_, false)
+	  << "\n";
+      ofs << "            <ul>\n";
+      // for (const std::string& cof: cc_info.faces_)
+      // 	ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, false)
+      // 	    << " is a face of " << simplex_format(cof, false) << "\n";
+      ofs << "            </ul>\n";
+      ofs << "            <ul>\n";
+      if (cc_info.status_ == CC_detail_info::Result_type::self) {
+	ofs << "            <p><span style=\"color:blue\">The simplex "
+	    << simplex_format(cc_info.simplex_, false)
+	    << " already exists in the cell complex!</span></p>\n";
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::face) {
+	ofs << "            <p><span style=\"color:red\">The simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is a face of the simplex "
+	    << simplex_format(cc_info.trigger_, false) << "!</span><br>\n";
+	ofs << "              <ul>\n";
+	for (const std::string post_face: cc_info.post_faces_)
+	  ofs << "                <li id = \"" << id_from_simplex("I" + post_face) << "\">"
+	      << "Post deleting " << simplex_format(post_face, false) << "</li>\n";
+	ofs << "              </ul>\n";
+	ofs << "            </p>\n";
+	ofs << "            <p id = \"" << id_from_simplex("I" + cc_info.trigger_) << "\">"
+	    << "Deleting " << simplex_format(cc_info.trigger_, false) << "</p>\n";
+      }
+      // for (const std::string& fac: cc_info.cofaces_)
+      // 	ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, false)
+      // 	    << " is a coface of " << simplex_format(fac, false) << "\n";
+      if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+	ofs << "            <p><span style=\"color:darkorange\">The simplex "
+	    << simplex_format(cc_info.simplex_, false) << " is a coface of the simplex "
+	    << simplex_format(cc_info.trigger_, false) << "!</span><p>\n";
+      }
+      if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+	ofs << "            <p><span style=\"color:green\">Successfully inserted "
+	    << simplex_format(cc_info.simplex_, false) << "!</span><p>\n";
+      }
+      ofs << "            </ul>\n";
+      ofs << "          </li>\n";
+    }
+    ofs << "        </ul>\n";
+
+    if (i < cc_boundary_summary_lists.size()) {
+      ofs << "      <h4 id=\"dim" << i << "b\"> Summary for boundary simplices</h4>\n";
+      ofs << "      <p><a href=\"#dim" << i << "i\">Go to interior</a></p>\n";
+      ofs << "        <ul>\n";
+      for (const CC_summary_info& cc_info: cc_boundary_summary_lists[i])
+	ofs << "          <li  id = \"" << id_from_simplex("B" + cc_info.face_) << "\">"
+	    << simplex_format(cc_info.face_, true)
+	    << " cell =" << cc_info.cell_ << "</li>\n";
+      ofs << "        </ul>\n";
+      ofs << "      <h4> Prejoin state of the boundary cells of dimension " << i << "</h4>\n";
+      auto prejoin_it = cc_boundary_prejoin_lists[i].begin();
+      while (prejoin_it != cc_boundary_prejoin_lists[i].end()) {
+	std::size_t j = prejoin_it->dimension_;
+	ofs << "        <h5>" << j << "-dimensional ambient simplices</h5>\n";
+	ofs << "          <ul>\n";
+	for (; prejoin_it->dimension_ == j; ++prejoin_it) {
+	  ofs << "          <li>" << simplex_format(prejoin_it->simplex_, true)
+	      << " join = " << simplex_format(prejoin_it->join_, true)
+	      << " boundary:\n"
+	      << "            <ul>\n";
+	  for (const auto& face: prejoin_it->faces_)
+	    ofs << "              <li>" << simplex_format(face) << "</li>";
+	  ofs << "            </ul>\n";
+	  switch (prejoin_it->status_) {
+	  case (CC_prejoin_info::Result_type::join_single):
+	    ofs << "            <p style=\"color: red\">Deleted "
+		<< simplex_format(prejoin_it->simplex_, true)
+		<< " as it has a single face.</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_is_face):
+	    ofs << "            <p style=\"color: red\">Deleted "
+		<< simplex_format(prejoin_it->simplex_, true)
+		<< " as its join " << simplex_format(prejoin_it->join_, true)
+		<< " is one of the faces.</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_different):
+	    ofs << "            <p style=\"color: magenta\">Deleted " << simplex_format(prejoin_it->simplex_, true)
+		<< " and replaced by its join " << simplex_format(prejoin_it->join_, true)
+		<< ".</p>";
+	    break;
+	  case (CC_prejoin_info::Result_type::join_same):
+	    ofs << "            <p style=\"color: green\">Kept " << simplex_format(prejoin_it->simplex_, true)
+		<< ".</p>";
+	  }
+	  ofs << "          </li>";
+	}
+	ofs << "          </ul>\n";
+      }
+    }
+    if (i < cc_boundary_detail_lists.size()) {
+      ofs << "      <h4> Details for boundary simplices</h4>\n"
+	  << "        <ul>\n";
+      for (const CC_detail_info& cc_info: cc_boundary_detail_lists[i]) {
+	if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	  ofs << "          <li style=\"color:magenta\" id = \""
+	      << id_from_simplex("B" + cc_info.simplex_) << "\"> Simplex "
+	      << simplex_format(cc_info.simplex_, true) << " has only one face ("
+	      << simplex_format(cc_info.trigger_, true) << ") and is deleted.";
+	  continue;
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::join_single) {
+	  ofs << "          <li style=\"color:darkmagenta\" id = \""
+	      << id_from_simplex("B" + cc_info.simplex_) << "\"> The join of the simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is one of its faces ("
+	      << simplex_format(cc_info.trigger_, true) << "), hence it is is deleted.";
+	  continue;
+	}
+	ofs << "          <li> Insert_simplex called on " << simplex_format(cc_info.simplex_, true);
+	ofs << "            <ul>\n";
+	// for (const std::string& cof: cc_info.faces_)
+	//   ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, true)
+	//       << " is a face of " << simplex_format(cof, true) << "\n";
+	ofs << "            </ul>\n";
+	ofs << "            <ul>\n";
+	if (cc_info.status_ == CC_detail_info::Result_type::self) {
+	  ofs << "            <p><span style=\"color:blue\">The simplex "
+	      << simplex_format(cc_info.simplex_, true)
+	      << " already exists in the cell complex!</span></p>\n";
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::face) {
+	  ofs << "            <p><span style=\"color:red\">The simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is a face of the simplex "
+	      << simplex_format(cc_info.trigger_, true) << "!</span><br>\n";
+	  ofs << "              <ul>\n";
+	  for (const std::string post_face: cc_info.post_faces_)
+	    ofs << "                <li id=\"" << id_from_simplex("B" + post_face)
+		<< "\">Post deleting " << simplex_format(post_face, true)
+		<< "</li>\n";
+	  ofs << "              </ul>\n";
+	  ofs << "            </p>\n";
+	  ofs << "            <p id=\"" << id_from_simplex(cc_info.trigger_)
+	      << "\">Deleting " << simplex_format(cc_info.trigger_, true) << "</p>\n";
+	}
+	// for (const std::string& fac: cc_info.cofaces_)
+	//   ofs << "              <li>Checking if " << simplex_format(cc_info.simplex_, true)
+	//       << " is a coface of " << simplex_format(fac, true) << "\n";
+	ofs << "            </ul>\n";
+	ofs << "          </li>\n";
+	if (cc_info.status_ == CC_detail_info::Result_type::coface) {
+	  ofs << "            <p><span style=\"color:darkorange\">The simplex "
+	      << simplex_format(cc_info.simplex_, true) << " is a coface of the simplex "
+	      << simplex_format(cc_info.trigger_, true) << "!</span><p>\n";
+	}
+	if (cc_info.status_ == CC_detail_info::Result_type::inserted) {
+	  ofs << "            <p><span style=\"color:green\">Successfully inserted "
+	      << simplex_format(cc_info.simplex_, true) << "!</span><p>\n";
+	}
+      }
+      ofs << "        </ul>\n";
+    }
+  }  
+  ofs << "    </div>\n";
+}
+
+void write_visu(std::ofstream& ofs) {
+  ofs << "    <div id=\"visu\">\n"
+      << "      <h2> Visualization details debug trace </h2>\n";
+  // std::vector<std::map<std::string, std::string> > vs_maps(cc_interior_summary_lists.size());
+  std::map<std::string, std::string> vs_map;
+  for (const auto& sv_pair: simplex_vlist_map)
+    vs_map.emplace(std::make_pair(sv_pair.second, sv_pair.first));
+  ofs << "      <ul>\n";
+  for (const auto& vs_pair: vs_map) {
+    std::string w_simplex = vs_pair.second.substr(1);
+    bool is_boundary = vs_pair.second[0] == 'B';
+    ofs << "      <li><b>" << vs_pair.first << "</b>: "
+	<< simplex_format(w_simplex, is_boundary) << "</li>\n";
+  }
+  ofs << "      </ul>\n";
+  ofs << "    </div>\n";
+}
+
+void write_to_html(std::string file_name) {
+  std::ofstream ofs(file_name + ".html", std::ofstream::out);
+  ofs << "<!DOCTYPE html>\n"
+      << "<html>\n";
+  write_head(ofs);
+  ofs << "  <body>\n";
+  write_nav(ofs);
+  ofs << "    <h1> Debug traces for " << file_name << " </h1>\n";
+  write_mt(ofs);
+  write_cc(ofs);
+  write_visu(ofs);  
+  ofs << "  </body>\n";
+  ofs << "</html>\n";
+
+  ofs.close();
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
+
diff --git a/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
new file mode 100644
index 00000000..fe379242
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/IO/output_meshes_to_medit.h
@@ -0,0 +1,181 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IO_OUTPUT_MESHES_TO_MEDIT_H_
+#define IO_OUTPUT_MESHES_TO_MEDIT_H_
+
+#include <gudhi/IO/Mesh_medit.h>
+#include <Eigen/Dense>
+#include <fstream>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+using Vertex_points = Mesh_medit::Vertex_points;
+using Mesh_elements = Mesh_medit::Mesh_elements;
+using Scalar_field_range = Mesh_medit::Scalar_field_range;
+
+template <std::size_t I = 0,
+	  typename... Meshes>
+typename std::enable_if<I == sizeof... (Meshes), void>::type
+fill_meshes (Vertex_points& vertex_points,
+	     Mesh_elements& edges,
+	     Mesh_elements& triangles,
+	     Mesh_elements& tetrahedra,
+	     Scalar_field_range& triangles_scalar_range,
+	     Scalar_field_range& tetrahedra_scalar_range,
+	     std::size_t index,
+	     const Meshes&... meshes) {
+}
+
+template <std::size_t I = 0,
+	  typename... Meshes>
+typename std::enable_if<I != sizeof... (Meshes), void>::type
+fill_meshes (Vertex_points& vertex_points,
+	     Mesh_elements& edges,
+	     Mesh_elements& triangles,
+	     Mesh_elements& tetrahedra,
+	     Scalar_field_range& triangles_scalar_range,
+	     Scalar_field_range& tetrahedra_scalar_range,
+	     std::size_t index,
+	     const Meshes&... meshes) {
+  auto mesh = std::get<I>(std::forward_as_tuple(meshes...));
+  for (const auto& v: mesh.vertex_points)
+    vertex_points.push_back(v);
+  for (const auto& e: mesh.edges) {
+    std::vector<std::size_t> edge;
+    for (const auto& v_i: e.first)
+      edge.push_back(v_i + index);
+    edges.emplace_back(std::make_pair(edge, e.second));
+  }
+  for (const auto& t: mesh.triangles) {
+    std::vector<std::size_t> triangle;
+    for (const auto& v_i: t.first)
+      triangle.push_back(v_i + index);
+    triangles.emplace_back(std::make_pair(triangle, t.second));
+  }
+  for (const auto& t: mesh.tetrahedra) {
+    std::vector<std::size_t> tetrahedron;
+    for (const auto& v_i: t.first)
+      tetrahedron.push_back(v_i + index);
+    tetrahedra.emplace_back(std::make_pair(tetrahedron, t.second));
+  }
+  for (const auto& b: mesh.triangles_scalar_range)
+    triangles_scalar_range.push_back(b);
+  for (const auto& b: mesh.tetrahedra_scalar_range)
+    tetrahedra_scalar_range.push_back(b);
+  fill_meshes<I+1, Meshes...>(vertex_points,
+			      edges,
+			      triangles,
+			      tetrahedra,
+			      triangles_scalar_range,
+			      tetrahedra_scalar_range,
+			      index + mesh.vertex_points.size(),
+			      meshes...);
+}
+
+/** \brief Outputs a text file with specified meshes that can be visualized in Medit.
+ *  
+ *  @param[in] amb_d Ambient dimension. Can be 2 or 3.
+ *  @param[in] file_name The name of the output file.
+ *  @param[in] meshes A pack of meshes to be specified separated by commas.
+ */
+template <typename... Meshes>
+void output_meshes_to_medit(std::size_t amb_d, std::string file_name, const Meshes&... meshes) {
+  Vertex_points vertex_points;
+  Mesh_elements edges, triangles, tetrahedra;
+  Scalar_field_range triangles_scalar_range, tetrahedra_scalar_range;
+  fill_meshes(vertex_points,
+	      edges,
+	      triangles,
+	      tetrahedra,
+	      triangles_scalar_range,
+	      tetrahedra_scalar_range,
+	      0,
+	      meshes...);
+  
+  std::ofstream ofs (file_name + ".mesh", std::ofstream::out);
+  std::ofstream ofs_bb (file_name + ".bb", std::ofstream::out);
+  
+  if (amb_d == 2) {
+    ofs << "MeshVersionFormatted 1\nDimension 2\n";
+    ofs_bb << "2 1 ";
+    ofs << "Vertices\n" << vertex_points.size() << "\n";
+    for (auto p: vertex_points) {
+      ofs << p[0] << " " << p[1] << " 2\n";
+    }
+    ofs << "Edges " << edges.size() << "\n";
+    for (auto e: edges) {
+      for (auto v: e.first)      
+    	ofs << v << " ";
+      ofs << e.second << std::endl;
+    }
+    ofs << "Triangles " << triangles.size() << "\n";
+    for (auto s: triangles) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+
+    ofs_bb << triangles_scalar_range.size() << " 1\n";
+    for (auto& b: triangles_scalar_range)
+      ofs_bb << b << "\n";
+
+  }
+  else {
+    ofs << "MeshVersionFormatted 1\nDimension 3\n";
+    ofs_bb << "3 1 ";
+    ofs << "Vertices\n" << vertex_points.size() << "\n";
+    for (auto p: vertex_points) {
+      ofs << p[0] << " " << p[1] << " " << p[2] << " 2\n";
+    }
+    ofs << "Edges " << edges.size() << "\n";
+    for (auto e: edges) {
+      for (auto v: e.first)      
+    	ofs << v << " ";
+      ofs << e.second << std::endl;
+    }
+    ofs << "Triangles " << triangles.size() << "\n";
+    for (auto s: triangles) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+    ofs << "Tetrahedra " << tetrahedra.size() << "\n";
+    for (auto s: tetrahedra) {
+      for (auto v: s.first) {
+    	ofs << v << " ";
+      }
+      ofs << s.second << std::endl;
+    }
+
+    ofs_bb << triangles_scalar_range.size() + tetrahedra_scalar_range.size() << " 1\n";
+    for (auto& b: triangles_scalar_range)
+      ofs_bb << b << "\n";
+    for (auto& b: tetrahedra_scalar_range)
+      ofs_bb << b << "\n";
+  }
+
+  
+  ofs.close();
+  ofs_bb.close();
+
+}
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
new file mode 100644
index 00000000..a522f691
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Implicit_manifold_intersection_oracle.h
@@ -0,0 +1,279 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
+#define IMPLICIT_MANIFOLD_INTERSECTION_ORACLE_H_
+
+#include <Eigen/Dense>
+
+#include <gudhi/Permutahedral_representation/face_from_indices.h>
+#include <gudhi/Functions/Constant_function.h>
+#include <gudhi/Functions/PL_approximation.h>
+#include <gudhi/Query_result.h>
+
+#include <vector>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+
+
+/** \class Implicit_manifold_intersection_oracle
+ *  \brief An oracle that supports the intersection query on an implicit manifold.
+ *
+ *  \tparam Function_ The function template parameter. Should be a model of 
+ *   the concept FunctionForImplicitManifold.
+ *  \tparam Domain_function_ The domain function template parameter. Should be a model of
+ *   the concept FunctionForImplicitManifold.
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template<class Function_,
+	 class Domain_function_ = Constant_function>
+class Implicit_manifold_intersection_oracle {
+
+  /* Computes the affine coordinates of the intersection point of the implicit manifold
+   * and the affine hull of the simplex. */
+  template <class Simplex_handle, 
+	    class Triangulation>
+  Eigen::VectorXd compute_lambda(const Simplex_handle& simplex,
+				 const Triangulation& triangulation) const {
+    std::size_t cod_d = this->cod_d();
+    Eigen::MatrixXd matrix(cod_d + 1, cod_d + 1);
+    for (std::size_t i = 0; i < cod_d + 1; ++i)
+      matrix(0, i) = 1;
+    std::size_t j = 0;
+    for (auto v: simplex.vertex_range()) {
+      Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
+      for (std::size_t i = 1; i < cod_d + 1; ++i)
+	matrix(i, j) = v_coords(i-1);
+      j++;
+    }
+    Eigen::VectorXd z(cod_d + 1);
+    z(0) = 1;
+    for (std::size_t i = 1; i < cod_d + 1; ++i)
+      z(i) = 0;
+    Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+    return lambda;
+  }
+
+  /* Computes the affine coordinates of the intersection point of the boundary
+   * of the implicit manifold and the affine hull of the simplex. */
+  template <class Simplex_handle, 
+	    class Triangulation>
+  Eigen::VectorXd compute_boundary_lambda(const Simplex_handle& simplex,
+					  const Triangulation& triangulation) const {
+    std::size_t cod_d = this->cod_d();
+    Eigen::MatrixXd matrix(cod_d + 2, cod_d + 2);
+    for (std::size_t i = 0; i < cod_d + 2; ++i)
+      matrix(0, i) = 1;
+    std::size_t j = 0;
+    for (auto v: simplex.vertex_range()) {
+      Eigen::VectorXd v_coords = fun_(triangulation.cartesian_coordinates(v));
+      for (std::size_t i = 1; i < cod_d + 1; ++i)
+	matrix(i, j) = v_coords(i-1);
+      Eigen::VectorXd bv_coords = domain_fun_(triangulation.cartesian_coordinates(v));
+      matrix(cod_d + 1, j) = bv_coords(0);
+      j++;
+    }
+    Eigen::VectorXd z(cod_d + 2);
+    z(0) = 1;
+    for (std::size_t i = 1; i < cod_d + 2; ++i)
+      z(i) = 0;
+    Eigen::VectorXd lambda = matrix.colPivHouseholderQr().solve(z);
+    return lambda;
+  }
+
+  /* Computes the intersection result for a given simplex in a triangulation. */
+  template <class Simplex_handle,
+	    class Triangulation>
+  Query_result<Simplex_handle> intersection_result(const Eigen::VectorXd& lambda,
+						   const Simplex_handle& simplex,
+						   const Triangulation& triangulation) const {
+    using QR = Query_result<Simplex_handle>;
+    std::size_t amb_d = triangulation.dimension();
+    std::size_t cod_d = simplex.dimension();
+
+    for (std::size_t i = 0; i < (std::size_t)lambda.size(); ++i)
+      if (lambda(i) < 0 || lambda(i) > 1)
+	return QR({Eigen::VectorXd(), false});
+
+    Eigen::MatrixXd vertex_matrix(cod_d + 1, amb_d);
+    auto v_range = simplex.vertex_range();
+    auto v_it = v_range.begin();
+    for (std::size_t i = 0; i < cod_d + 1 && v_it != v_range.end(); ++v_it, ++i) {
+      Eigen::VectorXd v_coords = triangulation.cartesian_coordinates(*v_it);
+      for (std::size_t j = 0; j < amb_d; ++j)
+	vertex_matrix(i, j) = v_coords(j);
+    }
+    Eigen::VectorXd intersection = lambda.transpose()*vertex_matrix;
+    return QR({intersection, true});
+  }
+  
+public:
+
+  /** \brief Ambient dimension of the implicit manifold. */
+  std::size_t amb_d() const {
+    return fun_.amb_d();
+  }
+  
+  /** \brief Codimension of the implicit manifold. */
+  std::size_t cod_d() const {
+    return fun_.cod_d();
+  }
+
+  /** \brief Intersection query with the relative interior of the manifold.
+   *  
+   *  \details The returned structure Query_result contains the boolean value
+   *   that is true only if the intersection point of the query simplex and
+   *   the relative interior of the manifold exists, the intersection point
+   *   and the face of the query simplex that contains 
+   *   the intersection point.
+   *   
+   *  \tparam Simplex_handle The class of the query simplex.
+   *   Needs to be a model of the concept SimplexInCoxeterTriangulation.
+   *  \tparam Triangulation The class of the triangulation.
+   *   Needs to be a model of the concept TriangulationForManifoldTracing.
+   *
+   *  @param[in] simplex The query simplex. The dimension of the simplex
+   *   should be the same as the codimension of the manifold 
+   *   (the codomain dimension of the function).
+   *  @param[in] triangulation The ambient triangulation. The dimension of 
+   *   the triangulation should be the same as the ambient dimension of the manifold 
+   *   (the domain dimension of the function).
+   */
+  template <class Simplex_handle,
+	    class Triangulation>
+  Query_result<Simplex_handle> intersects(const Simplex_handle& simplex,
+					  const Triangulation& triangulation) const {
+    Eigen::VectorXd lambda = compute_lambda(simplex, triangulation);
+    return intersection_result(lambda, simplex, triangulation);
+  }
+
+  /** \brief Intersection query with the boundary of the manifold.
+   *  
+   *  \details The returned structure Query_result contains the boolean value
+   *   that is true only if the intersection point of the query simplex and
+   *   the boundary of the manifold exists, the intersection point
+   *   and the face of the query simplex that contains 
+   *   the intersection point.
+   *   
+   *  \tparam Simplex_handle The class of the query simplex.
+   *   Needs to be a model of the concept SimplexInCoxeterTriangulation.
+   *  \tparam Triangulation The class of the triangulation.
+   *   Needs to be a model of the concept TriangulationForManifoldTracing.
+   *
+   *  @param[in] simplex The query simplex. The dimension of the simplex
+   *   should be the same as the codimension of the boundary of the manifold 
+   *   (the codomain dimension of the function + 1).
+   *  @param[in] triangulation The ambient triangulation. The dimension of 
+   *   the triangulation should be the same as the ambient dimension of the manifold 
+   *   (the domain dimension of the function).
+   */
+  template <class Simplex_handle,
+	    class Triangulation>
+  Query_result<Simplex_handle> intersects_boundary(const Simplex_handle& simplex,
+						   const Triangulation& triangulation) const {
+    Eigen::VectorXd lambda = compute_boundary_lambda(simplex, triangulation);
+    return intersection_result(lambda, simplex, triangulation);
+  }
+
+  
+  /** \brief Returns true if the input point lies inside the piecewise-linear
+   *   domain induced by the given ambient triangulation that defines the relative
+   *   interior of the piecewise-linear approximation of the manifold.
+   *
+   * @param p The input point. Needs to have the same dimension as the ambient
+   *  dimension of the manifold (the domain dimension of the function).
+   * @param triangulation The ambient triangulation. Needs to have the same
+   *  dimension as the ambient dimension of the manifold 
+   *  (the domain dimension of the function).
+   */
+  template <class Triangulation>
+  bool lies_in_domain(const Eigen::VectorXd& p,
+		      const Triangulation& triangulation) const {
+    Eigen::VectorXd pl_p = make_pl_approximation(domain_fun_, triangulation)(p);
+    return pl_p(0) < 0;
+  }
+
+  /** \brief Returns the function that defines the interior of the manifold */
+  const Function_& function() const {
+    return fun_;
+  }
+  
+  /** \brief Constructs an intersection oracle for an implicit manifold potentially 
+   *   with boundary from given function and domain.
+   *
+   *  @param function The input function that represents the implicit manifold
+   *   before the restriction with the domain.
+   *  @param domain_function The input domain function that can be used to define an implicit
+   *   manifold with boundary.
+   */
+  Implicit_manifold_intersection_oracle(const Function_& function,
+					const Domain_function_& domain_function)
+    : fun_(function), domain_fun_(domain_function) {}
+
+  /** \brief Constructs an intersection oracle for an implicit manifold 
+   *   without boundary from a given function.
+   *
+   *   \details To use this constructor, the template Domain_function_ needs to be left 
+   *   at its default value (Gudhi::coxeter_triangulation::Constant_function).
+   *
+   *  @param function The input function that represents the implicit manifold
+   *   without boundary.
+   */
+  Implicit_manifold_intersection_oracle(const Function_& function)
+    : fun_(function),
+      domain_fun_(function.amb_d(), 1, Eigen::VectorXd::Constant(1,-1)) {}
+  
+private:
+  Function_ fun_;
+  Domain_function_ domain_fun_;
+};
+
+/** \brief Static constructor of an intersection oracle from a function with a domain.
+ *
+ *  @param function The input function that represents the implicit manifold
+ *   before the restriction with the domain.
+ *  @param domain_function The input domain function that can be used to define an implicit
+ *   manifold with boundary.
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template<class Function_,
+	 class Domain_function_>
+Implicit_manifold_intersection_oracle<Function_, Domain_function_>
+make_oracle(const Function_& function,
+	    const Domain_function_& domain_function){
+  return Implicit_manifold_intersection_oracle<Function_, Domain_function_>(function,
+									    domain_function);
+}
+
+
+/** \brief Static constructor of an intersection oracle from a function without a domain.
+ *
+ *  @param function The input function that represents the implicit manifold
+ *   without boundary.
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template<class Function_>
+Implicit_manifold_intersection_oracle<Function_> make_oracle(const Function_& function){
+  return Implicit_manifold_intersection_oracle<Function_>(function);
+}
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h
new file mode 100644
index 00000000..fdb4f630
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Manifold_tracing.h
@@ -0,0 +1,307 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef MANIFOLD_TRACING_H_
+#define MANIFOLD_TRACING_H_
+
+#include <gudhi/Query_result.h>
+
+#include <boost/functional/hash.hpp>
+
+#include <Eigen/Dense>
+
+#include <queue>
+#include <unordered_map>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/**
+ *  \ingroup coxeter_triangulation
+ */
+
+/** \class Manifold_tracing 
+ *  \brief A class that assembles methods for manifold tracing algorithm.
+ *
+ *  \tparam Triangulation_ The type of the ambient triangulation.
+ *   Needs to be a model of the concept TriangulationForManifoldTracing.
+ */
+template <class Triangulation_>
+class Manifold_tracing {
+  
+  typedef typename Triangulation_::Simplex_handle Simplex_handle;
+  
+  struct Simplex_hash {
+    typedef Simplex_handle argument_type;
+    typedef std::size_t result_type;
+    result_type operator()(const argument_type& s) const noexcept {
+      return boost::hash<typename Simplex_handle::Vertex>()(s.vertex());
+    }
+  };
+  
+public:
+  /** \brief Type of the output simplex map with keys of type Triangulation_::Simplex_handle
+   *   and values of type Eigen::VectorXd.
+   *   This type should be used for the output in the method manifold_tracing_algorithm.
+   */
+  typedef std::unordered_map<Simplex_handle, Eigen::VectorXd, Simplex_hash> Out_simplex_map;
+
+  /**
+   * \brief Computes the set of k-simplices that intersect
+   * a boundaryless implicit manifold given by an intersection oracle, where k
+   * is the codimension of the manifold.
+   * The computation is based on the seed propagation --- it starts at the 
+   * given seed points and then propagates along the manifold.
+   *
+   * \tparam Point_range Range of points of type Eigen::VectorXd.
+   * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+   *  Needs to be a model of the concept IntersectionOracle.
+   *
+   * \param[in] seed_points The range of points on the manifold from which 
+   * the computation begins.
+   * \param[in] triangulation The ambient triangulation.
+   * \param[in] oracle The intersection oracle for the manifold.
+   * The ambient dimension needs to match the dimension of the
+   * triangulation.
+   * \param[out] out_simplex_map The output map, where the keys are k-simplices in
+   * the input triangulation that intersect the input manifold and the mapped values 
+   * are the intersection points.
+   */
+  template <class Point_range,
+	    class Intersection_oracle>
+  void manifold_tracing_algorithm(const Point_range& seed_points,
+				  const Triangulation_& triangulation,
+				  const Intersection_oracle& oracle,
+				  Out_simplex_map& out_simplex_map) {
+    std::size_t cod_d = oracle.cod_d();
+    std::queue<Simplex_handle> queue;
+
+    for (const auto& p: seed_points) {
+      Simplex_handle full_simplex = triangulation.locate_point(p);
+      for (Simplex_handle face: full_simplex.face_range(cod_d)) {
+	Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
+	if (qr.success &&
+	    out_simplex_map.emplace(std::make_pair(face, qr.intersection)).second) {
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+	  mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif
+	  queue.emplace(face);
+	  break;
+	}
+      }
+    }
+    
+
+    while (!queue.empty()) {
+      Simplex_handle s = queue.front();
+      queue.pop();
+      for (auto cof: s.coface_range(cod_d+1)) {
+	for (auto face: cof.face_range(cod_d)) {
+	  Query_result<Simplex_handle> qr = oracle.intersects(face, triangulation);
+	  if (qr.success &&
+	      out_simplex_map.emplace(std::make_pair(face, qr.intersection)).second)
+	    queue.emplace(face);
+	}
+      }
+    }
+  }
+
+  /**
+   * \brief Computes the set of k-simplices that intersect
+   * the dimensional manifold given by an intersection oracle, where k
+   * is the codimension of the manifold.
+   * The computation is based on the seed propagation --- it starts at the 
+   * given seed points and then propagates along the manifold.
+   *
+   * \tparam Point_range Range of points of type Eigen::VectorXd.
+   * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+   *  Needs to be a model of the concept IntersectionOracle.
+   *
+   * \param[in] seed_points The range of points on the manifold from which 
+   * the computation begins.
+   * \param[in] triangulation The ambient triangulation.
+   * \param[in] oracle The intersection oracle for the manifold.
+   * The ambient dimension needs to match the dimension of the
+   * triangulation.
+   * \param[out] interior_simplex_map The output map, where the keys are k-simplices in
+   * the input triangulation that intersect the relative interior of the input manifold  
+   * and the mapped values are the intersection points.
+   * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in
+   * the input triangulation that intersect the boundary of the input manifold  
+   * and the mapped values are the intersection points.
+   */
+  template <class Point_range,
+	    class Intersection_oracle>
+  void manifold_tracing_algorithm(const Point_range& seed_points,
+				  const Triangulation_& triangulation,
+				  const Intersection_oracle& oracle,
+				  Out_simplex_map& interior_simplex_map,
+				  Out_simplex_map& boundary_simplex_map) {
+    std::size_t cod_d = oracle.cod_d();
+    std::queue<Simplex_handle> queue;
+
+    for (const auto& p: seed_points) {
+      Simplex_handle full_simplex = triangulation.locate_point(p);
+      for (Simplex_handle face: full_simplex.face_range(cod_d)) {
+	auto qr = oracle.intersects(face, triangulation);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+	mt_seed_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif	
+	if (qr.success) {
+	  if (oracle.lies_in_domain(qr.intersection, triangulation)) {
+	    if (interior_simplex_map.emplace(std::make_pair(face, qr.intersection)).second)
+	      queue.emplace(face);
+	  }
+	  else {
+	    for (Simplex_handle cof: face.coface_range(cod_d+1)) {
+	      auto qrb = oracle.intersects_boundary(cof, triangulation);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+	      mt_seed_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
+#endif	
+	      if (qrb.success)
+		boundary_simplex_map.emplace(cof, qrb.intersection);
+	    }
+	  }
+	  // break;
+	}
+      }
+    }
+    
+    while (!queue.empty()) {
+      Simplex_handle s = queue.front();
+      queue.pop();
+      for (auto cof: s.coface_range(cod_d+1)) {
+	for (auto face: cof.face_range(cod_d)) {
+	  auto qr = oracle.intersects(face, triangulation);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+	  mt_inserted_list.push_back(MT_inserted_info(qr, face, false));
+#endif	
+	  if (qr.success) {
+	    if (oracle.lies_in_domain(qr.intersection, triangulation)) {
+	      if (interior_simplex_map.emplace(face, qr.intersection).second)
+		queue.emplace(face);
+	    }
+	    else {
+	      auto qrb = oracle.intersects_boundary(cof, triangulation);
+#ifdef GUDHI_COX_OUTPUT_TO_HTML
+	      mt_inserted_list.push_back(MT_inserted_info(qrb, cof, true));
+#endif	
+	      // assert (qrb.success); // always a success
+	      if (qrb.success)
+		boundary_simplex_map.emplace(cof, qrb.intersection);
+	    }
+	  }
+	}
+      }
+    }
+  }
+
+  /** \brief Empty constructor */
+  Manifold_tracing() {}
+  
+};
+
+/**
+ * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm
+ * that computes the set of k-simplices that intersect
+ * a boundaryless implicit manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the 
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ *  \tparam Triangulation_ The type of the ambient triangulation.
+ *   Needs to be a model of the concept TriangulationForManifoldTracing.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ *  Needs to be a model of the concept IntersectionOracle.
+ * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map.
+ *
+ * \param[in] seed_points The range of points on the manifold from which 
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] out_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the input manifold and the mapped values 
+ * are the intersection points.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Point_range,
+	  class Triangulation,
+	  class Intersection_oracle,
+	  class Out_simplex_map>
+void manifold_tracing_algorithm(const Point_range& seed_points,
+				const Triangulation& triangulation,
+				const Intersection_oracle& oracle,
+				Out_simplex_map& out_simplex_map) {
+  Manifold_tracing<Triangulation> mt;
+  mt.manifold_tracing_algorithm(seed_points,
+				triangulation,
+				oracle,
+				out_simplex_map);
+}
+
+/**
+ * \brief Static method for Manifold_tracing<Triangulation_>::manifold_tracing_algorithm
+ * the dimensional manifold given by an intersection oracle, where k
+ * is the codimension of the manifold.
+ * The computation is based on the seed propagation --- it starts at the 
+ * given seed points and then propagates along the manifold.
+ *
+ * \tparam Point_range Range of points of type Eigen::VectorXd.
+ *  \tparam Triangulation_ The type of the ambient triangulation.
+ *   Needs to be a model of the concept TriangulationForManifoldTracing.
+ * \tparam Intersection_oracle Intersection oracle that represents the manifold.
+ *  Needs to be a model of the concept IntersectionOracle.
+ * \tparam Out_simplex_map Needs to be Manifold_tracing<Triangulation_>::Out_simplex_map.
+ *
+ * \param[in] seed_points The range of points on the manifold from which 
+ * the computation begins.
+ * \param[in] triangulation The ambient triangulation.
+ * \param[in] oracle The intersection oracle for the manifold.
+ * The ambient dimension needs to match the dimension of the
+ * triangulation.
+ * \param[out] interior_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the relative interior of the input manifold  
+ * and the mapped values are the intersection points.
+ * \param[out] boundary_simplex_map The output map, where the keys are k-simplices in
+ * the input triangulation that intersect the boundary of the input manifold  
+ * and the mapped values are the intersection points.
+ *
+ * \ingroup coxeter_triangulation
+ */
+template <class Point_range,
+	  class Triangulation,
+	  class Intersection_oracle,
+	  class Out_simplex_map>
+void manifold_tracing_algorithm(const Point_range& seed_points,
+				const Triangulation& triangulation,
+				const Intersection_oracle& oracle,
+				Out_simplex_map& interior_simplex_map,
+				Out_simplex_map& boundary_simplex_map) {
+  Manifold_tracing<Triangulation> mt;
+  mt.manifold_tracing_algorithm(seed_points,
+				triangulation,
+				oracle,
+				interior_simplex_map,
+				boundary_simplex_map);
+}
+
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h
new file mode 100644
index 00000000..a3c2b2c7
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation.h
@@ -0,0 +1,257 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_H_
+#define PERMUTAHEDRAL_REPRESENTATION_H_
+
+#include <gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h>
+
+#include <unordered_set>
+#include <iostream>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** 
+ * \class Permutahedral_representation
+ * \brief A class that stores the permutahedral representation of a simplex
+ * in a Coxeter triangulation or a Freudenthal-Kuhn triangulation.
+ *
+ * \ingroup coxeter_triangulation
+ *
+ * \details The data structure is a record consisting of a range that
+ * represents the vertex and a range that represents the ordered set
+ * partition, both of which identify the simplex in the triangulation.
+ *
+ * \tparam Vertex_ needs to be a random-access range.
+ * \tparam Ordered_set_partition_ needs to be a a random-access range that consists of
+ * random-access ranges. 
+ */
+template <class Vertex_,
+	  class Ordered_set_partition_>
+class Permutahedral_representation {
+
+  typedef Permutahedral_representation<Vertex_, Ordered_set_partition_> Self;
+
+public:
+  /** \brief Type of the vertex. */
+  typedef Vertex_ Vertex;
+  
+  /** \brief Type of the ordered partition. */
+  typedef Ordered_set_partition_ OrderedSetPartition;
+
+  
+  /** \brief Permutahedral_representation constructor from a vertex and an ordered set partition.
+   *
+   * @param[in] vertex Vertex.
+   * @param[in] partition Ordered set partition.
+   * 
+   * \details If the size of vertex is d, the ranges in partition must consist
+   * of the integers 0,...,d without repetition or collision between the ranges.
+   */
+  Permutahedral_representation(const Vertex& vertex, const OrderedSetPartition& partition)
+    : vertex_(vertex), partition_(partition) {}
+
+  /** \brief Constructor for an empty permutahedral representation that does not correspond
+   *  to any simplex.
+   */
+  Permutahedral_representation() {}
+  
+  /** \brief Dimension of the simplex. */
+  std::size_t dimension() const {
+    return partition_.size() - 1;
+  }
+
+  /** \brief Lexicographically-minimal vertex. */
+  Vertex& vertex() {
+    return vertex_;
+  }
+
+  /** \brief Lexicographically-minimal vertex. */
+  const Vertex& vertex() const {
+    return vertex_;
+  }
+
+  /** \brief Ordered set partition. */
+  OrderedSetPartition& partition() {
+    return partition_;
+  }
+
+  /** \brief Identifying vertex. */
+  const OrderedSetPartition& partition() const {
+    return partition_;
+  }
+
+  /** \brief Equality operator.
+   * Returns true if an only if both vertex and the ordered set partition coincide.
+   */
+  bool operator==(const Permutahedral_representation& other) const {
+    if (dimension() != other.dimension())
+      return false;
+    if (vertex_ != other.vertex_)
+      return false;
+    for (std::size_t k = 0; k < partition_.size(); ++k)
+      if (partition_[k] != other.partition_[k])
+	return false;
+    return true;
+  }
+
+  /** \brief Inequality operator.
+   * Returns true if an only if either vertex or the ordered set partition are different.
+   */
+  bool operator!=(const Permutahedral_representation& other) const {
+    return !(*this == other);
+  }
+
+  typedef Gudhi::coxeter_triangulation::Vertex_iterator<Self> Vertex_iterator;
+  typedef boost::iterator_range<Vertex_iterator> Vertex_range;
+  /** \brief Returns a range of vertices of the simplex.
+   * The type of vertices is Vertex.
+   */
+  Vertex_range vertex_range() const {
+    return Vertex_range(Vertex_iterator(*this),
+  			Vertex_iterator());
+  }
+
+  typedef Gudhi::coxeter_triangulation::Face_iterator<Self> Face_iterator;
+  typedef boost::iterator_range<Face_iterator> Face_range;
+  /** \brief Returns a range of permutahedral representations of faces of the simplex.
+   * @param[in] value_dim The dimension of the faces. Must be between 0 and the dimension of the simplex.
+   */
+  Face_range face_range(std::size_t value_dim) const {
+    return Face_range(Face_iterator(*this, value_dim),
+		      Face_iterator());
+  }
+
+  /** \brief Returns a range of permutahedral representations of facets of the simplex.
+   * The dimension of the simplex must be strictly positive.
+   */
+  Face_range facet_range() const {
+    return Face_range(Face_iterator(*this, dimension()-1),
+		      Face_iterator());
+  }
+  
+  typedef Gudhi::coxeter_triangulation::Coface_iterator<Self> Coface_iterator;
+  typedef boost::iterator_range<Coface_iterator> Coface_range;
+  /** \brief Returns a range of permutahedral representations of cofaces of the simplex.
+   * @param[in] value_dim The dimension of the cofaces. Must be between the dimension of the simplex and the ambient dimension (the size of the vertex).
+   */
+  Coface_range coface_range(std::size_t value_dim) const {
+    return Coface_range(Coface_iterator(*this, value_dim),
+			Coface_iterator());
+  }
+
+  /** \brief Returns a range of permutahedral representations of cofacets of the simplex.
+   * The dimension of the simplex must be strictly different from the ambient dimension (the size of the vertex).
+   */
+  Coface_range cofacet_range() const {
+    return Coface_range(Coface_iterator(*this, dimension()+1),
+			Coface_iterator());
+  }
+
+  /** \brief Returns true, if the simplex is a face of other simplex.
+   *
+   * @param[in] other A simplex that is potential a coface of the current simplex.
+   */
+  bool is_face_of(const Permutahedral_representation& other) const {
+    using Part = typename OrderedSetPartition::value_type;
+      
+    if (other.dimension() < dimension())
+      return false;
+    if (other.vertex_.size() != vertex_.size())
+      std::cerr << "Error: Permutahedral_representation::is_face_of: incompatible ambient dimensions.\n";
+    
+    Vertex v_self = vertex_, v_other = other.vertex_;
+    auto self_partition_it = partition_.begin();
+    auto other_partition_it = other.partition_.begin();
+    while (self_partition_it != partition_.end()) {
+      while (other_partition_it != other.partition_.end() && v_self != v_other) {
+	const Part& other_part = *other_partition_it++;
+	if (other_partition_it == other.partition_.end())
+	  return false;
+	for (const auto& k: other_part)
+	  v_other[k]++;	  
+      }
+      if (other_partition_it == other.partition_.end())
+	return false;
+      const Part& self_part = *self_partition_it++;
+      if (self_partition_it == partition_.end())
+	return true;
+      for (const auto& k: self_part)
+	v_self[k]++;
+    }
+    return true;
+  }
+
+private:
+  Vertex vertex_;
+  OrderedSetPartition partition_;
+
+};
+
+/** \brief Print a permutahedral representation to a stream.
+ * \ingroup coxeter_triangulation
+ *
+ * @param[in] os The output stream.
+ * @param[in] simplex A simplex represented by its permutahedral representation.
+ */
+template <class Vertex,
+	  class OrderedSetPartition>
+std::ostream& operator<<(std::ostream& os,
+			 const Permutahedral_representation<Vertex, OrderedSetPartition>& simplex) {
+  // vertex part
+  os << "(";
+  if (simplex.vertex().empty()) {
+    os << ")";
+    return os;
+  }
+  auto v_it = simplex.vertex().begin();
+  os << *v_it++;
+  for (; v_it != simplex.vertex().end(); ++v_it)
+    os << ", " << *v_it;
+  os << ")";
+  
+  // ordered partition part
+  using Part = typename OrderedSetPartition::value_type;
+  auto print_part =
+    [&os](const Part& p) {
+      os << "{";
+      if (p.empty()) {
+	os << "}";
+      }
+      auto p_it = p.begin();
+      os << *p_it++;
+      for (; p_it != p.end(); ++p_it)
+	os << ", " << *p_it;  
+      os << "}";
+    };
+  os << " [";
+  if (simplex.partition().empty()) {
+    os << "]";
+    return os;
+  }
+  auto o_it = simplex.partition().begin();
+  print_part(*o_it++);
+  for (; o_it != simplex.partition().end(); ++o_it) {
+    os << ", ";
+    print_part(*o_it);
+  }
+  os << "]";
+  return os;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif // PERMUTAHEDRAL_REPRESENTATION_H_
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h
new file mode 100644
index 00000000..20a4a8ff
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Combination_iterator.h
@@ -0,0 +1,101 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_COMBINATION_ITERATOR_H_
+
+#include <vector>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate combinations of 
+ *   k elements in a set of n elements.
+ *  Based on the algorithm by Mifsud.
+*/
+class Combination_iterator : public boost::iterator_facade< Combination_iterator,
+							    std::vector<uint> const,
+							    boost::forward_traversal_tag> {
+  typedef std::vector<uint> value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+  
+  bool equal(Combination_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    if (value_[0] == n_ - k_) {
+      is_end_ = true;
+      return;
+    }
+    uint j = k_ - 1;
+    if (value_[j] < n_ - 1) {
+      value_[j]++;
+      return;
+    }
+    for (; j > 0; --j)
+      if (value_[j-1] < n_ - k_ + j-1) {
+	value_[j-1]++;
+	for (uint s = j; s < k_; s++)
+	  value_[s] = value_[j-1] + s - (j-1);
+	return;
+      }
+  }
+
+public:
+  
+  Combination_iterator(const uint& n, const uint& k)
+    :
+    value_(k),
+    is_end_(n == 0),
+    n_(n),
+    k_(k)
+  {
+    for (uint i = 0; i < k; ++i)
+      value_[i] = i;
+  }
+
+  // Used for the creating an end iterator
+  Combination_iterator() : is_end_(true), n_(0), k_(0) {}
+
+  void reinitialize() {
+    if (n_ > 0) {
+      is_end_ = false;
+      for (uint i = 0; i < n_; ++i)
+      value_[i] = i;
+    }
+  }
+  
+protected:
+  value_t value_; // the dereference value
+  bool is_end_;   // is true when the current permutation is the final one 
+
+  uint n_;
+  uint k_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h
new file mode 100644
index 00000000..47eb8b98
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Integer_combination_iterator.h
@@ -0,0 +1,128 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_INTEGER_COMBINATION_ITERATOR_H_
+
+#include <vector>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate combinations of 
+ *   k elements in a set of n elements.
+ *  Based on the algorithm by Mifsud.
+*/
+class Integer_combination_iterator : public boost::iterator_facade< Integer_combination_iterator,
+								    std::vector<uint> const,
+								    boost::forward_traversal_tag> {
+  typedef std::vector<uint> value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+  
+  bool equal(Integer_combination_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    uint j1 = 0;
+    uint s = 0;
+    while (value_[j1] == 0 && j1 < k_)
+      j1++;
+    uint j2 = j1+1;
+    while (value_[j2] == bounds_[j2]) {
+      if (bounds_[j2] != 0) {
+	s += value_[j1];
+	value_[j1] = 0;
+	j1 = j2;
+      }
+      j2++;
+    }
+    if (j2 >= k_) {
+      is_end_ = true;
+      return;
+    }
+    s += value_[j1] - 1;
+    value_[j1] = 0;
+    value_[j2]++;
+    uint i = 0;
+    while (s >= bounds_[i]) {
+      value_[i] = bounds_[i];
+      s -= bounds_[i];
+      i++;
+    }
+    value_[i++] = s;
+  }
+
+public:
+
+  template <class Bound_range>
+  Integer_combination_iterator(const uint& n, const uint& k, const Bound_range& bounds)
+    :
+    value_(k+2),
+    is_end_(n == 0 || k == 0),
+    n_(n),
+    k_(k)
+  {
+    bounds_.reserve(k+2);
+    uint sum_radices = 0;
+    for (auto b: bounds) {
+      bounds_.push_back(b);
+      sum_radices += b;      
+    }
+    bounds_.push_back(2);
+    bounds_.push_back(1);
+    if (n > sum_radices) {
+      is_end_ = true;
+      return;
+    }
+    uint i = 0;
+    uint s = n;
+    while (s >= bounds_[i]) {
+      value_[i] = bounds_[i];
+      s -= bounds_[i];
+      i++;
+    }
+    value_[i++] = s;
+
+    while (i < k_)
+      value_[i++] = 0;    
+    value_[k] = 1;
+    value_[k+1] = 0;
+  }
+
+  // Used for the creating an end iterator
+  Integer_combination_iterator() : is_end_(true), n_(0), k_(0) {}
+  
+protected:
+  value_t value_; // the dereference value
+  bool is_end_;   // is true when the current integer combination is the final one 
+
+  uint n_;
+  uint k_;
+  std::vector<uint> bounds_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h
new file mode 100644
index 00000000..55c32664
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h
@@ -0,0 +1,108 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_ORDERED_SET_PARTITION_ITERATOR_H_
+
+#include <vector>
+#include <limits>
+#include <gudhi/Permutahedral_representation/Permutation_iterator.h>
+#include <gudhi/Permutahedral_representation/Set_partition_iterator.h>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that represents an ordered set partition of a set {0,...,n-1} in k parts as
+ *         a pair of an unordered set partition given in lexicographic order and
+ *         a permutation of the parts.
+ */
+struct Ordered_set_partition {
+  Set_partition_iterator s_it_;
+  Permutation_iterator p_it_;
+
+  // Ordered_set_partition(const Set_partition_iterator& s_it, const Permutation_iterator& p_it)
+  //   : s_it_(s_it), p_it_(p_it) {}
+  
+  const std::vector<uint> operator[](const uint& i) const {
+    return (*s_it_)[(*p_it_)[i]];
+  }
+  
+  std::size_t size() const {
+    return s_it_->size();
+  }
+  
+};
+
+/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts.
+ *   
+*/
+class Ordered_set_partition_iterator : public boost::iterator_facade< Ordered_set_partition_iterator,
+								      Ordered_set_partition const,
+								      boost::forward_traversal_tag> {
+  typedef Ordered_set_partition value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+
+  bool equal(Ordered_set_partition_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    if (++value_.p_it_ == p_end_) {
+      if (++value_.s_it_ == s_end_) {
+	is_end_ = true;
+	return;
+      }
+      else
+	value_.p_it_.reinitialize();
+    }
+  }
+  
+public:
+  
+  Ordered_set_partition_iterator(const uint& n, const uint& k)
+    :
+    value_({Set_partition_iterator(n,k), Permutation_iterator(k)}),
+    is_end_(n == 0)
+  {
+  }
+
+  // Used for the creating an end iterator
+  Ordered_set_partition_iterator() : is_end_(true) {}
+
+  void reinitialize() {
+    is_end_ = false;
+    value_.p_it_.reinitialize();
+    value_.s_it_.reinitialize();
+  }
+  
+protected:
+  Set_partition_iterator s_end_; // Set partition iterator and the corresponding end iterator
+  Permutation_iterator p_end_; // Permutation iterator and the corresponding end iterator
+  value_t value_;         // the dereference value
+  bool is_end_;   // is true when the current permutation is the final one 
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h
new file mode 100644
index 00000000..4528ef20
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutahedral_representation_iterators.h
@@ -0,0 +1,288 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_
+#define PERMUTAHEDRAL_REPRESENTATION_PERMUTAHEDRAL_REPRESENTATION_ITERATORS_H_
+
+#include <gudhi/Permutahedral_representation/Size_range.h>
+#include <gudhi/Permutahedral_representation/Ordered_set_partition_iterator.h>
+#include <gudhi/Permutahedral_representation/Integer_combination_iterator.h>
+#include <gudhi/Permutahedral_representation/Combination_iterator.h>
+#include <gudhi/Permutahedral_representation/face_from_indices.h>
+#include <boost/iterator/iterator_facade.hpp>
+
+#include <vector>
+#include <iostream>
+#include <algorithm>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/* \addtogroup coxeter_triangulation
+ * Iterator types for Permutahedral_representation
+ * @{
+ */
+
+/* \brief Iterator over the vertices of a simplex
+ * represented by its permutahedral representation.
+ *
+ * Forward iterator, 'value_type' is Permutahedral_representation::Vertex.*/
+template <class Permutahedral_representation>
+class Vertex_iterator : public boost::iterator_facade< Vertex_iterator<Permutahedral_representation>,
+						       typename Permutahedral_representation::Vertex const,
+						       boost::forward_traversal_tag> {
+protected:
+  friend class boost::iterator_core_access;
+
+  using Vertex = typename Permutahedral_representation::Vertex;
+  using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+  typedef Vertex value_t;
+
+  
+  bool equal(Vertex_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void update_value() {
+    std::size_t d = value_.size();
+    for (auto i: *o_it_)
+      if (i != d)
+	value_[i]++;
+      else
+	for (std::size_t j = 0; j < d; j++)
+	  value_[j]--;
+  }
+		       
+  void increment() {
+    if (is_end_)
+      return;
+    update_value();
+    if (++o_it_ == o_end_)
+      is_end_ = true;
+  }
+
+public:  
+  Vertex_iterator(const Permutahedral_representation& simplex)
+    :
+    o_it_(simplex.partition().begin()),
+    o_end_(simplex.partition().end()),
+    value_(simplex.vertex()),
+    is_end_(o_it_ == o_end_)
+  {}
+
+  Vertex_iterator() : is_end_(true) {}
+  
+  
+protected:
+  typename Ordered_partition::const_iterator o_it_, o_end_; 
+  value_t value_;
+  bool is_end_;
+
+}; // Vertex_iterator
+
+/*---------------------------------------------------------------------------*/
+/* \brief Iterator over the k-faces of a simplex
+ *  given by its permutahedral representation.
+ *
+ * Forward iterator, value_type is Permutahedral_representation. */
+template <class Permutahedral_representation>
+class Face_iterator : public boost::iterator_facade< Face_iterator<Permutahedral_representation>,
+						     Permutahedral_representation const,
+						     boost::forward_traversal_tag> {
+  typedef Permutahedral_representation value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+
+  using Vertex = typename Permutahedral_representation::Vertex;
+  using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+  bool equal(Face_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    if (++c_it_ == c_end_) {
+      is_end_ = true;
+      return;
+    }
+    update_value();
+  }
+
+  void update_value() {
+    // Combination *c_it_ is supposed to be sorted in increasing order
+    value_ = face_from_indices<Permutahedral_representation>(simplex_, *c_it_);
+  }
+  
+public:
+  
+  Face_iterator(const Permutahedral_representation& simplex, const uint& k)
+    : simplex_(simplex),
+      k_(k),
+      l_(simplex.dimension()),
+      c_it_(l_+1, k_+1),
+      is_end_(k_ > l_),
+      value_({Vertex(simplex.vertex().size()), Ordered_partition(k+1)})
+  {
+    update_value();
+  }
+
+  // Used for the creating an end iterator
+  Face_iterator() : is_end_(true) {}
+  
+protected:
+  Permutahedral_representation simplex_; // Input simplex
+  uint k_;
+  uint l_; // Dimension of the input simplex
+  Combination_iterator c_it_, c_end_; // indicates the vertices in the current face
+
+  bool is_end_;   // is true when the current permutation is the final one
+  value_t value_; // the dereference value
+
+}; // Face_iterator
+
+/*---------------------------------------------------------------------------*/
+/* \brief Iterator over the k-cofaces of a simplex
+ *  given by its permutahedral representation.
+ *
+ * Forward iterator, value_type is Permutahedral_representation. */
+template <class Permutahedral_representation>
+class Coface_iterator : public boost::iterator_facade< Coface_iterator<Permutahedral_representation>,
+						       Permutahedral_representation const,
+						       boost::forward_traversal_tag> {
+  typedef Permutahedral_representation value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+
+  using Vertex = typename Permutahedral_representation::Vertex;
+  using Ordered_partition = typename Permutahedral_representation::OrderedSetPartition;
+
+  bool equal(Coface_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    uint i = 0;
+    for (; i < k_+1; i++) {
+      if (++(o_its_[i]) != o_end_)
+	break;
+    }
+    if (i == k_+1) {
+      if (++i_it_ == i_end_) {
+	is_end_ = true;
+	return;
+      }
+      o_its_.clear();
+      for (uint j = 0; j < k_ + 1; j++)
+	o_its_.emplace_back(Ordered_set_partition_iterator(simplex_.partition()[j].size(), (*i_it_)[j]+1));
+    }
+    else
+      for (uint j = 0; j < i; j++)
+	o_its_[j].reinitialize();
+    update_value();
+  }
+
+  void update_value() {
+    value_.vertex() = simplex_.vertex();
+    for (auto& p: value_.partition())
+      p.clear();
+    uint u_ = 0; // the part in o_its_[k_] that contains t_
+    for (; u_ <= (*i_it_)[k_]; u_++) {
+      auto range = (*o_its_[k_])[u_];
+      if (std::find(range.begin(), range.end(), t_) != range.end())
+	break;
+    }
+    uint i = 0;
+    for (uint j = u_+1; j <= (*i_it_)[k_]; j++, i++)
+      for (uint b: (*o_its_[k_])[j]) {
+	uint c = simplex_.partition()[k_][b];
+        value_.partition()[i].push_back(c);
+	value_.vertex()[c]--;
+      }
+    for (uint h = 0; h < k_; h++)
+      for (uint j = 0; j <= (*i_it_)[h]; j++, i++) {
+	for (uint b: (*o_its_[h])[j])
+	  value_.partition()[i].push_back(simplex_.partition()[h][b]);
+      }
+    for (uint j = 0; j <= u_; j++, i++)
+      for (uint b: (*o_its_[k_])[j])
+	value_.partition()[i].push_back(simplex_.partition()[k_][b]);
+    // sort the values in each part (probably not needed)
+    for (auto& part: value_.partition())
+      std::sort(part.begin(), part.end());
+  }
+  
+public:
+  
+  Coface_iterator(const Permutahedral_representation& simplex, const uint& l)
+    : simplex_(simplex),
+      d_(simplex.vertex().size()),
+      l_(l),
+      k_(simplex.dimension()),
+      i_it_(l_-k_ , k_+1, Size_range<Ordered_partition>(simplex.partition())),
+      is_end_(k_ > l_),
+      value_({Vertex(d_), Ordered_partition(l_+1)})
+  {
+    uint j = 0;
+    for (; j < simplex_.partition()[k_].size(); j++)
+      if (simplex_.partition()[k_][j] == d_) {
+	t_ = j;
+	break;
+      }
+    if (j == simplex_.partition()[k_].size()) {
+      std::cerr << "Coface iterator: the argument simplex is not a permutahedral representation\n";
+      is_end_ = true;
+      return;
+    }
+    for (uint i = 0; i < k_+1; i++)
+      o_its_.emplace_back(Ordered_set_partition_iterator(simplex_.partition()[i].size(), (*i_it_)[i]+1));
+    update_value();
+  }
+
+  // Used for the creating an end iterator
+  Coface_iterator() : is_end_(true) {}
+  
+protected:
+  Permutahedral_representation simplex_; // Input simplex
+  uint d_; // Ambient dimension
+  uint l_; // Dimension of the coface
+  uint k_; // Dimension of the input simplex
+  uint t_; // The position of d in simplex_.partition()[k_]
+  Integer_combination_iterator i_it_, i_end_; // indicates in how many parts each simplex_[i] is subdivided
+  std::vector<Ordered_set_partition_iterator> o_its_; // indicates subdivision for each simplex_[i]
+  Ordered_set_partition_iterator o_end_;      // one end for all o_its_
+
+  bool is_end_;   // is true when the current permutation is the final one
+  value_t value_; // the dereference value
+
+}; // Coface_iterator
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+  
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h
new file mode 100644
index 00000000..7cf6158b
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Permutation_iterator.h
@@ -0,0 +1,137 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_PERMUTATION_ITERATOR_H_
+
+#include <vector>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate permutations.
+ *   Based on the optimization of the Heap's algorithm by Sedgewick.
+*/
+class Permutation_iterator : public boost::iterator_facade< Permutation_iterator,
+							    std::vector<uint> const,
+							    boost::forward_traversal_tag> {
+  typedef std::vector<uint> value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+
+  bool equal(Permutation_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void swap_two_indices(std::size_t i, std::size_t j) {
+    uint t = value_[i];
+    value_[i] = value_[j];
+    value_[j] = t;
+  }
+  
+  void elementary_increment() {
+    uint j = 0;
+    while (d_[j] == j+1) {
+      d_[j] = 0;
+      ++j;
+    }
+    if (j == n_ - 1) {
+      is_end_ = true;
+      return;
+    }
+    uint k = j+1;
+    uint x = (k%2 ? d_[j] : 0);
+    swap_two_indices(k, x);
+    ++d_[j];
+  }
+
+  void elementary_increment_optim_3() {
+    if (ct_ != 0) {
+      --ct_;
+      swap_two_indices(1 + (ct_%2), 0);
+    }
+    else {
+      ct_ = 5;
+      uint j = 2;
+      while (d_[j] == j+1) {
+	d_[j] = 0;
+	++j;
+      }
+      if (j == n_ - 1) {
+	is_end_ = true;
+	return;
+      }
+      uint k = j+1;
+      uint x = (k%2 ? d_[j] : 0);
+      swap_two_indices(k, x);
+      ++d_[j];
+    }
+  }
+    
+  void increment() {
+    if (optim_3_)
+      elementary_increment_optim_3();
+    else
+      elementary_increment();      
+  }
+
+public:
+  
+  Permutation_iterator(const uint& n)
+    :
+    value_(n),
+    is_end_(n == 0),
+    optim_3_(n >= 3),
+    n_(n),
+    d_(n),
+    ct_(5)
+  {
+    for (uint i = 0; i < n; ++i) {
+      value_[i] = i;
+      d_[i] = 0;
+    }
+    if (n > 0)
+      d_[n-1] = -1;
+  }
+
+  // Used for the creating an end iterator
+  Permutation_iterator() : is_end_(true), n_(0) {}
+
+  void reinitialize() {
+    if (n_ > 0)
+      is_end_ = false;
+  }
+  
+protected:
+  value_t value_; // the dereference value
+  bool is_end_;   // is true when the current permutation is the final one 
+  bool optim_3_;  // true if n>=3. for n >= 3, the algorithm is optimized
+
+  uint n_;
+  std::vector<uint> d_; // mix radix digits with radix [2 3 4 ... n-1 (sentinel=-1)]
+  uint ct_;             // counter with values in {0,...,5} used in the n>=3 optimization.
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h
new file mode 100644
index 00000000..46f67752
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Set_partition_iterator.h
@@ -0,0 +1,137 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SET_PARTITION_ITERATOR_H_
+
+#include <vector>
+#include <limits>
+#include <boost/range/iterator_range.hpp>
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+typedef unsigned uint;
+
+/** \brief Class that allows the user to generate set partitions of a set {0,...,n-1} in k parts.
+ *   
+*/
+class Set_partition_iterator : public boost::iterator_facade< Set_partition_iterator,
+							      std::vector<std::vector<uint> > const,
+							      boost::forward_traversal_tag> {
+  typedef std::vector<std::vector<uint> > value_t;
+  
+protected:
+  friend class boost::iterator_core_access;
+
+  bool equal(Set_partition_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  value_t const& dereference() const {
+    return value_;
+  }
+
+  void update_value() {
+    for (uint i = 0; i < k_; i++)
+      value_[i].clear();
+    for (uint i = 0; i < n_; i++)
+      value_[rgs_[i]].push_back(i);
+  }
+  
+  void increment() {
+    if (k_ <= 1) {
+      is_end_ = true;
+      return;
+    }
+    uint i = n_ - 1;
+    while (rgs_[i] + 1 > max_[i] ||
+	   rgs_[i] + 1 >= k_)
+      i--;
+    if (i == 0) {
+      is_end_ = true;
+      return;
+    }
+    rgs_[i]++;
+    uint mm = max_[i];
+    mm += (rgs_[i] >= mm);
+    max_[i+1] = mm;
+    while (++i < n_) {
+      rgs_[i] = 0;
+      max_[i+1] = mm;
+    }
+    uint p = k_;
+    if (mm < p)
+      do {
+	max_[i] = p;
+	--i;
+	--p;
+	rgs_[i] = p;
+      } while (max_[i] < p);
+    update_value();
+  }
+  
+public:
+  
+  Set_partition_iterator(const uint& n, const uint& k)
+    :
+    value_(k),
+    rgs_(n, 0),
+    max_(n+1),
+    is_end_(n == 0),
+    n_(n),
+    k_(k)
+  {
+    max_[0] = std::numeric_limits<uint>::max();
+    for (uint i = 0; i <= n-k; ++i)
+      value_[0].push_back(i);
+    for (uint i = n-k+1, j = 1; i < n; ++i, ++j) {
+      rgs_[i] = j;
+      value_[j].push_back(i);
+    }
+    for (uint i = 1; i <= n; i++)
+      max_[i] = rgs_[i-1] + 1;
+    update_value();
+  }
+
+  // Used for creating an end iterator
+  Set_partition_iterator() : is_end_(true), n_(0), k_(0) {}
+
+  void reinitialize() {
+    if (n_ > 0)
+      is_end_ = false;
+    for (uint i = 0; i <= n_-k_; ++i)
+      rgs_[i] = 0;
+    for (uint i = n_-k_+1, j = 1; i < n_; ++i, ++j)
+      rgs_[i] = j;
+    for (uint i = 1; i <= n_; i++)
+      max_[i] = rgs_[i-1] + 1;
+    update_value();
+  }
+  
+protected:
+  value_t value_;         // the dereference value
+  std::vector<uint> rgs_; // restricted growth string
+  std::vector<uint> max_; // max_[i] = max(rgs_[0],...,rgs[i-1]) + 1
+  bool is_end_;   // is true when the current permutation is the final one 
+
+  uint n_;
+  uint k_;
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h
new file mode 100644
index 00000000..5b3c29e9
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Simplex_comparator.h
@@ -0,0 +1,64 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SIMPLEX_COMPARATOR_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Simplex_comparator 
+ *  \brief A comparator class for Permutahedral_representation.
+ *   The comparison is in lexicographic order first on
+ *   vertices and then on ordered partitions with sorted parts.
+ *   The lexicographic order forces that any face is larger than
+ *   a coface.
+ *
+ *  \tparam Permutahdral_representation_ Needs to be 
+ *   Permutahedral_representation<Vertex_, Ordered_set_partition_>
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template <class Permutahedral_representation_>
+struct Simplex_comparator {
+
+  /** \brief Comparison between two permutahedral representations.
+   *  Both permutahedral representations need to be valid and
+   *  the vertices of both permutahedral representations need to be of the same size.
+   */
+  bool operator() (const Permutahedral_representation_& lhs,
+		   const Permutahedral_representation_& rhs) const {
+    if (lhs.vertex() < rhs.vertex())
+      return true;
+    if (lhs.vertex() > rhs.vertex())
+      return false;
+    
+    if (lhs.partition().size() > rhs.partition().size())
+      return true;
+    if (lhs.partition().size() < rhs.partition().size())
+      return false;
+    
+    if (lhs.partition() < rhs.partition())
+      return true;
+    if (lhs.partition() > rhs.partition())
+      return false;
+    
+    return false;
+  }
+};
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h
new file mode 100644
index 00000000..dd9c20dc
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/Size_range.h
@@ -0,0 +1,89 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_
+#define PERMUTAHEDRAL_REPRESENTATION_SIZE_RANGE_H_
+
+#include <boost/range/iterator_range.hpp>
+#include <cstdlib>         // std::size_t
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \brief Auxillary iterator class for sizes of parts in an ordered set partition.
+ */
+template <class T_it>
+class Size_iterator : public boost::iterator_facade< Size_iterator<T_it>,
+						     std::size_t const,
+						     boost::forward_traversal_tag> {
+  friend class boost::iterator_core_access;
+
+protected:
+  bool equal(Size_iterator const& other) const {
+    return (is_end_ && other.is_end_);
+  }
+
+  std::size_t const& dereference() const {
+    return value_;
+  }
+
+  void increment() {
+    if (++t_it_ == t_end_) {
+      is_end_ = true;
+      return;
+    }
+    value_ = t_it_->size()-1;
+  }
+
+public:
+
+  Size_iterator(const T_it& t_begin, const T_it& t_end)
+    : t_it_(t_begin), t_end_(t_end), is_end_(t_begin == t_end) {
+    if (!is_end_)
+      value_ = t_it_->size()-1;
+  }
+
+protected:
+  T_it t_it_, t_end_;
+  bool is_end_;
+  std::size_t value_;
+};
+
+template <class T>
+class Size_range {
+  const T& t_;
+
+public:
+  typedef Size_iterator<typename T::const_iterator> iterator;
+  
+  Size_range(const T& t) : t_(t) {}
+  
+  std::size_t operator[](std::size_t i) const {
+    return t_[i].size()-1;
+  }
+
+  iterator begin() const {
+    return iterator(t_.begin(), t_.end());
+  }
+
+  iterator end() const {
+    return iterator(t_.end(), t_.end());
+  }
+
+};
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h
new file mode 100644
index 00000000..461a01e7
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Permutahedral_representation/face_from_indices.h
@@ -0,0 +1,71 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_
+#define PERMUTAHEDRAL_REPRESENTATION_FACE_FROM_INDICES_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+ /** \brief Computes the permutahedral representation of a face of a given simplex
+  *  and a range of the vertex indices that compose the face.
+  *
+  * \tparam Permutahedral_representation has to be Permutahedral_representation
+  * \tparam Index_range is a range of unsigned integers taking values in 0,...,k,
+  * where k is the dimension of the simplex simplex.
+  *
+  * @param[in] simplex Input simplex.
+  * @param[in] indices Input range of indices.
+  */
+template <class Permutahedral_representation,
+	  class Index_range>
+Permutahedral_representation face_from_indices(const Permutahedral_representation& simplex,
+					       const Index_range& indices) {
+  using range_index = typename Index_range::value_type;
+  using Ordered_set_partition = typename Permutahedral_representation::OrderedSetPartition;
+  using Part = typename Ordered_set_partition::value_type;
+  using part_index = typename Part::value_type;
+  Permutahedral_representation value;
+  std::size_t d = simplex.vertex().size();
+  value.vertex() = simplex.vertex();
+  std::size_t k = indices.size()-1;
+  value.partition().resize(k+1);
+  std::size_t l = simplex.partition().size()-1;
+  for (std::size_t h = 1; h < k+1; h++)
+    for (range_index i = indices[h-1]; i < indices[h]; i++)
+      for (part_index j: simplex.partition()[i])
+	value.partition()[h-1].push_back(j);
+  for (range_index i = indices[k]; i < l+1; i++)
+    for (part_index j: simplex.partition()[i])
+      value.partition()[k].push_back(j);
+  for (range_index i = 0; i < indices[0]; i++)
+    for (part_index j: simplex.partition()[i]) {
+      if (j != d)
+	value.vertex()[j]++;
+      else
+	for (std::size_t l = 0; l < d; l++)
+	  value.vertex()[l]--;
+      value.partition()[k].push_back(j);
+    }
+  // sort the values in each part (probably not needed)
+  for (auto& part: value.partition())
+    std::sort(part.begin(), part.end());
+  return value;
+}
+
+} // namespace coxeter_triangulation
+
+} // namespace Gudhi
+
+
+#endif
diff --git a/src/Coxeter_triangulation/include/gudhi/Query_result.h b/src/Coxeter_triangulation/include/gudhi/Query_result.h
new file mode 100644
index 00000000..c7384c0b
--- /dev/null
+++ b/src/Coxeter_triangulation/include/gudhi/Query_result.h
@@ -0,0 +1,42 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#ifndef QUERY_RESULT_H_
+#define QUERY_RESULT_H_
+
+namespace Gudhi {
+
+namespace coxeter_triangulation {
+
+/** \class Query_result
+ *  \brief The result of a query by an oracle such as Implicit_manifold_intersection_oracle. 
+ *
+ *  \tparam Simplex_handle The class of the query simplex.
+ *
+ *  \ingroup coxeter_triangulation
+ */
+template <class Simplex_handle>
+struct Query_result {
+  /** \brief The potentially lower-dimensional face of the query simplex
+   *   that contains the intersection point. OBSOLETE: as the snapping is removed. */
+  // Simplex_handle face;
+  /** \brief The intersection point. */
+  Eigen::VectorXd intersection;
+  /** \brief True if the query simplex intersects the manifold. */
+  bool success;
+};
+
+} // namespace coxeter_triangulation 
+
+} // namespace Gudhi
+
+#endif
diff --git a/src/Coxeter_triangulation/test/CMakeLists.txt b/src/Coxeter_triangulation/test/CMakeLists.txt
new file mode 100644
index 00000000..663768b1
--- /dev/null
+++ b/src/Coxeter_triangulation/test/CMakeLists.txt
@@ -0,0 +1,33 @@
+project(Coxeter_triangulation_test)
+
+include(GUDHI_test_coverage)
+
+
+add_executable ( Coxeter_triangulation_permutahedral_representation_test perm_rep_test.cpp )
+target_link_libraries( Coxeter_triangulation_permutahedral_representation_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} )
+gudhi_add_coverage_test(Coxeter_triangulation_permutahedral_representation_test)
+
+add_executable ( Coxeter_triangulation_freudenthal_triangulation_test freud_triang_test.cpp )
+target_link_libraries( Coxeter_triangulation_freudenthal_triangulation_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Coxeter_triangulation_freudenthal_triangulation_test)
+
+add_executable ( Coxeter_triangulation_functions_test function_test.cpp )
+target_link_libraries( Coxeter_triangulation_functions_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Coxeter_triangulation_functions_test)
+
+add_executable ( Coxeter_triangulation_oracle_test oracle_test.cpp )
+target_link_libraries( Coxeter_triangulation_oracle_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Coxeter_triangulation_oracle_test)
+
+add_executable ( Coxeter_triangulation_manifold_tracing_test manifold_tracing_test.cpp )
+target_link_libraries( Coxeter_triangulation_manifold_tracing_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Coxeter_triangulation_manifold_tracing_test)
+
+add_executable ( Coxeter_triangulation_cell_complex_test cell_complex_test.cpp )
+target_link_libraries( Coxeter_triangulation_cell_complex_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Coxeter_triangulation_cell_complex_test)
+
+if (TBB_FOUND)
+   target_link_libraries(Coxeter_triangulation_cell_complex_test ${TBB_LIBRARIES})
+endif()
+
diff --git a/src/Coxeter_triangulation/test/cell_complex_test.cpp b/src/Coxeter_triangulation/test/cell_complex_test.cpp
new file mode 100644
index 00000000..a48696a6
--- /dev/null
+++ b/src/Coxeter_triangulation/test/cell_complex_test.cpp
@@ -0,0 +1,62 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "cell_complex"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <gudhi/Debug_utils.h>
+#include <gudhi/IO/output_debug_traces_to_html.h>
+#include <iostream>
+
+#include <gudhi/Coxeter_triangulation.h>
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Function_torus_in_R3.h>
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+#include <gudhi/Manifold_tracing.h>
+#include <gudhi/Cell_complex.h>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(cell_complex) {
+
+  double radius = 1.1111;
+  Function_torus_in_R3 fun_torus(radius, 3*radius);
+  Eigen::VectorXd seed = fun_torus.seed();
+  Function_Sm_in_Rd fun_bound(2.5*radius, 2, seed);
+    
+  auto oracle = make_oracle(fun_torus, fun_bound);
+  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());
+  
+  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);
+  
+  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);
+  
+  std::size_t interior_sc_map_size0 = cell_complex.interior_simplex_cell_map(0).size();
+  std::size_t interior_sc_map_size1 = cell_complex.interior_simplex_cell_map(1).size();
+  std::size_t interior_sc_map_size2 = cell_complex.interior_simplex_cell_map(2).size();
+  std::size_t boundary_sc_map_size0 = cell_complex.boundary_simplex_cell_map(0).size();
+  std::size_t boundary_sc_map_size1 = cell_complex.boundary_simplex_cell_map(1).size();
+  BOOST_CHECK (interior_simplex_map.size() == interior_sc_map_size0 );
+  BOOST_CHECK ( boundary_sc_map_size0 - boundary_sc_map_size1 == 0 );
+  BOOST_CHECK ( interior_sc_map_size0 - interior_sc_map_size1 + interior_sc_map_size2 == 0 );
+}
diff --git a/src/Coxeter_triangulation/test/freud_triang_test.cpp b/src/Coxeter_triangulation/test/freud_triang_test.cpp
new file mode 100644
index 00000000..6e7bf865
--- /dev/null
+++ b/src/Coxeter_triangulation/test/freud_triang_test.cpp
@@ -0,0 +1,104 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "freudenthal_triangulation"
+#include <boost/test/unit_test.hpp>
+
+#include <gudhi/Unitary_tests_utils.h>
+#include <gudhi/Freudenthal_triangulation.h>
+#include <gudhi/Coxeter_triangulation.h>
+
+BOOST_AUTO_TEST_CASE(freudenthal_triangulation) {
+  
+  // Point location check
+  typedef std::vector<double> Point;
+  typedef Gudhi::coxeter_triangulation::Freudenthal_triangulation<> FK_triangulation;
+  typedef typename FK_triangulation::Simplex_handle Simplex_handle;
+  typedef typename FK_triangulation::Vertex_handle Vertex_handle;
+  typedef typename Simplex_handle::OrderedSetPartition Ordered_set_partition; 
+  typedef typename Ordered_set_partition::value_type Part;
+
+  FK_triangulation tr(3);
+
+  // Point location check
+  {
+    Point point({3, -1, 0});
+    Simplex_handle s = tr.locate_point(point);
+    BOOST_CHECK( s.vertex() == Vertex_handle({3, -1, 0}) );
+    BOOST_CHECK( s.partition() == Ordered_set_partition({Part({0, 1, 2, 3})}) );
+  }
+
+  {
+    Point point({3.5, -1.5, 0.5});
+    Simplex_handle s = tr.locate_point(point);
+    BOOST_CHECK( s.vertex() == Vertex_handle({3, -2, 0}) );
+    BOOST_CHECK( s.partition() == Ordered_set_partition({Part({0, 1, 2}), Part({3})}) );
+  }
+
+  {
+    Point point({3.5, -1.8, 0.5});
+    Simplex_handle s = tr.locate_point(point);
+    BOOST_CHECK( s.vertex() == Vertex_handle({3, -2, 0}) );
+    BOOST_CHECK( s.partition() == Ordered_set_partition({Part({0, 2}), Part({1}), Part({3})}) );
+  }
+
+  {
+    Point point({3.5, -1.8, 0.3});
+    Simplex_handle s = tr.locate_point(point);
+    BOOST_CHECK( s.vertex() == Vertex_handle({3, -2, 0}) );
+    BOOST_CHECK( s.partition() == Ordered_set_partition({Part({0}), Part({2}), Part({1}), Part({3})}) );
+  }
+
+  // Dimension check
+  BOOST_CHECK( tr.dimension() == 3 );
+  // Matrix check
+  Eigen::MatrixXd default_matrix = Eigen::MatrixXd::Identity(3, 3);  
+  BOOST_CHECK( tr.matrix() == default_matrix );
+  // Vector check
+  Eigen::MatrixXd default_offset = Eigen::VectorXd::Zero(3);  
+  BOOST_CHECK( tr.offset() == default_offset );
+
+  // Barycenter check
+  Point point({3.5, -1.8, 0.3});
+  Simplex_handle s = tr.locate_point(point);
+  Eigen::Vector3d barycenter_cart = Eigen::Vector3d::Zero();
+  for (auto v: s.vertex_range())
+    for (std::size_t i = 0; i < v.size(); i++)
+      barycenter_cart(i) += v[i];
+  barycenter_cart /= 4.; // simplex is three-dimensional
+  Eigen::Vector3d barycenter = tr.barycenter(s);
+  for (std::size_t i = 0; (long int)i < barycenter.size(); i++)
+    GUDHI_TEST_FLOAT_EQUALITY_CHECK(barycenter(i), barycenter_cart(i), 1e-7);
+
+  // Barycenter check for twice the scale
+  s = tr.locate_point(point, 2);
+  barycenter_cart = Eigen::Vector3d::Zero();
+  for (auto v: s.vertex_range())
+    for (std::size_t i = 0; i < v.size(); i++)
+      barycenter_cart(i) += v[i];
+  barycenter_cart /= 3.; // simplex is now a two-dimensional face
+  barycenter_cart /= 2.; // scale
+  barycenter = tr.barycenter(s, 2);
+  for (std::size_t i = 0; (long int)i < barycenter.size(); i++)
+    GUDHI_TEST_FLOAT_EQUALITY_CHECK(barycenter(i), barycenter_cart(i), 1e-7);
+  
+  // Matrix and offset change check
+  Eigen::MatrixXd new_matrix(3,3);
+  new_matrix << 1, 0, 0, -1, 1, 0, -1, 0, 1;
+  Eigen::Vector3d new_offset(1.5, 1, 0.5);
+  tr.change_matrix(new_matrix);
+  tr.change_offset(new_offset);
+  
+  BOOST_CHECK( tr.matrix() == new_matrix );
+  BOOST_CHECK( tr.offset() == new_offset );
+}
diff --git a/src/Coxeter_triangulation/test/function_test.cpp b/src/Coxeter_triangulation/test/function_test.cpp
new file mode 100644
index 00000000..518bef9b
--- /dev/null
+++ b/src/Coxeter_triangulation/test/function_test.cpp
@@ -0,0 +1,177 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+// workaround for the annoying boost message in boost 1.69
+#define BOOST_PENDING_INTEGER_LOG2_HPP
+#include <boost/integer/integer_log2.hpp>
+// end workaround 
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "function"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Function_affine_plane_in_Rd.h>
+#include <gudhi/Functions/Constant_function.h>
+#include <gudhi/Functions/Function_chair_in_R3.h>
+#include <gudhi/Functions/Function_torus_in_R3.h>
+#include <gudhi/Functions/Function_whitney_umbrella_in_R3.h>
+#include <gudhi/Functions/Function_lemniscate_revolution_in_R3.h>
+#include <gudhi/Functions/Function_iron_in_R3.h>
+#include <gudhi/Functions/Function_moment_curve_in_Rd.h>
+#include <gudhi/Functions/random_orthogonal_matrix.h>
+#include <gudhi/Functions/Embed_in_Rd.h>
+#include <gudhi/Functions/Translate.h>
+#include <gudhi/Functions/Linear_transformation.h>
+#include <gudhi/Functions/Negation.h>
+#include <gudhi/Functions/Cartesian_product.h>
+#include <gudhi/Functions/PL_approximation.h>
+
+#include <gudhi/Coxeter_triangulation.h>
+
+#include <string>
+
+#include <random>
+#include <cstdlib>
+
+using namespace Gudhi::coxeter_triangulation;
+
+template <class Function>
+void test_function(const Function& fun) {
+  Eigen::VectorXd seed = fun.seed();
+  Eigen::VectorXd res_seed = fun(fun.seed());
+  BOOST_CHECK( seed.size() == (long int)fun.amb_d() );
+  BOOST_CHECK( res_seed.size() == (long int)fun.cod_d() );
+  for (std::size_t i = 0; i < fun.cod_d(); i++)
+    GUDHI_TEST_FLOAT_EQUALITY_CHECK(res_seed(i), 0., 1e-10);
+}
+
+BOOST_AUTO_TEST_CASE(function) {
+
+  {
+    // the sphere testing part
+    std::size_t m = 3, d = 5;
+    Eigen::VectorXd center(d); center << 2, 1.5, -0.5, 4.5, -1;
+    double radius = 5;
+    typedef Function_Sm_in_Rd Function_sphere;
+    Function_sphere fun_sphere(radius, m, d, center);
+    test_function(fun_sphere);
+  }
+  {
+    // the affine plane testing part
+    std::size_t m = 0, d = 5;
+    Eigen::MatrixXd normal_matrix = Eigen::MatrixXd::Zero(d, d-m);
+    for (std::size_t i = 0; i < d-m; ++i)
+      normal_matrix(i,i) = 1;
+    typedef Function_affine_plane_in_Rd Function_plane;
+    Function_plane fun_plane(normal_matrix);
+    test_function(fun_plane);
+  }
+  {
+    // the constant function testing part
+    std::size_t k = 2, d = 5;
+    auto x = Eigen::VectorXd::Constant(k, 1);
+    Constant_function fun_const(d, k, x);
+    Eigen::VectorXd res_zero = fun_const(Eigen::VectorXd::Zero(d));
+    for (std::size_t i = 0; i < k; ++i)
+      GUDHI_TEST_FLOAT_EQUALITY_CHECK(res_zero(i), x(i), 1e-10);
+  }
+  {
+    // the chair function
+    Function_chair_in_R3 fun_chair;
+    test_function(fun_chair);
+  }
+  {
+    // the torus function
+    Function_torus_in_R3 fun_torus;
+    test_function(fun_torus);
+  }
+  {
+    // the whitney umbrella function
+    Function_whitney_umbrella_in_R3 fun_umbrella;
+    test_function(fun_umbrella);
+  }
+  {
+    // the lemniscate revolution function
+    Function_lemniscate_revolution_in_R3 fun_lemniscate;
+    test_function(fun_lemniscate);
+  }
+  {
+    // the iron function
+    Function_iron_in_R3 fun_iron;
+    test_function(fun_iron);
+  }
+  {
+    Function_moment_curve_in_Rd fun_moment_curve(3, 5);
+    test_function(fun_moment_curve);
+  }
+  {
+    // random orthogonal matrix
+    Eigen::MatrixXd matrix = random_orthogonal_matrix(5);
+    Eigen::MatrixXd id_matrix = matrix.transpose() * matrix;
+    for (std::size_t i = 0; i < 5; ++i)
+      for (std::size_t j = 0; j < 5; ++j)
+	if (i == j)
+	  GUDHI_TEST_FLOAT_EQUALITY_CHECK(id_matrix(i,j), 1.0, 1e-10);
+	else
+	  GUDHI_TEST_FLOAT_EQUALITY_CHECK(id_matrix(i,j), 0.0, 1e-10);
+  }
+  {
+    // function embedding
+    Function_iron_in_R3 fun_iron;
+    auto fun_embed = make_embedding(fun_iron, 5);
+    test_function(fun_iron);
+
+    // function translation
+    Eigen::VectorXd off = Eigen::VectorXd::Random(5);
+    auto fun_trans = translate(fun_embed, off);
+    test_function(fun_trans);
+
+    // function linear transformation
+    Eigen::MatrixXd matrix = Eigen::MatrixXd::Random(5, 5);
+    BOOST_CHECK( matrix.determinant() != 0. );
+    auto fun_lin = make_linear_transformation(fun_trans, matrix);
+    test_function(fun_lin);
+
+    // function negative
+    auto fun_neg = negation(fun_lin);
+    test_function(fun_neg);
+
+    // function product
+    typedef Function_Sm_in_Rd Function_sphere;
+    Function_sphere fun_sphere(1, 1);
+    auto fun_prod = make_product_function(fun_sphere, fun_sphere, fun_sphere);
+    test_function(fun_prod);
+    
+    // function PL approximation
+    Coxeter_triangulation<> cox_tr(6);
+    typedef Coxeter_triangulation<>::Vertex_handle Vertex_handle;
+    auto fun_pl = make_pl_approximation(fun_prod, cox_tr);
+    Vertex_handle v0 = Vertex_handle(cox_tr.dimension(), 0);
+    Eigen::VectorXd x0 = cox_tr.cartesian_coordinates(v0);
+    Eigen::VectorXd value0 = fun_prod(x0);
+    Eigen::VectorXd pl_value0 = fun_pl(x0);
+    for (std::size_t i = 0; i < fun_pl.cod_d(); i++)
+      GUDHI_TEST_FLOAT_EQUALITY_CHECK(value0(i), pl_value0(i), 1e-10);
+    Vertex_handle v1 = v0;
+    v1[0] += 1;
+    Eigen::VectorXd x1 = cox_tr.cartesian_coordinates(v1);
+    Eigen::VectorXd value1 = fun_prod(x1);
+    Eigen::VectorXd pl_value1 = fun_pl(x1);
+    for (std::size_t i = 0; i < fun_pl.cod_d(); i++)
+      GUDHI_TEST_FLOAT_EQUALITY_CHECK(value1(i), pl_value1(i), 1e-10);
+    Eigen::VectorXd pl_value_mid = fun_pl(0.5*x0 + 0.5*x1);
+    for (std::size_t i = 0; i < fun_pl.cod_d(); i++)
+      GUDHI_TEST_FLOAT_EQUALITY_CHECK(0.5*value0(i) + 0.5*value1(i), pl_value_mid(i), 1e-10);
+  }
+}
diff --git a/src/Coxeter_triangulation/test/manifold_tracing_test.cpp b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
new file mode 100644
index 00000000..cff6ddfa
--- /dev/null
+++ b/src/Coxeter_triangulation/test/manifold_tracing_test.cpp
@@ -0,0 +1,66 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "manifold_tracing"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#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>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(manifold_tracing) {
+  // manifold without boundary
+  Function_Sm_in_Rd fun_sph(5.1111, 2);
+  auto oracle = make_oracle(fun_sph);
+  Coxeter_triangulation<> cox_tr(oracle.amb_d());
+  // cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+  
+  using MT = Manifold_tracing<Coxeter_triangulation<> >;
+  Eigen::VectorXd seed = fun_sph.seed();
+  std::vector<Eigen::VectorXd> seed_points(1, seed);
+  typename MT::Out_simplex_map out_simplex_map;
+  manifold_tracing_algorithm(seed_points, cox_tr, oracle, out_simplex_map);
+
+  for (auto si_pair: out_simplex_map) {
+    BOOST_CHECK ( si_pair.first.dimension() == oracle.function().cod_d() );
+    BOOST_CHECK ( si_pair.second.size() == (long int)oracle.function().amb_d() );
+  }
+  std::cout << "out_simplex_map.size() = " << out_simplex_map.size() << "\n";
+  BOOST_CHECK ( out_simplex_map.size() == 1140 );
+  
+
+  // manifold with boundary
+  Function_Sm_in_Rd fun_boundary(3.0, 2, fun_sph.seed());
+  auto oracle_with_boundary = make_oracle(fun_sph, fun_boundary);
+  typename MT::Out_simplex_map interior_simplex_map, boundary_simplex_map;
+  manifold_tracing_algorithm(seed_points, cox_tr, oracle_with_boundary, interior_simplex_map, boundary_simplex_map);
+  for (auto si_pair: interior_simplex_map) {
+    BOOST_CHECK ( si_pair.first.dimension() == oracle.function().cod_d() );
+    BOOST_CHECK ( si_pair.second.size() == (long int)oracle.function().amb_d() );
+  }
+  std::cout << "interior_simplex_map.size() = " << interior_simplex_map.size() << "\n";
+  BOOST_CHECK ( interior_simplex_map.size() == 96 );
+  for (auto si_pair: boundary_simplex_map) {
+    BOOST_CHECK ( si_pair.first.dimension() == oracle.function().cod_d()+1 );
+    BOOST_CHECK ( si_pair.second.size() == (long int)oracle.function().amb_d() );
+  }
+  std::cout << "boundary_simplex_map.size() = " << boundary_simplex_map.size() << "\n";
+  BOOST_CHECK ( boundary_simplex_map.size() == 55 );
+  
+}
diff --git a/src/Coxeter_triangulation/test/oracle_test.cpp b/src/Coxeter_triangulation/test/oracle_test.cpp
new file mode 100644
index 00000000..8d371656
--- /dev/null
+++ b/src/Coxeter_triangulation/test/oracle_test.cpp
@@ -0,0 +1,60 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "oracle"
+#include <boost/test/unit_test.hpp>
+#include <gudhi/Unitary_tests_utils.h>
+
+#include <string>
+
+#include <gudhi/Implicit_manifold_intersection_oracle.h>
+
+#include <gudhi/Functions/Function_Sm_in_Rd.h>
+#include <gudhi/Functions/Cartesian_product.h>
+
+#include <gudhi/Coxeter_triangulation.h>
+
+#include <random>
+#include <cstdlib>
+
+using namespace Gudhi::coxeter_triangulation;
+
+BOOST_AUTO_TEST_CASE(oracle) {
+
+  Function_Sm_in_Rd fun_sph(5.1111, 2);
+  auto oracle = make_oracle(fun_sph);
+  Coxeter_triangulation<> cox_tr(oracle.amb_d());
+  // cox_tr.change_offset(Eigen::VectorXd::Random(oracle.amb_d()));
+  
+  Eigen::VectorXd seed = fun_sph.seed();
+  auto s = cox_tr.locate_point(seed);
+
+  std::size_t num_intersected_edges = 0;
+  for (auto f: s.face_range(oracle.cod_d())) {
+    auto qr = oracle.intersects(f, cox_tr);
+    if (qr.success)
+      num_intersected_edges++;
+    auto vertex_it = f.vertex_range().begin();
+    Eigen::Vector3d p1 = cox_tr.cartesian_coordinates(*vertex_it++);
+    Eigen::Vector3d p2 = cox_tr.cartesian_coordinates(*vertex_it++);
+    BOOST_CHECK( vertex_it == f.vertex_range().end() );
+    Eigen::MatrixXd m(3,3);
+    if (qr.success) {      
+      m.col(0) = qr.intersection;
+      m.col(1) = p1;
+      m.col(2) = p2;
+      GUDHI_TEST_FLOAT_EQUALITY_CHECK(m.determinant(), 0.0, 1e-10);
+    }
+  }
+  BOOST_CHECK( num_intersected_edges == 3 || num_intersected_edges == 4 );
+}
diff --git a/src/Coxeter_triangulation/test/perm_rep_test.cpp b/src/Coxeter_triangulation/test/perm_rep_test.cpp
new file mode 100644
index 00000000..f2cdbce1
--- /dev/null
+++ b/src/Coxeter_triangulation/test/perm_rep_test.cpp
@@ -0,0 +1,66 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Siargey Kachanovich
+ *
+ *    Copyright (C) 2019 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "permutahedral_representation"
+#include <boost/test/unit_test.hpp>
+
+#include <gudhi/Permutahedral_representation.h>
+
+BOOST_AUTO_TEST_CASE(permutahedral_representation) {
+  typedef std::vector<int> Vertex;
+  typedef std::vector<std::size_t> Part;
+  typedef std::vector<Part> Partition;
+  typedef Gudhi::coxeter_triangulation::Permutahedral_representation<Vertex, Partition> Simplex_handle;
+  Vertex v0(10, 0);
+  Partition omega = {Part({5}), Part({2}), Part({3,7}), Part({4,9}), Part({0,6,8}), Part({1,10})};
+  Simplex_handle s(v0, omega);
+
+  // Dimension check
+  BOOST_CHECK(s.dimension() == 5);
+
+  // Vertex number check
+  std::vector<Vertex> vertices;
+  for (auto& v: s.vertex_range())
+    vertices.push_back(v);
+  BOOST_CHECK(vertices.size() == 6);
+  
+  // Facet number check
+  std::vector<Simplex_handle> facets;
+  for (auto& f: s.facet_range())
+    facets.push_back(f);
+  BOOST_CHECK(facets.size() == 6);
+
+  // Face of dim 3 number check
+  std::vector<Simplex_handle> faces3;
+  for (auto& f: s.face_range(3))
+    faces3.push_back(f);
+  BOOST_CHECK(faces3.size() == 15);
+
+  // Cofacet number check
+  std::vector<Simplex_handle> cofacets;
+  for (auto& f: s.cofacet_range())
+    cofacets.push_back(f);
+  BOOST_CHECK(cofacets.size() == 12);
+
+  // Is face check
+  Vertex v1(10, 0);
+  Partition omega1 = {Part({5}), Part({0,1,2,3,4,6,7,8,9,10})};
+  Simplex_handle s1(v1, omega1);
+  Vertex v2(10, 0); v2[1] = -1;
+  Partition omega2 = {Part({1}), Part({5}), Part({2}), Part({3,7}), Part({4,9}), Part({0,6,8}), Part({10})};
+  Simplex_handle s2(v2, omega2);
+  BOOST_CHECK(s.is_face_of(s));
+  BOOST_CHECK(s1.is_face_of(s));
+  BOOST_CHECK(!s2.is_face_of(s));
+  BOOST_CHECK(s.is_face_of(s2));
+}