Adding arbitrary-dimension simplex shape (#22)

* Adding NSimplex shape

docs/doxygen/ShivaConfig.hpp

@@ -8,7 +8,7 @@
 
 #define SHIVA_VERSION_PATCHLEVEL 0
 
-#define SHIVA_USE_CUDA
+/* #undef SHIVA_USE_CUDA */
 
 /* #undef SHIVA_USE_HIP */
 

src/common/MathUtilities.hpp

@@ -62,7 +62,23 @@ struct factorial< T, N, std::integer_sequence< T, ISEQ... > >
   static constexpr T value = (static_cast< T >(1) * ... *(ISEQ + 1));
 };
 
+/**
+ * @brief Struct that computes the binomial coefficient N over K
+ * @tparam T The type of the number
+ * @tparam N The number of objects in the set
+ * @tparam K The number of objects to be chose (irrespective of the order) from the set
+ */
+template< typename T, T N, T K >
+struct binomialCoefficient
+{
+  static_assert( 0 <= K, "K must be greater or equal than 0" );
+  static_assert( K <= N, "N must be greater or equal than K" );
 
+  /**
+   * @brief The binomial coefficient N over K
+   */
+  static constexpr T value = factorial< T, N >::value / ( factorial< T, K >::value * factorial< T, N - K >::value );
+};
 
 /**
  * @brief Inverse of a 3x3 matrix

src/geometry/CMakeLists.txt

@@ -6,6 +6,7 @@ set( geometry_headers
      mapping/UniformScaling.hpp
      shapes/InterpolatedShape.hpp
      shapes/NCube.hpp
+     shapes/NSimplex.hpp
      )
 
 set( geometry_sources

src/geometry/shapes/NCube.hpp

@@ -107,6 +107,22 @@ class NCube
     }
   }
 
+  /**
+   * @brief Returns the number of cells (3-cube) in the n-cube.
+   * @return The number of cells (3-cube) in the n-cube.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numCells()
+  {
+    if constexpr ( N > 2 )
+    {
+      return numMCubes< 3 >();
+    }
+    else
+    {
+      return 0;
+    }
+  }
+
   /**
    * @brief Returns the number of hyperfaces (n-1-cube) in the n-cube.
    * @return The number of hyperfaces (n-1-cube) in the n-cube.

src/geometry/shapes/NSimplex.hpp

@@ -0,0 +1,214 @@
+/**
+ * @file NSimplex.hpp
+ */
+
+#pragma once
+
+#include "common/ShivaMacros.hpp"
+#include "common/IndexTypes.hpp"
+#include "common/types.hpp"
+#include "common/MathUtilities.hpp"
+
+namespace shiva
+{
+namespace geometry
+{
+
+/**
+ * @brief NSimplex represents a n-simplex.
+ * @tparam REAL_TYPE The floating point type.
+ * @tparam N The number of dimensions of the n-simplex.
+ * @tparam MIN The minimum coordinate of the n-simplex.
+ * @tparam MAX The maximum coordinate of the n-simplex.
+ * @tparam DIVISOR The divisor of the coordinates of the n-simplex. This is
+ * required because the coordinates of the n-simplex are integers, but the
+ * coordinates of the n-simplex are floating point numbers.
+ *
+ * A n-simplex is a generalization of a triangle (n = 2) or tetrahedron (n = 3)
+ * to arbitrary dimensions
+ * <a href="https://en.wikipedia.org/wiki/Simplex"> Wikipedia Simplex</a>
+ */
+template< typename REAL_TYPE, int N, int MIN, int MAX, int DIVISOR >
+class NSimplex
+{
+public:
+
+  static_assert( MIN < MAX, "MIN must be less than MAX" );
+  static_assert( DIVISOR > 0, "DIVISOR must be greater than 0" );
+
+  /**
+   * @brief The number of dimension of the simplex.
+   * @return The number dimension of the simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numDims() {return N;}
+
+  /// Alias for the floating point type.
+  using RealType = REAL_TYPE;
+
+  /// Alias for the floating point type for the coordinates.
+  using CoordType = REAL_TYPE[N];
+
+  /**
+   * @brief Returns the number of m-simplexes in the n-simplex.
+   * @tparam M The number of dimensions of the m-simplex.
+   * @return The number of m-simplexes in the n-simplex. An m-simplex is the lower
+   * dimensional object contained in the n-simplex. For instance, the 0-simplex is a
+   * vertex, the 1-simplex is a line, the 2-simplex is a triangle, the 3-simplex is
+   * a tetrahedron, etc. M must be less than or equal to N
+   */
+  template< int M >
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numMSimplexes()
+  {
+    static_assert( M <= N, "M must be less than or equal to N" );
+    return mathUtilities::binomialCoefficient< int, N + 1, M + 1 >::value;
+  }
+
+  /**
+   * @brief Returns the number of vertices (0-simplex) in the n-simplex.
+   * @return The number of vertices (0-simplex) in the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numVertices()
+  {
+    return numMSimplexes< 0 >();
+  }
+
+  /**
+   * @brief Returns the number of edges (1-simplex) in the n-simplex.
+   * @return The number of edges (1-simplex) in the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numEdges()
+  {
+    if constexpr ( N > 0 )
+    {
+      return numMSimplexes< 1 >();
+    }
+    else
+    {
+      return 0;
+    }
+  }
+
+  /**
+   * @brief Returns the number of faces (2-simplex) in the n-simplex.
+   * @return The number of faces (2-simplex) in the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numFaces()
+  {
+    if constexpr ( N > 1 )
+    {
+      return numMSimplexes< 2 >();
+    }
+    else
+    {
+      return 0;
+    }
+  }
+
+  /**
+   * @brief Returns the number of cells (3-simplex) in the n-simplex.
+   * @return The number of cells (3-simplex) in the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numCells()
+  {
+    if constexpr ( N > 2 )
+    {
+      return numMSimplexes< 3 >();
+    }
+    else
+    {
+      return 0;
+    }
+  }
+
+  /**
+   * @brief Returns the number of hyperfaces (n-1-simplex) in the n-simplex.
+   * @return The number of hyperfaces (n-1-simplex) in the n-simplex.
+   *
+   * The hyperfaces can be considered the number of n-1 dimensional objects in
+   * an n-simplex. For instance, the hyperfaces of a tetrahedron are the faces of
+   * the tetrahedron. The hyperfaces of a triangle are the edges of the triangle.
+   * The hyperfaces of a line are the vertices of the line.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE int numHyperFaces()
+  {
+    return numMSimplexes< N - 1 >();
+  }
+
+  /**
+   * @brief returns the minimum coordinate of the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE REAL_TYPE minCoord()
+  {
+    return static_cast< RealType >(MIN) / DIVISOR;
+  }
+
+  /**
+   * @brief returns the maximum coordinate of the n-simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE REAL_TYPE maxCoord()
+  {
+    return static_cast< RealType >(MAX) / DIVISOR;
+  }
+
+
+  /**
+   * @brief Returns the length dimension of the simplex.
+   * @return The length of the simplex.
+   */
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE REAL_TYPE length()
+  {
+    return static_cast< RealType >( MAX - MIN ) / DIVISOR;
+  }
+
+
+  /**
+   * @brief Returns the volume of the n-simplex.
+   * @return The volume of the n-simplex.
+   */
+  template< int DIM = N >
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE
+  typename std::enable_if< ( DIM > 0 ), REAL_TYPE >::type
+  volume()
+  {
+    return volumeHelper();
+  }
+
+
+private:
+  /**
+   * @brief Helper recusive function for volume calculation.
+   * @tparam DIM The current dimension for the recursion.
+   * @return The volume of the n-simplex.
+   */
+  template< int DIM = N >
+  SHIVA_STATIC_CONSTEXPR_HOSTDEVICE_FORCEINLINE REAL_TYPE
+  volumeHelper()
+  {
+    if constexpr ( DIM == 0 )
+    {
+      return 1;
+    }
+    else
+    {
+      return length() / DIM * volumeHelper< DIM - 1 >();
+    }
+  }
+
+};
+
+/**
+ * @brief Alias for a 2-simplex.
+ * @tparam REAL_TYPE The floating point type.
+ */
+template< typename REAL_TYPE >
+using Triangle = NSimplex< REAL_TYPE, 2, 0, 1, 1 >;
+
+/**
+ * @brief Alias for a 3-simplex
+ * @tparam REAL_TYPE The floating point type.
+ */
+template< typename REAL_TYPE >
+using Tetrahedron = NSimplex< REAL_TYPE, 3, 0, 1, 1 >;
+
+} // namespace geometry
+} // namespace shiva

