Add triangles with pseudonormals and faces.

.github/workflows/CI.yml

@@ -85,8 +85,8 @@ jobs:
         working-directory: ${{ matrix.directory }}
         run: |
           source /opt/intel/oneapi/setvars.sh
-          icpc -std=c++17 -O0 \
-          -Wall -Warray-bounds -Wchar-subscripts -Wcomment -Wenum-compare -Wformat -Wuninitialized -Wmaybe-uninitialized -Wmain -Wnarrowing -Wnonnull \
+          icpx -std=c++17 -O0 \
+          -Wall -Warray-bounds -Wcomment -Wenum-compare -Wformat -Wuninitialized -Wmaybe-uninitialized -Wmain -Wnarrowing -Wnonnull \
           -Wparentheses -Wpointer-sign -Wreorder -Wreturn-type -Wsign-compare -Wsequence-point -Wtrigraphs -Wunused-function -Wunused-but-set-variable \
           -Wunused-variable -Wwrite-strings -Werror \
           main.cpp -lstdc++fs

Docs/Sphinx/source/Concepts.rst

@@ -155,6 +155,13 @@ where the sum runs over all faces which share :math:`v` as a vertex, and where :
 
    Edge and vertex pseudonormals.
 
+Triangles
+=========
+
+``EBGeometry`` also supports using direct triangles rather than DCEL polygons.
+In this case the user can parse his/her files and ask for triangle description rather than a DCEL description.
+Internally, triangles are represented as *pseudo-triangles* that contain normal vectors on both the edges and vertices, completely similar to a DCEL polygon.
+
 .. _Chap:BVH:
 
 Bounding volume hierarchies

Docs/Sphinx/source/ImplemBVH.rst

@@ -85,6 +85,7 @@ To use bottom-up construction, one may use the member function
 .. literalinclude:: ../../../Source/EBGeometry_BVH.hpp
    :language: c++
    :lines: 298-309
+   :dedent: 4
 
 The template argument is the space-filling curve that the user wants to apply.
 Currently, we support Morton codes and nested indices.
@@ -273,6 +274,7 @@ If the traversal decides to visit a node, it immediately computes the specified
 .. literalinclude:: ../../../Source/EBGeometry_BVHImplem.hpp
    :language: c++
    :lines: 284-322
+   :dedent: 2
    :caption: Tree traversal algorithm for the BVH tree.
 
 Traversal examples
@@ -294,7 +296,7 @@ These rules are given below.
 
 .. literalinclude:: ../../../Source/EBGeometry_MeshDistanceFunctionsImplem.hpp
    :language: c++
-   :lines: 97-132
+   :lines: 145-181
    :caption: Tree traversal criterion for computing the signed distance to a DCEL mesh using the BVH accelerator.
 	     See :file:`Source/EBGeometry_MeshDistanceFunctionsImplem.hpp` for details.
 

Docs/Sphinx/source/Parsers.rst

@@ -25,6 +25,18 @@ If you have one or multiple STL files, you can quickly turn them into signed dis
 This will build DCEL meshes for each input file, and wrap the meshes in BVHs.
 See :ref:`Chap:LinearSTL` for further details.
 
+.. tip::
+
+   If the input files consist only of triangles, use the version
+
+   .. code-block:: c++
+
+      std::vector<std::string> files; // <---- List of file names.
+
+      const auto distanceFields = EBGeometry::Parser::readIntoTriangleBVH<float>(files);
+
+   This version will convert all DCEL polygons to triangles, and usually provides a nice code speedup.
+
 Reading STL files
 -----------------
 
@@ -42,7 +54,8 @@ To read one or multiple STL files and turn it into DCEL meshes, use
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 53-67
+   :lines: 54-68
+   :dedent: 2	   
 
 Note that this will only expose the DCEL mesh, but not include any signed distance functionality.
 
@@ -53,7 +66,8 @@ To read one or multiple STL files and also turn it into signed distance represen
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 69-83
+   :lines: 70-84
+   :dedent: 2	   
 
 DCEL mesh SDF with full BVH
 ___________________________
@@ -62,7 +76,8 @@ To read one or multiple STL files and turn it into signed distance representatio
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 85-105
+   :lines: 86-106
+   :dedent: 2	   
 
 .. _Chap:LinearSTL:
 