src/geometry/shapes/unitTests/CMakeLists.txt

@@ -4,6 +4,7 @@
 
 set( unit_tests_sources
      testNCube.cpp
+     testNSimplex.cpp
      testInterpolatedShape.cpp
      )
 

src/geometry/shapes/unitTests/testNCube.cpp

@@ -10,7 +10,7 @@ using namespace shiva::geometry;
 template< typename NCUBE >
 struct NCubeSolutions;
 
-#define NCUBE( NCUBE_NAME, N, MIN, MAX, DIVISOR, NUM_VERTICES, NUM_EDGES, NUM_FACES, NUM_HYPERFACES, MIN_COORD, MAX_COORD, LENGTH, VOLUME ) \
+#define NCUBE( NCUBE_NAME, N, MIN, MAX, DIVISOR, NUM_VERTICES, NUM_EDGES, NUM_FACES, NUM_CELLS, NUM_HYPERFACES, MIN_COORD, MAX_COORD, LENGTH, VOLUME ) \
         using NCUBE_NAME = NCube< double, N, MIN, MAX, DIVISOR >; \
         template<> \
         struct NCubeSolutions< NCUBE_NAME > \
@@ -19,6 +19,7 @@ struct NCubeSolutions;
           static constexpr int numVertices() { return NUM_VERTICES;} \
           static constexpr int numEdges() { return NUM_EDGES;} \
           static constexpr int numFaces() { return NUM_FACES;} \
+          static constexpr int numCells() { return NUM_CELLS;} \
           static constexpr int numHyperFaces() { return NUM_HYPERFACES;} \
           static constexpr double minCoord() { return MIN_COORD;} \
           static constexpr double maxCoord() { return MAX_COORD;} \
@@ -29,57 +30,57 @@ struct NCubeSolutions;
 //*****************************************************************************
 NCUBE( NCube_1_0_1_1,
        1, 0, 1, 1,
-       2, 1, 0, 2, 0.0, 1.0, 1.0, 1.0 )
+       2, 1, 0, 0, 2, 0.0, 1.0, 1.0, 1.0 )
 
 NCUBE( NCube_1_m1_1_2,
        1, -1, 1, 2,
-       2, 1, 0, 2, -0.5, 0.5, 1.0, 1.0 )
+       2, 1, 0, 0, 2, -0.5, 0.5, 1.0, 1.0 )
 
 NCUBE( NCube_1_m1_1_1,
        1, -1, 1, 1,
-       2, 1, 0, 2, -1.0, 1.0, 2.0, 2.0 )
+       2, 1, 0, 0, 2, -1.0, 1.0, 2.0, 2.0 )
 
 
 //*****************************************************************************
 NCUBE( NCube_2_0_1_1,
        2, 0, 1, 1,
-       4, 4, 1, 4, 0.0, 1.0, 1.0, 1.0 )
+       4, 4, 1, 0, 4, 0.0, 1.0, 1.0, 1.0 )
 
 NCUBE( NCube_2_m1_1_2,
        2, -1, 1, 2,
-       4, 4, 1, 4, -0.5, 0.5, 1.0, 1.0 )
+       4, 4, 1, 0, 4, -0.5, 0.5, 1.0, 1.0 )
 
 NCUBE( NCube_2_m1_1_1,
        2, -1, 1, 1,
-       4, 4, 1, 4, -1.0, 1.0, 2.0, 4.0 )
+       4, 4, 1, 0, 4, -1.0, 1.0, 2.0, 4.0 )
 
 
 //*****************************************************************************
 NCUBE( NCube_3_0_1_1,
        3, 0, 1, 1,
-       8, 12, 6, 6, 0.0, 1.0, 1.0, 1.0 )
+       8, 12, 6, 1, 6, 0.0, 1.0, 1.0, 1.0 )
 
 NCUBE( NCube_3_m1_1_2,
        3, -1, 1, 2,
-       8, 12, 6, 6, -0.5, 0.5, 1.0, 1.0 )
+       8, 12, 6, 1, 6, -0.5, 0.5, 1.0, 1.0 )
 
 NCUBE( NCube_3_m1_1_1,
        3, -1, 1, 1,
-       8, 12, 6, 6, -1.0, 1.0, 2.0, 8.0 )
+       8, 12, 6, 1, 6, -1.0, 1.0, 2.0, 8.0 )
 
 
 //*****************************************************************************
 NCUBE( NCube_4_0_1_1,
        4, 0, 1, 1,
-       16, 32, 24, 8, 0.0, 1.0, 1.0, 1.0 )
+       16, 32, 24, 8, 8, 0.0, 1.0, 1.0, 1.0 )
 
 NCUBE( NCube_4_m1_1_2,
        4, -1, 1, 2,
-       16, 32, 24, 8, -0.5, 0.5, 1.0, 1.0 )
+       16, 32, 24, 8, 8, -0.5, 0.5, 1.0, 1.0 )
 
 NCUBE( NCube_4_m1_1_1,
        4, -1, 1, 1,
-       16, 32, 24, 8, -1.0, 1.0, 2.0, 16.0 )
+       16, 32, 24, 8, 8, -1.0, 1.0, 2.0, 16.0 )
 
 
 template< typename NCUBE >
@@ -90,6 +91,7 @@ testNCubeHelper()
   static_assert( NCUBE::numVertices() == NCubeSolutions< NCUBE >::numVertices() );
   static_assert( NCUBE::numEdges() == NCubeSolutions< NCUBE >::numEdges() );
   static_assert( NCUBE::numFaces() == NCubeSolutions< NCUBE >::numFaces() );
+  static_assert( NCUBE::numCells() == NCubeSolutions< NCUBE >::numCells() );
   static_assert( NCUBE::numHyperFaces() == NCubeSolutions< NCUBE >::numHyperFaces() );
   static_assert( pmpl::check( NCUBE::minCoord(), NCubeSolutions< NCUBE >::minCoord(), 1.0e-12 ) );
   static_assert( pmpl::check( NCUBE::maxCoord(), NCubeSolutions< NCUBE >::maxCoord(), 1.0e-12 ) );