@@ -73,7 +88,22 @@ To read one or multiple STL files and turn it into signed distance representatio
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 107-127
+   :lines: 107-128
+   :dedent: 2
+
+
+Triangle meshes with BVH
+________________________
+
+To read one or multiple STL files and turn it into signed distance representations using a compact BVH representation, use
+
+.. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
+   :language: c++
+   :lines: 130-147
+   :dedent: 2	   
+
+This version differs from the DCEL meshes in that each DCEL polygon is converted into triangles after parsing.
+The code will throw an error if not all DCEL polygon are not actual triangles.
 
 .. _Chap:PolySoups:
 
@@ -94,7 +124,8 @@ To turn this into a DCEL mesh, one should compress the triangle soup (get rid of
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 146-165
+   :lines: 182-201
+   :dedent: 2
 
 The ``compress`` function will discard duplicate vertices from the soup, while the ``soupToDCEL`` will tie the remaining polygons into a DCEL mesh.
 This function will also compute the vertex and edge normal vectors.

EBGeometry.hpp

@@ -14,6 +14,8 @@
 #include "Source/EBGeometry_SignedDistanceFunction.hpp"
 #include "Source/EBGeometry_SimpleTimer.hpp"
 #include "Source/EBGeometry_Transform.hpp"
+#include "Source/EBGeometry_Triangle.hpp"
+#include "Source/EBGeometry_Triangle2D.hpp"
 
 /*!
   @brief Name space for all of EBGeometry

Examples/AMReX_DCEL/main.cpp

@@ -31,7 +31,7 @@ class AMReXSDF
   AMReXSDF(const std::string a_filename, const bool a_use_bvh)
   {
     if (a_use_bvh) {
-      m_sdf = EBGeometry::Parser::readIntoLinearBVH<T, Meta, BV, K>(a_filename);
+      m_sdf = EBGeometry::Parser::readIntoTriangleBVH<T, Meta, BV, K>(a_filename);
     }
     else {
       m_sdf = EBGeometry::Parser::readIntoMesh<T, Meta>(a_filename);

Examples/AMReX_F18/main.cpp

@@ -25,7 +25,7 @@ using Vec3         = EBGeometry::Vec3T<T>;
 using BV           = EBGeometry::BoundingVolumes::AABBT<T>;
 using Face         = EBGeometry::DCEL::FaceT<T, MetaData>;
 using Mesh         = EBGeometry::DCEL::MeshT<T, MetaData>;
-using FastSDF      = EBGeometry::FastCompactMeshSDF<T, MetaData, BV, K>;
+using FastSDF      = EBGeometry::FastTriMeshSDF<T, MetaData, BV, K>;
 
 // F18 geometry, using nifty EBGeometry bindings and accelerators.
 class F18

Examples/AMReX_PaintEB/main.cpp

@@ -37,13 +37,14 @@ class AMReXSDF
   {
     // Read in the mesh into a DCEL mesh and partition it into a bounding volume hierarchy
     auto mesh = EBGeometry::Parser::readIntoDCEL<T, Meta>(a_filename);
-    m_sdf     = std::make_shared<EBGeometry::FastCompactMeshSDF<T, Meta, BV, K>>(mesh);
 
     // Set the meta-data for all facets to their "index", i.e. position in the list of facets
     auto& faces = mesh->getFaces();
     for (size_t i = 0; i < faces.size(); i++) {
       faces[i]->getMetaData() = 1.0 * i;
     }
+
+    m_sdf = std::make_shared<EBGeometry::FastTriMeshSDF<T, Meta, BV, K>>(mesh);
   }
 
   /*!

Examples/Chombo3_DCEL/main.cpp

@@ -27,7 +27,7 @@ class ChomboSDF : public BaseIF
 
   ChomboSDF(const std::string a_filename)
   {
-    m_implicitFunction = EBGeometry::Parser::readIntoLinearBVH<T, Meta, BV, K>(a_filename);
+    m_implicitFunction = EBGeometry::Parser::readIntoTriangleBVH<T, Meta, BV, K>(a_filename);
     m_implicitFunction = EBGeometry::Complement<T>(m_implicitFunction);
   }
 

Examples/Chombo3_F18/main.cpp

@@ -22,7 +22,7 @@ using MetaData     = EBGeometry::DCEL::DefaultMetaData;
 using Face         = EBGeometry::DCEL::FaceT<T>;
 using Mesh         = EBGeometry::DCEL::MeshT<T>;
 using BV           = EBGeometry::BoundingVolumes::AABBT<T>;
-using FastSDF      = EBGeometry::FastCompactMeshSDF<T, MetaData, BV, K>;
+using FastSDF      = EBGeometry::FastTriMeshBVH<T, MetaData, BV, K>;
 using Vec3         = EBGeometry::Vec3T<T>;
 
 class F18 : public BaseIF

Examples/EBGeometry_DCEL/main.cpp

@@ -21,6 +21,7 @@ using T    = float;
 using Meta = short;
 using BV   = EBGeometry::BoundingVolumes::AABBT<T>;
 using Vec3 = EBGeometry::Vec3T<T>;
+using Tri  = EBGeometry::Triangle<T, Meta>;
 
 int
 main(int argc, char* argv[])
@@ -48,6 +49,7 @@ main(int argc, char* argv[])
   const auto dcelSDF = EBGeometry::Parser::readIntoMesh<T, Meta>(file);
   const auto bvhSDF  = EBGeometry::Parser::readIntoFullBVH<T, Meta, BV, K>(file);
   const auto linSDF  = EBGeometry::Parser::readIntoLinearBVH<T, Meta, BV, K>(file);
+  const auto triSDF  = EBGeometry::Parser::readIntoTriangleBVH<T, Meta, BV, K>(file);
 
   // Sample some random points around the object.
   constexpr size_t Nsamp = 100;
@@ -74,6 +76,7 @@ main(int argc, char* argv[])
   T dcelSum = 0.0;
   T bvhSum  = 0.0;
   T linSum  = 0.0;
+  T triSum  = 0.0;
 
   const auto t0 = std::chrono::high_resolution_clock::now();
   for (const auto& x : ranPoints) {
@@ -88,24 +91,33 @@ main(int argc, char* argv[])
     linSum += linSDF->signedDistance(x);
   }
   const auto t3 = std::chrono::high_resolution_clock::now();
+  for (const auto& x : ranPoints) {
+    triSum += triSDF->signedDistance(x);
+  }
+  const auto t4 = std::chrono::high_resolution_clock::now();
 
   const std::chrono::duration<T, std::micro> dcelTime = t1 - t0;
   const std::chrono::duration<T, std::micro> bvhTime  = t2 - t1;
   const std::chrono::duration<T, std::micro> linTime  = t3 - t2;
+  const std::chrono::duration<T, std::micro> triTime  = t4 - t3;
 
   if (std::abs(bvhSum - dcelSum) > std::numeric_limits<T>::epsilon()) {
     std::cerr << "Full BVH did not give same distance! Diff = " << bvhSum - dcelSum << "\n";
   }
   if (std::abs(linSum - dcelSum) > std::numeric_limits<T>::epsilon()) {
     std::cerr << "Linearized BVH did not give same distance! Diff = " << linSum - dcelSum << "\n";
   }
+  if (std::abs(triSum - dcelSum) > std::numeric_limits<T>::epsilon()) {
+    std::cerr << "TriMesh BVH did not give same distance! Diff = " << triSum - dcelSum << "\n";
+  }
 
   // clang-format off
   std::cout << "Bounding box = " << lo << "\t" << hi << "\n";
   std::cout << "Accumulated distance and time using direct DCEL = " << dcelSum << ", which took " << dcelTime.count() / Nsamp << " us\n";
   std::cout << "Accumulated distance and time using full BVH    = " << bvhSum  << ", which took " << bvhTime.count()  / Nsamp << " us\n";
   std::cout << "Accumulated distance and time using compact BVH = " << linSum  << ", which took " << linTime.count()  / Nsamp << " us\n";
-  std::cout << "Relative speedup using BVH vs direct DCL        = " << dcelTime.count()/linTime.count() << "\n";
+  std::cout << "Accumulated distance and time using trimesh BVH = " << triSum  << ", which took " << triTime.count()  / Nsamp << " us\n";  
+  std::cout << "Relative speedup using BVH vs direct DCEL        = " << dcelTime.count()/triTime.count() << "\n";
   // clang-format on  
 
   return 0;

Examples/EBGeometry_F18/main.cpp

@@ -21,7 +21,7 @@ using MetaData     = std::map<int, std::vector<unsigned long long>>; // Attach s
 using Vec3         = EBGeometry::Vec3T<T>;
 using BV           = EBGeometry::BoundingVolumes::AABBT<T>;
 using SlowSDF      = EBGeometry::MeshSDF<T, MetaData>;
-using FastSDF      = EBGeometry::FastCompactMeshSDF<T, MetaData, BV, K>;
+using FastSDF      = EBGeometry::FastTriMeshSDF<T, MetaData, BV, K>;
 
 int
 main(int argc, char* argv[])
@@ -137,6 +137,9 @@ main(int argc, char* argv[])
   if (std::abs(sumSlowSlow) - std::abs(sumFastFast) > std::numeric_limits<T>::min()) {
     std::cerr << "Got wrong distance! Diff = " << std::abs(sumSlowFast) - std::abs(sumFastFast) << "\n";
   }
+  if (std::abs(sumSlowSlow) - std::abs(sumFastSlow) > std::numeric_limits<T>::min()) {
+    std::cerr << "Got wrong distance! Diff = " << std::abs(sumSlowFast) - std::abs(sumFastSlow) << "\n";
+  }
 
   const std::chrono::duration<T, std::micro> slowSlowTime = (t1 - t0);
   const std::chrono::duration<T, std::micro> slowFastTime = (t2 - t1);

Source/EBGeometry_CSGImplem.hpp

@@ -330,11 +330,10 @@ FastUnionIF<T, P, BV, K>::value(const Vec3T<T>& a_point) const noexcept
   };
 
   BVH::Sorter<Node, T, K> sorter =
-    [&a_point](std::array<std::pair<std::shared_ptr<const Node>, T>, K>& a_leaves) noexcept -> void {
-    std::sort(
-      a_leaves.begin(),
-      a_leaves.end(),
-      [&a_point](const std::pair<std::shared_ptr<const Node>, T>& n1,
+    [](std::array<std::pair<std::shared_ptr<const Node>, T>, K>& a_leaves) noexcept -> void {
+    std::sort(a_leaves.begin(),
+              a_leaves.end(),
+              [](const std::pair<std::shared_ptr<const Node>, T>& n1,
                  const std::pair<std::shared_ptr<const Node>, T>& n2) -> bool { return n1.second > n2.second; });
   };
 

Source/EBGeometry_DCEL_EdgeImplem.hpp

@@ -234,10 +234,11 @@ namespace DCEL {
   inline T
   EdgeT<T, Meta>::signedDistance(const Vec3& a_x0) const noexcept
   {
+    T retval = std::numeric_limits<T>::max();
+
     // Project point to edge.
     const T t = this->projectPointToEdge(a_x0);
 
-    T retval;
     if (t <= 0.0) {
       // Closest point is the starting vertex
       retval = this->getVertex()->signedDistance(a_x0);

Source/EBGeometry_DCEL_Face.hpp

@@ -259,6 +259,13 @@ namespace DCEL {
     inline std::vector<VertexPtr>
     gatherVertices() const noexcept;
 
+    /*!
+      @brief Return all the half-edges on this polygon
+      @details This builds a list of all the edges and returns it.
+    */
+    inline std::vector<EdgePtr>
+    gatherEdges() const noexcept;
+
     /*!
       @brief Get the lower-left-most coordinate of this polygon face
     */

Source/EBGeometry_DCEL_FaceImplem.hpp

@@ -244,6 +244,21 @@ namespace DCEL {
     return vertices;
   }
 
+  template <class T, class Meta>
+  inline std::vector<std::shared_ptr<EdgeT<T, Meta>>>
+  FaceT<T, Meta>::gatherEdges() const noexcept
+  {
+    std::vector<EdgePtr> edges;
+
+    for (EdgeIterator iter(*this); iter.ok(); ++iter) {
+      EdgePtr& edge = iter();
+
+      edges.emplace_back(edge);
+    }
+
+    return edges;
+  }
+
   template <class T, class Meta>
   inline std::vector<Vec3T<T>>
   FaceT<T, Meta>::getAllVertexCoordinates() const noexcept