From 219030bafbe895ca281276006f15e3468c647d6a Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 12:18:16 +0200 Subject: [PATCH 01/16] Add triangle classes --- EBGeometry.hpp | 2 + Source/EBGeometry_Triangle.hpp | 272 +++++++++++++++++++++++++ Source/EBGeometry_Triangle2D.hpp | 204 +++++++++++++++++++ Source/EBGeometry_Triangle2DImplem.hpp | 215 +++++++++++++++++++ Source/EBGeometry_TriangleImplem.hpp | 235 +++++++++++++++++++++ 5 files changed, 928 insertions(+) create mode 100644 Source/EBGeometry_Triangle.hpp create mode 100644 Source/EBGeometry_Triangle2D.hpp create mode 100644 Source/EBGeometry_Triangle2DImplem.hpp create mode 100644 Source/EBGeometry_TriangleImplem.hpp diff --git a/EBGeometry.hpp b/EBGeometry.hpp index 4adddcae..32e0494d 100644 --- a/EBGeometry.hpp +++ b/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 diff --git a/Source/EBGeometry_Triangle.hpp b/Source/EBGeometry_Triangle.hpp new file mode 100644 index 00000000..2c7e8aa3 --- /dev/null +++ b/Source/EBGeometry_Triangle.hpp @@ -0,0 +1,272 @@ +/* EBGeometry + * Copyright © 2024 Robert Marskar + * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. + */ + +/*! + @file EBGeometry_Triangle.hpp + @brief Declaration of a triangle class with signed distance functionality. + @author Robert Marskar +*/ + +#ifndef EBGeometry_Triangle +#define EBGeometry_Triangle + +// Std includes +#include + +// Our includes +#include "EBGeometry_Vec.hpp" +#include "EBGeometry_Triangle2D.hpp" + +namespace EBGeometry { + + /*! + @brief Triangle class with signed distance functionality. + @details This class represents a planar triangle and has a signed distance functionality. It is + self-contained such that it can be directly copied to GPUs. The class contains a triangle face normal + vector; three vertex positions, and normal vectors for the three vertices and three edges. + + This class assumes that the vertices are organized with the right-hand rule. I.e., edges are enumerated + as follows: + + Edge 1 points from vertex 1 to vertex 2 + Edge 2 points from vertex 2 to vertex 3 + Edge 3 points from vertex 3 to vertex 0 + + This class can compute its own normal vector from the vertex positions, and the triangle orientation + is then implicitly given by the vertex order. + + To compute the distance from a point to the triangle, one must determine if the point projects to the + "inside" or "outside" of the triangle. This class contains a 2D embedding of the triangle that can perform + this project. If the query point projects to the inside of the triangle, the distance is simply the + projected distance onto the triangle plane. If it projects to the outside of the triangle, we check the + distance against the triangle edges and vertices. + the ed + */ + template + class Triangle + { + public: + /*! + @brief Alias for 2D vector type + */ + using Vec2 = Vec2T; + + /*! + @brief Alias for 3D vector type + */ + using Vec3 = Vec3T; + + /*! + @brief Default constructor. Does not put the triangle in a usable state. + */ + Triangle() noexcept = default; + + /*! + @brief Copy constructor. + @param[in] a_otherTriangle Other triangle. + */ + Triangle(const Triangle& a_otherTriangle) noexcept = default; + + /*! + @brief Full constructor. + @param[in] a_vertexPositions Triangle vertex positions + */ + Triangle(const std::array& a_vertexPositions) noexcept; + + /*! + @brief Destructor (does nothing). + */ + ~Triangle() noexcept = default; + + /*! + @brief Copy assignment. + @param[in] a_otherTriangle Other triangle. + */ + Triangle& + operator=(const Triangle& a_otherTriangle) noexcept = default; + + /*! + @brief Move assignment. + @param[in, out] a_otherTriangle Other triangle. + */ + Triangle& + operator=(Triangle&& a_otherTriangle) noexcept = default; + + /*! + @brief Set the triangle normal vector. + @param[in] a_normal Normal vector (should be consistent with the vertex ordering!). + */ + void + setNormal(const Vec3& a_normal) noexcept; + + /*! + @brief Set the triangle vertex positions + @param[in] a_vertexPositions Vertex positions + */ + void + setVertexPositions(const std::array& a_vertexPositions) noexcept; + + /*! + @brief Set the triangle vertex normals + @param[in] a_vertexNormals Vertex normals + */ + void + setVertexNormals(const std::array& a_vertexNormals) noexcept; + + /*! + @brief Set the triangle edge normals + @param[in] a_edgeNormals Edge normals + */ + void + setEdgeNormals(const std::array& a_edgeNormals) noexcept; + + /*! + @brief Set the triangle meta-data + @param[in] a_metaData Triangle metadata. + */ + void + setMetaData(const Meta& a_metaData) noexcept; + + /*! + @brief Compute the triangle normal vector. + @details This computes the normal vector from two of the triangle edges. + */ + void + computeNormal() noexcept; + + /*! + @brief Get the triangle normal vector. + @return m_triangleNormal + */ + Vec3& + getNormal() noexcept; + + /*! + @brief Get the triangle normal vector. + @return m_triangleNormal + */ + const Vec3& + getNormal() const noexcept; + + /*! + @brief Get the vertex positions + @return m_vertexPositions + */ + std::array& + getVertexPositions() noexcept; + + /*! + @brief Get the vertex positions + @return m_vertexPositions + */ + const std::array& + getVertexPositions() const noexcept; + + /*! + @brief Get the vertex normals + @return m_vertexNormals + */ + std::array& + getVertexNormals() noexcept; + + /*! + @brief Get the vertex normals + @return m_vertexNormals + */ + const std::array& + getVertexNormals() const noexcept; + + /*! + @brief Get the edge normals + @return m_edgeNormals + */ + std::array& + getEdgeNormals() noexcept; + + /*! + @brief Get the edge normals + @return m_edgeNormals + */ + const std::array& + getEdgeNormals() const noexcept; + + /*! + @brief Get the triangle meta-data + @return m_metaData + */ + Meta& + getMetaData() noexcept; + + /*! + @brief Get the triangle meta-data + @return m_metaData + */ + const Meta& + getMetaData() const noexcept; + + /*! + @brief Compute the signed distance from the input point x to the triangle + @param[in] a_point Point + */ + T + signedDistance(const Vec3& a_point) const noexcept; + + /*! + @brief Check if a line passes through the triangle. + @details Returns true if the line segment passes through the triangle, edges of the triangle, + or through one of the vertices. + @param[in] a_x0 One endpoint of the line + @param[in] a_x1 Other endpoint of the line + */ + bool + intersects(const Vec3& a_x0, const Vec3& a_x1) const noexcept; + + protected: + /*! + @brief Triangle face normal + */ + Vec3 m_triangleNormal = Vec3::max(); + + /*! + @brief Triangle vertex positions + */ + std::array m_vertexPositions{Vec3::max(), Vec3::max(), Vec3::max()}; + + /*! + @brief Triangle vertex normals + */ + std::array m_vertexNormals{Vec3::max(), Vec3::max(), Vec3::max()}; + /*! + @brief Triangle edge normals + */ + std::array m_edgeNormals{Vec3::max(), Vec3::max(), Vec3::max()}; + + /*! + @brief 2D projection of the triangle to one of the Cartesian coordinate directions + */ + Triangle2D m_triangle2D; + + /*! + @brief Triangle meta-data normals + */ + Meta m_metaData; + + /*! + @brief Returns the "projection" of a point to an edge. + @details This function parametrizes the edge as x(t) = x0 + (x1-x0)*t and + returns where on the this edge the point a_x0 projects. If projects onto the + edge if t = [0,1] and to one of the start/end vertices otherwise. + @param[in] a_point Query point + @param[in] a_x0 Starting edge coordinate + @param[in] a_x1 Final edge coordinate + */ + T + projectPointToEdge(const Vec3& a_point, const Vec3& a_x0, const Vec3& a_x1) const noexcept; + }; +} // namespace EBGeometry + +#include "EBGeometry_TriangleImplem.hpp" // NOLINT + +#endif diff --git a/Source/EBGeometry_Triangle2D.hpp b/Source/EBGeometry_Triangle2D.hpp new file mode 100644 index 00000000..83c75501 --- /dev/null +++ b/Source/EBGeometry_Triangle2D.hpp @@ -0,0 +1,204 @@ +/* EBGeometry + * Copyright © 2024 Robert Marskar + * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. + */ +/*! + @file EBGeometry_Triangle2D.hpp + @brief Declaration of a class that encapsulates the projection of a 3D triangle + into a 2D plane. + @author Robert Marskar +*/ + +#ifndef EBGeometry_Triangle2D +#define EBGeometry_Triangle2D + +// Std includes +#include + +// Our includes +#include "EBGeometry_Vec.hpp" + +namespace EBGeometry { + /*! + @brief Class for embedding a triangle face into 2D. + @details This class is required for determining whether or not a 3D point + projected to the plane of a triangle lies inside or outside the + triangle. To do this we compute the 2D embedding of the triangle face, + reducing the problem to a tractable dimension where we can use well-tested + algorithm. Once the 2D embedding is computed, we can use well-known + algorithms for checking if a point lies inside or outside. The supported + algorithms are 1) The winding number algorithm (computing the winding number), + 2) Computing the subtended angle of the point with the edges of the polygon + (sums to 360 degrees if the point is inside), or computing the crossing number + which checks how many times a ray cast from the point crosses the edges of the + triangle. + */ + template + class Triangle2D + { + public: + /*! + @brief Alias for 2D vector type + */ + using Vec2 = Vec2T; + + /*! + @brief Alias for 3D vector type + */ + using Vec3 = Vec3T; + + /*! + @brief Supported algorithms for performing inside/outside tests when + checking if a point projects to the inside or outside of a polygon face. + */ + enum class InsideOutsideAlgorithm // NOLINT + { + SubtendedAngle, + CrossingNumber, + WindingNumber + }; + + /*! + @brief Default constructor. Must subsequently call the define function. + */ + Triangle2D() noexcept = default; + + /*! + @brief Full constructor + @param[in] a_normal Normal vector of the 3D triangle + @param[in] a_vertices Vertex coordinates of the 3D triangle + */ + Triangle2D(const Vec3& a_normal, const std::array& a_vertices) noexcept; + + /*! + @brief Copy constructor. + @param[in] a_triangle2D Other triangle + */ + Triangle2D(const Triangle2D& a_triangle2D) noexcept = default; + + /*! + @brief Move constructor. + @param[in, out] a_triangle2D Other triangle + */ + Triangle2D(Triangle2D&& a_triangle2D) noexcept = default; + + /*! + @brief Destructor (does nothing) + */ + ~Triangle2D() noexcept = default; + + /*! + @brief Copy assignment. + @param[in] a_triangle2D Other triangle + */ + Triangle2D& + operator=(const Triangle2D& a_triangle2D) noexcept = default; + + /*! + @brief Move assignment. + @param[in, out] a_triangle2D Other triangle + */ + Triangle2D& + operator=(Triangle2D&& a_triangle2D) noexcept = default; + + /*! + @brief Define function. Puts object in usable state. + @param[in] a_normal Normal vector of the 3D triangle + @param[in] a_vertices Vertex coordinates of the 3D triangle + */ + void + define(const Vec3& a_normal, const std::array& a_vertices) noexcept; + + /*! + @brief Check if a point is inside or outside the 2D polygon + @param[in] a_point 3D point coordinates + @param[in] a_algorithm Inside/outside algorithm + @details This will call the function corresponding to a_algorithm. + */ + bool + isPointInside(const Vec3& a_point, InsideOutsideAlgorithm a_algorithm) const noexcept; + + /*! + @brief Check if a point is inside a 2D polygon, using the winding number + algorithm + @param[in] a_point 3D point coordinates + @return Returns true if the 3D point projects to the inside of the 2D + polygon + */ + bool + isPointInsideWindingNumber(const Vec3& a_point) const noexcept; + + /*! + @brief Check if a point is inside a 2D polygon, by computing the number of + times a ray crosses the polygon edges. + @param[in] a_point 3D point coordinates + @return Returns true if the 3D point projects to the inside of the 2D + polygon + */ + bool + isPointInsideCrossingNumber(const Vec3& a_point) const noexcept; + + /*! + @brief Check if a point is inside a 2D polygon, using the subtended angles + @param[in] a_point 3D point coordinates + @return Returns true if the 3D point projects to the inside of the 2D + polygon + */ + bool + isPointInsideSubtend(const Vec3& a_point) const noexcept; + + private: + /*! + @brief The corresponding 2D x-direction (one direction is ignored) + */ + int m_xDir = -1; + + /*! + @brief The corresponding 2D y-direction (one direction is ignored) + */ + int m_yDir = -1; + + /*! + @brief List of points in 2D. + @details This is the position of the vertices, projected into a 2D plane. + */ + std::array m_vertices{Vec2::max(), Vec2::max()}; + + /*! + @brief Project a 3D point onto the 2D polygon plane (this ignores one of the + vector components) + @param[in] a_poitn 3D point + @return 2D point, ignoring one of the coordinate directions. + */ + Vec2 + projectPoint(const Vec3& a_point) const noexcept; + + /*! + @brief Compute the winding number for a point P with the 2D polygon + @param[in] P 2D point + @return Returns winding number. + */ + int + computeWindingNumber(const Vec2& a_point) const noexcept; + + /*! + @brief Compute the crossing number for a point P with the 2D polygon + @param[in] a_point 2D point + @return Returns crossing number. + */ + int + computeCrossingNumber(const Vec2& a_point) const noexcept; + + /*! + @brief Compute the subtended angle for a point a_point with the 2D polygon + @param[in] a_point 2D point + @return Returns subtended angle. + */ + T + computeSubtendedAngle(const Vec2& a_point) const noexcept; + }; +} // namespace EBGeometry + +#include "EBGeometry_Triangle2DImplem.hpp" // NOLINT + +#endif diff --git a/Source/EBGeometry_Triangle2DImplem.hpp b/Source/EBGeometry_Triangle2DImplem.hpp new file mode 100644 index 00000000..590b26f1 --- /dev/null +++ b/Source/EBGeometry_Triangle2DImplem.hpp @@ -0,0 +1,215 @@ +/* EBGeometry + * Copyright © 2024 Robert Marskar + * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. + */ +/*! + @file EBGeometry_Triangle2DImplem.hpp + @brief Implementation of EBGeometry_Triangle2DImplem.hpp + @author Robert Marskar +*/ + +#ifndef EBGeometry_Triangle2DImplem +#define EBGeometry_Triangle2DImplem + +// Std includes +#include + +// Our includes +#include "EBGeometry_Triangle2D.hpp" + +namespace EBGeometry { + + template + Triangle2D::Triangle2D(const Vec3T& a_normal, const std::array, 3>& a_vertices) noexcept + { + this->define(a_normal, a_vertices); + } + + template + void + Triangle2D::define(const Vec3T& a_normal, const std::array, 3>& a_vertices) noexcept + { + int projectDir = 0; + + for (int dir = 1; dir < 3; dir++) { + if (std::abs(a_normal[dir]) > std::abs(a_normal[projectDir])) { + projectDir = dir; + } + } + + m_xDir = 3; + m_yDir = 0; + + for (int dir = 0; dir < 3; dir++) { + if (dir != projectDir) { + m_xDir = std::min(m_xDir, dir); + m_yDir = std::max(m_yDir, dir); + } + } + + for (int i = 0; i < 3; i++) { + m_vertices[i] = this->projectPoint(a_vertices[i]); + } + } + + template + bool + Triangle2D::isPointInside(const Vec3T& a_point, const InsideOutsideAlgorithm a_algorithm) const noexcept + { + bool ret = false; + + switch (a_algorithm) { + case InsideOutsideAlgorithm::SubtendedAngle: { + ret = this->isPointInsideSubtend(a_point); + + break; + } + case InsideOutsideAlgorithm::CrossingNumber: { + ret = this->isPointInsideCrossingNumber(a_point); + + break; + } + case InsideOutsideAlgorithm::WindingNumber: { + ret = this->isPointInsideWindingNumber(a_point); + + break; + } + } + + return ret; + } + + template + bool + Triangle2D::isPointInsideWindingNumber(const Vec3T& a_point) const noexcept + { + const Vec2T projectedPoint = this->projectPoint(a_point); + + const int windingNumber = this->computeWindingNumber(projectedPoint); + + return windingNumber != 0; + } + + template + bool + Triangle2D::isPointInsideCrossingNumber(const Vec3T& a_point) const noexcept + { + const Vec2T projectedPoint = this->projectPoint(a_point); + + const int crossingNumber = this->computeCrossingNumber(projectedPoint); + + return (crossingNumber & 1); + } + + template + bool + Triangle2D::isPointInsideSubtend(const Vec3T& a_point) const noexcept + { + const Vec2T projectedPoint = this->projectPoint(a_point); + + T sumTheta = this->computeSubtendedAngle(projectedPoint); + + sumTheta = std::abs(sumTheta) / (T(2.0) * M_PI); // NOLINT + + return (round(sumTheta) == 1); + } + + template + Vec2T + Triangle2D::projectPoint(const Vec3T& a_point) const noexcept + { + + return Vec2T(a_point[m_xDir], a_point[m_yDir]); + } + + template + int + Triangle2D::computeWindingNumber(const Vec2T& a_point) const noexcept + { + int windingNumber = 0; + + constexpr T zero = T(0.0); + + auto isLeft = [](const Vec2T& p0, const Vec2T& p1, const Vec2T& p2) { + return (p1.x() - p0.x()) * (p2.y() - p0.y()) - (p2.x() - p0.x()) * (p1.y() - p0.y()); + }; + + // Loop through all edges of the polygon + for (int i = 0; i < 3; i++) { + + const Vec2T& P = a_point; + const Vec2T& p1 = m_vertices[i]; + const Vec2T& p2 = m_vertices[(i + 1) % 3]; + + const T res = isLeft(p1, p2, P); + + if (p1.y() <= P.y()) { + if (p2.y() > P.y() && res > zero) { + windingNumber += 1; + } + } + else { + if (p2.y() <= P.y() && res < zero) { + windingNumber -= 1; + } + } + } + + return windingNumber; + } + + template + int + Triangle2D::computeCrossingNumber(const Vec2T& a_point) const noexcept + { + int crossingNumber = 0; + + for (int i = 0; i < 3; i++) { + const Vec2T& p1 = m_vertices[i]; + const Vec2T& p2 = m_vertices[(i + 1) % 3]; + + const bool upwardCrossing = (p1.y() <= a_point.y()) && (p2.y() > a_point.y()); + const bool downwardCrossing = (p1.y() > a_point.y()) && (p2.y() <= a_point.y()); + + if (upwardCrossing || downwardCrossing) { + const T t = (a_point.y() - p1.y()) / (p2.y() - p1.y()); + + if (a_point.x() < p1.x() + t * (p2.x() - p1.x())) { + crossingNumber += 1; + } + } + } + + return crossingNumber; + } + + template + T + Triangle2D::computeSubtendedAngle(const Vec2T& a_point) const noexcept + { + T sumTheta = 0.0; + + for (int i = 0; i < 3; i++) { + const Vec2T p1 = m_vertices[i] - a_point; + const Vec2T p2 = m_vertices[(i + 1) % 3] - a_point; + + const T theta1 = static_cast(atan2(p1.y(), p1.x())); + const T theta2 = static_cast(atan2(p2.y(), p2.x())); + + T dTheta = theta2 - theta1; + + while (dTheta > M_PI) { + dTheta -= 2.0 * M_PI; + } + while (dTheta < -M_PI) { + dTheta += 2.0 * M_PI; + } + + sumTheta += dTheta; + } + + return sumTheta; + } +} // namespace EBGeometry + +#endif diff --git a/Source/EBGeometry_TriangleImplem.hpp b/Source/EBGeometry_TriangleImplem.hpp new file mode 100644 index 00000000..f2203845 --- /dev/null +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -0,0 +1,235 @@ +/* EBGeometry + * Copyright © 2024 Robert Marskar + * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. + */ + +/*! + @file EBGeometry_TriangleImplem.hpp + @brief Implementation of EBGeometry_Triangle.hpp + @author Robert Marskar +*/ + +#ifndef EBGeometry_TriangleImplem +#define EBGeometry_TriangleImplem + +// Our includes +#include "EBGeometry_Triangle.hpp" + +namespace EBGeometry { + + template + Triangle::Triangle(const std::array, 3>& a_vertexPositions) noexcept + { + this->setVertexPositions(a_vertexPositions); + } + + template + void + Triangle::setNormal(const Vec3T& a_normal) noexcept + { + this->m_triangleNormal = a_normal; + } + + template + void + Triangle::setVertexPositions(const std::array, 3>& a_vertexPositions) noexcept + { + for (int i = 0; i < 3; i++) { + this->m_vertexPositions[i] = a_vertexPositions[i]; + } + + this->computeNormal(); + + m_triangle2D.define(m_triangleNormal, m_vertexPositions); + } + + template + void + Triangle::setVertexNormals(const std::array, 3>& a_vertexNormals) noexcept + { + for (int i = 0; i < 3; i++) { + this->m_vertexNormals[i] = a_vertexNormals[i]; + } + } + + template + void + Triangle::setEdgeNormals(const std::array, 3>& a_edgeNormals) noexcept + { + for (int i = 0; i < 3; i++) { + this->m_edgeNormals[i] = a_edgeNormals[i]; + } + } + + template + void + Triangle::setMetaData(const Meta& a_metaData) noexcept + { + this->m_metaData = a_metaData; + } + + template + void + Triangle::computeNormal() noexcept + { + const Vec3T x1x0 = m_vertexPositions[1] - m_vertexPositions[0]; + const Vec3T x2x1 = m_vertexPositions[2] - m_vertexPositions[1]; + + m_triangleNormal = cross(x1x0, x2x1); + m_triangleNormal = m_triangleNormal / m_triangleNormal.length(); + } + + template + Vec3T& + Triangle::getNormal() noexcept + { + return (this->m_triangleNormal); + } + + template + const Vec3T& + Triangle::getNormal() const noexcept + { + return (this->m_triangleNormal); + } + + template + std::array, 3>& + Triangle::getVertexPositions() noexcept + { + return (this->m_vertexPositions); + } + + template + const std::array, 3>& + Triangle::getVertexPositions() const noexcept + { + return (this->m_vertexPositions); + } + + template + std::array, 3>& + Triangle::getVertexNormals() noexcept + { + return (this->m_vertexNormals); + } + + template + const std::array, 3>& + Triangle::getVertexNormals() const noexcept + { + return (this->m_vertexNormals); + } + template + std::array, 3>& + Triangle::getEdgeNormals() noexcept + { + return (this->m_edgeNormals); + } + + template + const std::array, 3>& + Triangle::getEdgeNormals() const noexcept + { + return (this->m_edgeNormals); + } + + template + Meta& + Triangle::getMetaData() noexcept + { + return (this->m_metaData); + } + + template + const Meta& + Triangle::getMetaData() const noexcept + { + return (this->m_metaData); + } + + template + T + Triangle::signedDistance(const Vec3T& a_point) const noexcept + { + // Check if projection falls inside. + const bool isPointInside = m_triangle2D.isPointInsideWindingNumber(a_point); + + T ret = std::numeric_limits::max(); + + if (isPointInside) { + ret = m_triangleNormal.dot(a_point - m_vertexPositions[0]); + } + else { + // If this triggers then + auto sgn = [](const T x) -> int { return (x > 0.0) ? 1 : -1; }; + + // Check distances to vertices. + for (int i = 0; i < 3; i++) { + const Vec3T delta = a_point - m_vertexPositions[i]; + + ret = (delta.length() > std::abs(ret)) ? ret : delta.length() * sgn(m_vertexNormals[i].dot(delta)); + } + + // Check distances to edges + for (int i = 0; i < 3; i++) { + const Vec3T& a = m_vertexPositions[i]; + const Vec3T& b = m_vertexPositions[(i + 1) % 3]; + + const T t = this->projectPointToEdge(a_point, a, b); + + if (t > 0.0 && t < 1.0) { + const Vec3T delta = a_point - (a + t * (b - a)); + + ret = (delta.length() > std::abs(ret)) ? ret : delta.length() * sgn(m_edgeNormals[i].dot(delta)); + } + } + } + + return ret; + } + + template + T + Triangle::projectPointToEdge(const Vec3T& a_point, + const Vec3T& a_x0, + const Vec3T& a_x1) const noexcept + { + const Vec3T a = a_point - a_x0; + const Vec3T b = a_x1 - a_x0; + + EBGEOMETRY_EXPECT(b.length() > std::numeric_limits::min()); + + return a.dot(b) / (b.dot(b)); + } + + template + bool + Triangle::intersects(const Vec3T& a_x0, const Vec3T& a_x1) const noexcept + { + const T epsilon = std::numeric_limits::eps(); + + const Vec3T edge1 = m_vertexPositions[1] - m_vertexPositions[0]; + const Vec3T edge2 = m_vertexPositions[2] - m_vertexPositions[0]; + const Vec3T ray = a_x1 - a_x0; + + const T det = -dot(ray, m_triangleNormal); + const T invDet = T(1.0) / det; + + const Vec3T AO = a_x0 - m_vertexPositions[0]; + const Vec3T DAO = cross(AO, ray); + + const T u = dot(edge2, DAO) * invDet; + const T v = -dot(edge1, DAO) * invDet; + const T t = dot(AO, m_triangleNormal) * invDet; + + const bool a = std::abs(det) > epsilon; + const bool b = (t >= 0.0) && (t <= 1.0); + const bool c = (u >= 0.0) && (std::abs(u - 1.0) >= 0.0); + const bool d = (v >= 0.0) && (std::abs(u + v - 1.0) >= 0.0); + + return (a && b && c && d); + } +} // namespace EBGeometry + +#endif From 8dfb150feb6aec0ac981ef39d360b4b694dd5358 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 13:26:11 +0200 Subject: [PATCH 02/16] Update with trimesh. Seems to run faster --- Examples/EBGeometry_DCEL/main.cpp | 14 +++ Source/EBGeometry_DCEL_Face.hpp | 7 ++ Source/EBGeometry_DCEL_FaceImplem.hpp | 15 +++ Source/EBGeometry_MeshDistanceFunctions.hpp | 75 +++++++++++ ...EBGeometry_MeshDistanceFunctionsImplem.hpp | 116 ++++++++++++++++++ Source/EBGeometry_Parser.hpp | 38 +++++- Source/EBGeometry_ParserImplem.hpp | 66 ++++++++++ Source/EBGeometry_Triangle2DImplem.hpp | 20 +-- Source/EBGeometry_TriangleImplem.hpp | 6 +- 9 files changed, 342 insertions(+), 15 deletions(-) diff --git a/Examples/EBGeometry_DCEL/main.cpp b/Examples/EBGeometry_DCEL/main.cpp index da0f9b85..f6bad6c0 100644 --- a/Examples/EBGeometry_DCEL/main.cpp +++ b/Examples/EBGeometry_DCEL/main.cpp @@ -10,6 +10,8 @@ #include "../../EBGeometry.hpp" +#warning "When the TriMesh class is done, all example files should be updated to use the new functionality." + using namespace EBGeometry; using namespace EBGeometry::DCEL; @@ -21,6 +23,7 @@ using T = float; using Meta = short; using BV = EBGeometry::BoundingVolumes::AABBT; using Vec3 = EBGeometry::Vec3T; +using Tri = EBGeometry::Triangle; int main(int argc, char* argv[]) @@ -48,6 +51,7 @@ main(int argc, char* argv[]) const auto dcelSDF = EBGeometry::Parser::readIntoMesh(file); const auto bvhSDF = EBGeometry::Parser::readIntoFullBVH(file); const auto linSDF = EBGeometry::Parser::readIntoLinearBVH(file); + const auto triSDF = EBGeometry::Parser::readIntoTriangleBVH(file); // Sample some random points around the object. constexpr size_t Nsamp = 100; @@ -74,6 +78,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,10 +93,15 @@ 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 dcelTime = t1 - t0; const std::chrono::duration bvhTime = t2 - t1; const std::chrono::duration linTime = t3 - t2; + const std::chrono::duration triTime = t4 - t3; if (std::abs(bvhSum - dcelSum) > std::numeric_limits::epsilon()) { std::cerr << "Full BVH did not give same distance! Diff = " << bvhSum - dcelSum << "\n"; @@ -99,12 +109,16 @@ main(int argc, char* argv[]) if (std::abs(linSum - dcelSum) > std::numeric_limits::epsilon()) { std::cerr << "Linearized BVH did not give same distance! Diff = " << linSum - dcelSum << "\n"; } + if (std::abs(triSum - dcelSum) > std::numeric_limits::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 << "Accumulated distance and time using trimesh BVH = " << triSum << ", which took " << triTime.count() / Nsamp << " us\n"; std::cout << "Relative speedup using BVH vs direct DCL = " << dcelTime.count()/linTime.count() << "\n"; // clang-format on diff --git a/Source/EBGeometry_DCEL_Face.hpp b/Source/EBGeometry_DCEL_Face.hpp index 50fae6cb..5296493f 100644 --- a/Source/EBGeometry_DCEL_Face.hpp +++ b/Source/EBGeometry_DCEL_Face.hpp @@ -259,6 +259,13 @@ namespace DCEL { inline std::vector 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 + gatherEdges() const noexcept; + /*! @brief Get the lower-left-most coordinate of this polygon face */ diff --git a/Source/EBGeometry_DCEL_FaceImplem.hpp b/Source/EBGeometry_DCEL_FaceImplem.hpp index fe10fcf8..a4a30066 100644 --- a/Source/EBGeometry_DCEL_FaceImplem.hpp +++ b/Source/EBGeometry_DCEL_FaceImplem.hpp @@ -244,6 +244,21 @@ namespace DCEL { return vertices; } + template + inline std::vector>> + FaceT::gatherEdges() const noexcept + { + std::vector edges; + + for (EdgeIterator iter(*this); iter.ok(); ++iter) { + EdgePtr& edge = iter(); + + edges.emplace_back(edge); + } + + return edges; + } + template inline std::vector> FaceT::getAllVertexCoordinates() const noexcept diff --git a/Source/EBGeometry_MeshDistanceFunctions.hpp b/Source/EBGeometry_MeshDistanceFunctions.hpp index 18a3a0dc..5da09927 100644 --- a/Source/EBGeometry_MeshDistanceFunctions.hpp +++ b/Source/EBGeometry_MeshDistanceFunctions.hpp @@ -17,8 +17,11 @@ #include "EBGeometry_SignedDistanceFunction.hpp" #include "EBGeometry_DCEL_Mesh.hpp" #include "EBGeometry_BVH.hpp" +#include "EBGeometry_Triangle.hpp" #include "EBGeometry_NamespaceHeader.hpp" +#warning "Must add the getClosestTriangles functions to TriMesh" + namespace DCEL { /*! @@ -256,6 +259,78 @@ class FastCompactMeshSDF : public SignedDistanceFunction std::shared_ptr m_bvh; }; +/*! + @brief Signed distance function for a triangle mesh. This class uses the full BVH representation. +*/ +template +class FastTriMeshSDF : public SignedDistanceFunction +{ +public: + /*! + @brief Alias for DCEL face type + */ + using Tri = typename EBGeometry::Triangle; + + /*! + @brief Alias for which BVH root node + */ + using Root = typename EBGeometry::BVH::LinearBVH; + + /*! + @brief Alias for linearized BVH + */ + using Node = typename Root::LinearNode; + + /*! + @brief Default disallowed constructor + */ + FastTriMeshSDF() = delete; + + /*! + @brief Full constructor. Takes the input triangles and creates the BVH. + @param[in] a_triangles Input triangle soup + @param[in] a_build Specification of build method. Must be TopDown, Morton, or Nested. + */ + FastTriMeshSDF(const std::vector>& a_triangles, + const BVH::Build a_build = BVH::Build::TopDown) noexcept; + + /*! + @brief Destructor + */ + virtual ~FastTriMeshSDF() = default; + + /*! + @brief Value function + @param[in] a_point Input point. + */ + virtual T + signedDistance(const Vec3T& a_point) const noexcept override; + + /*! + @brief Get the bounding volume hierarchy enclosing the mesh + */ + virtual std::shared_ptr& + getRoot() noexcept; + + /*! + @brief Get the bounding volume hierarchy enclosing the mesh + */ + virtual const std::shared_ptr& + getRoot() const noexcept; + + /*! + @brief Compute bounding volume for this mesh. + */ + BV + computeBoundingVolume() const noexcept; + +protected: + /*! + @brief Bounding volume hierarchy + */ + std::shared_ptr m_bvh; +}; + #include "EBGeometry_NamespaceFooter.hpp" #include "EBGeometry_MeshDistanceFunctionsImplem.hpp" diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index e03a8a28..193bf255 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -60,6 +60,54 @@ DCEL::buildFullBVH(const std::shared_ptr>& a_dc return bvh; } +template +std::shared_ptr, BV, K>> +buildTriMeshFullBVH(const std::vector>>& a_triangles, + const BVH::Build a_build = BVH::Build::TopDown) noexcept +{ + using Prim = EBGeometry::Triangle; + using PrimAndBVList = std::vector, BV>>; + + // Create a pair-wise list of DCEL faces and their bounding volumes. + PrimAndBVList primsAndBVs; + + for (const auto& tri : a_triangles) { + const auto& vertexPositions = tri->getVertexPositions(); + + std::vector> vertices{vertexPositions[0], vertexPositions[1], vertexPositions[2]}; + + primsAndBVs.emplace_back(std::make_pair(tri, BV(vertices))); + } + + // Partition the BVH using the default input arguments. + auto bvh = std::make_shared>(primsAndBVs); + + switch (a_build) { + case BVH::Build::TopDown: { + bvh->topDownSortAndPartition(); + + break; + } + case BVH::Build::Morton: { + bvh->template bottomUpSortAndPartition(); + + break; + } + case BVH::Build::Nested: { + bvh->template bottomUpSortAndPartition(); + + break; + } + default: { + std::cerr << "EBGeometry::DCEL::buildFullBVH - unsupported build method requested" << std::endl; + + break; + } + } + + return bvh; +} + template MeshSDF::MeshSDF(const std::shared_ptr& a_mesh) noexcept { @@ -353,6 +401,74 @@ FastCompactMeshSDF::computeBoundingVolume() const noexcept return m_bvh->getBoundingVolume(); }; +template +FastTriMeshSDF::FastTriMeshSDF(const std::vector>& a_triangles, + const BVH::Build a_build) noexcept +{ + auto bvh = EBGeometry::buildTriMeshFullBVH(a_triangles, a_build); + + m_bvh = bvh->flattenTree(); +} + +template +T +FastTriMeshSDF::signedDistance(const Vec3T& a_point) const noexcept +{ + T minDist = std::numeric_limits::infinity(); + + BVH::Updater updater = [&minDist, + &a_point](const std::vector>& triangles) noexcept -> void { + for (const auto& tri : triangles) { + const T curDist = tri->signedDistance(a_point); + + minDist = std::abs(curDist) < std::abs(minDist) ? curDist : minDist; + } + }; + + BVH::Visiter visiter = [&minDist, &a_point](const Node& a_node, const T& a_bvDist) noexcept -> bool { + return a_bvDist <= std::abs(minDist); + }; + + BVH::Sorter sorter = + [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort( + a_leaves.begin(), + a_leaves.end(), + [&a_point](const std::pair, T>& n1, + const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); + }; + + BVH::MetaUpdater metaUpdater = [&a_point](const Node& a_node) noexcept -> T { + return a_node.getDistanceToBoundingVolume(a_point); + }; + + // Traverse the tree. + m_bvh->traverse(updater, visiter, sorter, metaUpdater); + + return minDist; +} + +template +std::shared_ptr, BV, K>>& +FastTriMeshSDF::getRoot() noexcept +{ + return (m_bvh); +} + +template +const std::shared_ptr, BV, K>>& +FastTriMeshSDF::getRoot() const noexcept +{ + return (m_bvh); +} + +template +BV +FastTriMeshSDF::computeBoundingVolume() const noexcept +{ + return m_bvh->getBoundingVolume(); +}; + #include "EBGeometry_NamespaceFooter.hpp" #endif diff --git a/Source/EBGeometry_Parser.hpp b/Source/EBGeometry_Parser.hpp index 1a9d9854..09b06461 100644 --- a/Source/EBGeometry_Parser.hpp +++ b/Source/EBGeometry_Parser.hpp @@ -19,9 +19,10 @@ #include // Our includes +#include "EBGeometry_BoundingVolumes.hpp" #include "EBGeometry_DCEL_Mesh.hpp" #include "EBGeometry_MeshDistanceFunctions.hpp" -#include "EBGeometry_BoundingVolumes.hpp" +#include "EBGeometry_Triangle.hpp" #include "EBGeometry_NamespaceHeader.hpp" /*! @@ -82,6 +83,22 @@ namespace Parser { inline static std::vector>> readIntoMesh(const std::vector a_files) noexcept; + /*! + @brief Read a file containing a single watertight object and return it as an implicit function. + @param[in] a_filename File name + */ + template + inline static std::vector>> + readIntoTriangles(const std::string a_filename) noexcept; + + /*! + @brief Read multiple files containing single watertight objects and return them as an implicit functions. + @param[in] a_files File names + */ + template + inline static std::vector>>> + readIntoTriangles(const std::vector a_files) noexcept; + /*! @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a full BVH. @param[in] a_filename File name @@ -104,6 +121,25 @@ namespace Parser { inline static std::vector>> readIntoFullBVH(const std::vector a_files) noexcept; + /*! + @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a full BVH. + @param[in] a_filename File name + */ + template , size_t K = 4> + inline static std::shared_ptr> + readIntoTriangleBVH(const std::string a_filename) noexcept; + + /*! + @brief Read multiple files containing single watertight objects and return them as DCEL meshes enclosed in BVHs. + @param[in] a_files File names + */ + template , + size_t K = 4> + inline static std::vector>> + readIntoTriangleBVH(const std::vector a_files) noexcept; + /*! @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a linearized BVH @param[in] a_filename File name diff --git a/Source/EBGeometry_ParserImplem.hpp b/Source/EBGeometry_ParserImplem.hpp index 9c55e6f0..7ebe6cb5 100644 --- a/Source/EBGeometry_ParserImplem.hpp +++ b/Source/EBGeometry_ParserImplem.hpp @@ -93,6 +93,50 @@ Parser::readIntoMesh(const std::vector a_files) noexcept return implicitFunctions; } +template +std::vector>> +Parser::readIntoTriangles(const std::string a_filename) noexcept +{ + const auto mesh = Parser::readIntoDCEL(a_filename); + + std::vector>> triangles; + + for (const auto& f : mesh->getFaces()) { + const auto normal = f->getNormal(); + const auto vertices = f->gatherVertices(); + const auto edges = f->gatherEdges(); + + if (vertices.size() != 3) { + std::cerr << "Parser::readIntoTriangles -- DCEL mesh not composed of only triangles!" << "\n"; + } + + // Create the triangle + auto tri = std::make_shared>(); + + tri->setNormal(normal); + tri->setVertexPositions({vertices[0]->getPosition(), vertices[1]->getPosition(), vertices[2]->getPosition()}); + tri->setVertexNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); + tri->setEdgeNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); + + triangles.emplace_back(tri); + } + + return triangles; +} + +template +std::vector>>> +Parser::readIntoTriangles(const std::vector a_files) noexcept +{ + std::vector>>> triangles; + + for (const auto& file : a_files) { + triangles.emplace_back(Parser::readIntoTriangles(file)); + } + + return triangles; +} + template inline std::shared_ptr> Parser::readIntoFullBVH(const std::string a_filename) noexcept @@ -115,6 +159,28 @@ Parser::readIntoFullBVH(const std::vector a_files) noexcept return implicitFunctions; } +template +inline std::shared_ptr> +Parser::readIntoTriangleBVH(const std::string a_filename) noexcept +{ + auto triangles = EBGeometry::Parser::readIntoTriangles(a_filename); + + return std::make_shared>(triangles); +} + +template +inline std::vector>> +Parser::readIntoTriangleBVH(const std::vector a_files) noexcept +{ + std::vector>> implicitFunctions; + + for (const auto& file : a_files) { + implicitFunctions.emplace_back(Parser::readIntoTriangleBVH(file)); + } + + return implicitFunctions; +} + template inline std::shared_ptr> Parser::readIntoLinearBVH(const std::string a_filename) noexcept diff --git a/Source/EBGeometry_Triangle2DImplem.hpp b/Source/EBGeometry_Triangle2DImplem.hpp index 590b26f1..fea518d4 100644 --- a/Source/EBGeometry_Triangle2DImplem.hpp +++ b/Source/EBGeometry_Triangle2DImplem.hpp @@ -131,7 +131,7 @@ namespace EBGeometry { constexpr T zero = T(0.0); auto isLeft = [](const Vec2T& p0, const Vec2T& p1, const Vec2T& p2) { - return (p1.x() - p0.x()) * (p2.y() - p0.y()) - (p2.x() - p0.x()) * (p1.y() - p0.y()); + return (p1.x - p0.x) * (p2.y - p0.y) - (p2.x - p0.x) * (p1.y - p0.y); }; // Loop through all edges of the polygon @@ -143,13 +143,13 @@ namespace EBGeometry { const T res = isLeft(p1, p2, P); - if (p1.y() <= P.y()) { - if (p2.y() > P.y() && res > zero) { + if (p1.y <= P.y) { + if (p2.y > P.y && res > zero) { windingNumber += 1; } } else { - if (p2.y() <= P.y() && res < zero) { + if (p2.y <= P.y && res < zero) { windingNumber -= 1; } } @@ -168,13 +168,13 @@ namespace EBGeometry { const Vec2T& p1 = m_vertices[i]; const Vec2T& p2 = m_vertices[(i + 1) % 3]; - const bool upwardCrossing = (p1.y() <= a_point.y()) && (p2.y() > a_point.y()); - const bool downwardCrossing = (p1.y() > a_point.y()) && (p2.y() <= a_point.y()); + const bool upwardCrossing = (p1.y <= a_point.y) && (p2.y > a_point.y); + const bool downwardCrossing = (p1.y > a_point.y) && (p2.y <= a_point.y); if (upwardCrossing || downwardCrossing) { - const T t = (a_point.y() - p1.y()) / (p2.y() - p1.y()); + const T t = (a_point.y - p1.y) / (p2.y - p1.y); - if (a_point.x() < p1.x() + t * (p2.x() - p1.x())) { + if (a_point.x < p1.x + t * (p2.x - p1.x)) { crossingNumber += 1; } } @@ -193,8 +193,8 @@ namespace EBGeometry { const Vec2T p1 = m_vertices[i] - a_point; const Vec2T p2 = m_vertices[(i + 1) % 3] - a_point; - const T theta1 = static_cast(atan2(p1.y(), p1.x())); - const T theta2 = static_cast(atan2(p2.y(), p2.x())); + const T theta1 = static_cast(atan2(p1.y, p1.x)); + const T theta2 = static_cast(atan2(p2.y, p2.x)); T dTheta = theta2 - theta1; diff --git a/Source/EBGeometry_TriangleImplem.hpp b/Source/EBGeometry_TriangleImplem.hpp index f2203845..5f934462 100644 --- a/Source/EBGeometry_TriangleImplem.hpp +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -75,7 +75,7 @@ namespace EBGeometry { const Vec3T x1x0 = m_vertexPositions[1] - m_vertexPositions[0]; const Vec3T x2x1 = m_vertexPositions[2] - m_vertexPositions[1]; - m_triangleNormal = cross(x1x0, x2x1); + m_triangleNormal = x1x0.cross(x2x1); m_triangleNormal = m_triangleNormal / m_triangleNormal.length(); } @@ -198,8 +198,6 @@ namespace EBGeometry { const Vec3T a = a_point - a_x0; const Vec3T b = a_x1 - a_x0; - EBGEOMETRY_EXPECT(b.length() > std::numeric_limits::min()); - return a.dot(b) / (b.dot(b)); } @@ -217,7 +215,7 @@ namespace EBGeometry { const T invDet = T(1.0) / det; const Vec3T AO = a_x0 - m_vertexPositions[0]; - const Vec3T DAO = cross(AO, ray); + const Vec3T DAO = AO.cross(ray); const T u = dot(edge2, DAO) * invDet; const T v = -dot(edge1, DAO) * invDet; From 13266b87acd6085461159e7613e37760e07ceabd Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 13:36:22 +0200 Subject: [PATCH 03/16] WIP --- Examples/AMReX_DCEL/main.cpp | 2 + Examples/AMReX_F18/main.cpp | 2 + Examples/AMReX_PaintEB/main.cpp | 2 + Examples/Chombo3_DCEL/main.cpp | 2 + Examples/Chombo3_F18/main.cpp | 2 + Examples/EBGeometry_DCEL/main.cpp | 2 +- Examples/EBGeometry_F18/main.cpp | 2 + Source/EBGeometry_MeshDistanceFunctions.hpp | 13 +++- ...EBGeometry_MeshDistanceFunctionsImplem.hpp | 67 ++++++++++++++++++- 9 files changed, 90 insertions(+), 4 deletions(-) diff --git a/Examples/AMReX_DCEL/main.cpp b/Examples/AMReX_DCEL/main.cpp index d5b25029..18c76ef8 100644 --- a/Examples/AMReX_DCEL/main.cpp +++ b/Examples/AMReX_DCEL/main.cpp @@ -13,6 +13,8 @@ // Our include #include "../../EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + using namespace amrex; /*! diff --git a/Examples/AMReX_F18/main.cpp b/Examples/AMReX_F18/main.cpp index b0c14dc2..082eb2e1 100644 --- a/Examples/AMReX_F18/main.cpp +++ b/Examples/AMReX_F18/main.cpp @@ -16,6 +16,8 @@ // Our include #include "../../EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + using namespace amrex; constexpr size_t K = 4; diff --git a/Examples/AMReX_PaintEB/main.cpp b/Examples/AMReX_PaintEB/main.cpp index 3eb56961..18dd09a0 100644 --- a/Examples/AMReX_PaintEB/main.cpp +++ b/Examples/AMReX_PaintEB/main.cpp @@ -13,6 +13,8 @@ // Our include #include "../../EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + using namespace amrex; /*! diff --git a/Examples/Chombo3_DCEL/main.cpp b/Examples/Chombo3_DCEL/main.cpp index 5a842ec8..7ef1c0f9 100644 --- a/Examples/Chombo3_DCEL/main.cpp +++ b/Examples/Chombo3_DCEL/main.cpp @@ -18,6 +18,8 @@ // Our includes #include "EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + // Binding for exposing EBGeometry's signed distance functions to Chombo template class ChomboSDF : public BaseIF diff --git a/Examples/Chombo3_F18/main.cpp b/Examples/Chombo3_F18/main.cpp index e87922a8..18264906 100644 --- a/Examples/Chombo3_F18/main.cpp +++ b/Examples/Chombo3_F18/main.cpp @@ -16,6 +16,8 @@ // Our includes #include "EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + constexpr size_t K = 4; using T = Real; using MetaData = EBGeometry::DCEL::DefaultMetaData; diff --git a/Examples/EBGeometry_DCEL/main.cpp b/Examples/EBGeometry_DCEL/main.cpp index f6bad6c0..240fcb87 100644 --- a/Examples/EBGeometry_DCEL/main.cpp +++ b/Examples/EBGeometry_DCEL/main.cpp @@ -119,7 +119,7 @@ main(int argc, char* argv[]) 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 << "Accumulated distance and time using trimesh BVH = " << triSum << ", which took " << triTime.count() / Nsamp << " us\n"; - std::cout << "Relative speedup using BVH vs direct DCL = " << dcelTime.count()/linTime.count() << "\n"; + std::cout << "Relative speedup using BVH vs direct DCEL = " << dcelTime.count()/triTime.count() << "\n"; // clang-format on return 0; diff --git a/Examples/EBGeometry_F18/main.cpp b/Examples/EBGeometry_F18/main.cpp index 6b7ea1d0..cf8d3d57 100644 --- a/Examples/EBGeometry_F18/main.cpp +++ b/Examples/EBGeometry_F18/main.cpp @@ -12,6 +12,8 @@ #include "../../EBGeometry.hpp" +#warning "Example must be updated to use the trimesh class" + using namespace EBGeometry; using namespace EBGeometry::DCEL; diff --git a/Source/EBGeometry_MeshDistanceFunctions.hpp b/Source/EBGeometry_MeshDistanceFunctions.hpp index 5da09927..30b4d38c 100644 --- a/Source/EBGeometry_MeshDistanceFunctions.hpp +++ b/Source/EBGeometry_MeshDistanceFunctions.hpp @@ -20,8 +20,6 @@ #include "EBGeometry_Triangle.hpp" #include "EBGeometry_NamespaceHeader.hpp" -#warning "Must add the getClosestTriangles functions to TriMesh" - namespace DCEL { /*! @@ -306,6 +304,17 @@ class FastTriMeshSDF : public SignedDistanceFunction virtual T signedDistance(const Vec3T& a_point) const noexcept override; + /*! + @brief Get the closest triangles to the input point + @details This returns a list of candidate triangles that are close to the input point. The returned + argment consists of the triangles and the unsigned distance to the triangles. + @param[in] a_point Input point + @param[in] a_sorted Sort the output vector by distance or not. Closest go first. + @return List of candidate triangles (potentially sorted) + */ + virtual std::vector, T>> + getClosestTriangles(const Vec3T& a_point, const bool a_sorted) const noexcept; + /*! @brief Get the bounding volume hierarchy enclosing the mesh */ diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index 193bf255..7405ec40 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -184,7 +184,6 @@ template std::vector>, T>> FastMeshSDF::getClosestFaces(const Vec3T& a_point, const bool a_sorted) const noexcept { - using FaceAndDist = std::pair, T>; // List of candidate faces. @@ -448,6 +447,72 @@ FastTriMeshSDF::signedDistance(const Vec3T& a_point) const no return minDist; } +template +std::vector>, T>> +FastTriMeshSDF::getClosestTriangles(const Vec3T& a_point, const bool a_sorted) const noexcept +{ + using TriAndDist = std::pair, T>; + + // List of candidate faces. + std::vector candidateTriangles; + + // Declaration of the BVH metadata attached to each node - this will be the distance to the node itself. + using BVHMeta = T; + + // Shortest distance so far. + BVHMeta shortestDistanceSoFar = std::numeric_limits::infinity(); + + // Visitation pattern - go into the node if the point is inside or the distance to the BV is shorter than + // the shortest distance so far. + EBGeometry::BVH::Visiter visiter = [&shortestDistanceSoFar](const Node& a_node, + const BVHMeta& a_bvDist) noexcept -> bool { + return a_bvDist <= 0.0 || a_bvDist <= shortestDistanceSoFar; + }; + + // Sorter for BVH nodes, visit closest nodes first + EBGeometry::BVH::Sorter sorter = + [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort( + a_leaves.begin(), + a_leaves.end(), + [&a_point](const std::pair, T>& n1, + const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); + }; + + // Meta-data updater - this meta-data enters into the visitor pattern. + EBGeometry::BVH::MetaUpdater metaUpdater = [&a_point](const Node& a_node) noexcept -> BVHMeta { + return a_node.getDistanceToBoundingVolume(a_point); + }; + + // Update rule for the BVH. Go through the faces and check + EBGeometry::BVH::Updater updater = [&shortestDistanceSoFar, &a_point, &candidateTriangles]( + const std::vector>& a_faces) noexcept -> void { + // Calculate the distance to each face in the leaf node. If it is shorter than the shortest distance so far, add this face + // to the list of faces and update the shortest distance. + for (const auto& f : a_faces) { + const T distToTri = std::abs(f->signedDistance(a_point)); + + if (distToTri <= shortestDistanceSoFar) { + candidateTriangles.emplace_back(f, distToTri); + + shortestDistanceSoFar = distToTri; + } + } + }; + + // Traverse the tree + m_bvh->traverse(updater, visiter, sorter, metaUpdater); + + // Sort if the user asks for it + if (a_sorted) { + std::sort(candidateTriangles.begin(), candidateTriangles.end(), [](const TriAndDist& a, const TriAndDist& b) { + return a.second < b.second; + }); + } + + return candidateTriangles; +} + template std::shared_ptr, BV, K>>& FastTriMeshSDF::getRoot() noexcept From 7658136f7521a307de39a0850db34b092759fd73 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 13:54:05 +0200 Subject: [PATCH 04/16] update --- Docs/Sphinx/source/Concepts.rst | 7 ++++ Docs/Sphinx/source/Parsers.rst | 40 ++++++++++++++++--- Source/EBGeometry_Parser.hpp | 58 ++++++++++++++-------------- Source/EBGeometry_Triangle.hpp | 10 ----- Source/EBGeometry_TriangleImplem.hpp | 28 -------------- 5 files changed, 71 insertions(+), 72 deletions(-) diff --git a/Docs/Sphinx/source/Concepts.rst b/Docs/Sphinx/source/Concepts.rst index 092dd4b4..cf9a86ed 100644 --- a/Docs/Sphinx/source/Concepts.rst +++ b/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 diff --git a/Docs/Sphinx/source/Parsers.rst b/Docs/Sphinx/source/Parsers.rst index 154fd0f6..0badd6c1 100644 --- a/Docs/Sphinx/source/Parsers.rst +++ b/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 files; // <---- List of file names. + + const auto distanceFields = EBGeometry::Parser::readIntoTriangleBVH(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,7 @@ 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 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. diff --git a/Source/EBGeometry_Parser.hpp b/Source/EBGeometry_Parser.hpp index 09b06461..bed4ac9e 100644 --- a/Source/EBGeometry_Parser.hpp +++ b/Source/EBGeometry_Parser.hpp @@ -84,42 +84,48 @@ namespace Parser { readIntoMesh(const std::vector a_files) noexcept; /*! - @brief Read a file containing a single watertight object and return it as an implicit function. + @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a full BVH. @param[in] a_filename File name */ - template - inline static std::vector>> - readIntoTriangles(const std::string a_filename) noexcept; + template , + size_t K = 4> + inline static std::shared_ptr> + readIntoFullBVH(const std::string a_filename) noexcept; /*! - @brief Read multiple files containing single watertight objects and return them as an implicit functions. + @brief Read multiple files containing single watertight objects and return them as DCEL meshes enclosed in BVHs. @param[in] a_files File names */ - template - inline static std::vector>>> - readIntoTriangles(const std::vector a_files) noexcept; + template , + size_t K = 4> + inline static std::vector>> + readIntoFullBVH(const std::vector a_files) noexcept; /*! - @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a full BVH. + @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a linearized BVH @param[in] a_filename File name */ template , size_t K = 4> - inline static std::shared_ptr> - readIntoFullBVH(const std::string a_filename) noexcept; + inline static std::shared_ptr> + readIntoLinearBVH(const std::string a_filename) noexcept; /*! - @brief Read multiple files containing single watertight objects and return them as DCEL meshes enclosed in BVHs. + @brief Read multiple files containing single watertight objects and return them as DCEL meshes enclosed in linearized BVHs. @param[in] a_files File names */ template , size_t K = 4> - inline static std::vector>> - readIntoFullBVH(const std::vector a_files) noexcept; + inline static std::vector>> + readIntoLinearBVH(const std::vector a_files) noexcept; /*! @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a full BVH. @@ -138,29 +144,23 @@ namespace Parser { typename BV = EBGeometry::BoundingVolumes::AABBT, size_t K = 4> inline static std::vector>> - readIntoTriangleBVH(const std::vector a_files) noexcept; + readIntoTriangleBVH(const std::vector a_files) noexcept; /*! - @brief Read a file containing a single watertight object and return it as a DCEL mesh enclosed in a linearized BVH + @brief Read a file containing a single watertight object and return it as an implicit function. @param[in] a_filename File name */ - template , - size_t K = 4> - inline static std::shared_ptr> - readIntoLinearBVH(const std::string a_filename) noexcept; + template + inline static std::vector>> + readIntoTriangles(const std::string a_filename) noexcept; /*! - @brief Read multiple files containing single watertight objects and return them as DCEL meshes enclosed in linearized BVHs. + @brief Read multiple files containing single watertight objects and return them as an implicit functions. @param[in] a_files File names */ - template , - size_t K = 4> - inline static std::vector>> - readIntoLinearBVH(const std::vector a_files) noexcept; + template + inline static std::vector>>> + readIntoTriangles(const std::vector a_files) noexcept; /*! @brief Get file type diff --git a/Source/EBGeometry_Triangle.hpp b/Source/EBGeometry_Triangle.hpp index 2c7e8aa3..efb71f0b 100644 --- a/Source/EBGeometry_Triangle.hpp +++ b/Source/EBGeometry_Triangle.hpp @@ -213,16 +213,6 @@ namespace EBGeometry { T signedDistance(const Vec3& a_point) const noexcept; - /*! - @brief Check if a line passes through the triangle. - @details Returns true if the line segment passes through the triangle, edges of the triangle, - or through one of the vertices. - @param[in] a_x0 One endpoint of the line - @param[in] a_x1 Other endpoint of the line - */ - bool - intersects(const Vec3& a_x0, const Vec3& a_x1) const noexcept; - protected: /*! @brief Triangle face normal diff --git a/Source/EBGeometry_TriangleImplem.hpp b/Source/EBGeometry_TriangleImplem.hpp index 5f934462..c5f5e875 100644 --- a/Source/EBGeometry_TriangleImplem.hpp +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -200,34 +200,6 @@ namespace EBGeometry { return a.dot(b) / (b.dot(b)); } - - template - bool - Triangle::intersects(const Vec3T& a_x0, const Vec3T& a_x1) const noexcept - { - const T epsilon = std::numeric_limits::eps(); - - const Vec3T edge1 = m_vertexPositions[1] - m_vertexPositions[0]; - const Vec3T edge2 = m_vertexPositions[2] - m_vertexPositions[0]; - const Vec3T ray = a_x1 - a_x0; - - const T det = -dot(ray, m_triangleNormal); - const T invDet = T(1.0) / det; - - const Vec3T AO = a_x0 - m_vertexPositions[0]; - const Vec3T DAO = AO.cross(ray); - - const T u = dot(edge2, DAO) * invDet; - const T v = -dot(edge1, DAO) * invDet; - const T t = dot(AO, m_triangleNormal) * invDet; - - const bool a = std::abs(det) > epsilon; - const bool b = (t >= 0.0) && (t <= 1.0); - const bool c = (u >= 0.0) && (std::abs(u - 1.0) >= 0.0); - const bool d = (v >= 0.0) && (std::abs(u + v - 1.0) >= 0.0); - - return (a && b && c && d); - } } // namespace EBGeometry #endif From d5223e3393e76a956619faf75dc5d5229cc98fb0 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 13:58:30 +0200 Subject: [PATCH 05/16] Update docs --- Docs/Sphinx/source/ImplemBVH.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Docs/Sphinx/source/ImplemBVH.rst b/Docs/Sphinx/source/ImplemBVH.rst index 836fff33..e5280961 100644 --- a/Docs/Sphinx/source/ImplemBVH.rst +++ b/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. From 5801d7100cf636fd430a61773bf63dae40c5d01d Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 14:14:34 +0200 Subject: [PATCH 06/16] Update --- Docs/Sphinx/source/Parsers.rst | 1 + Source/EBGeometry_MeshDistanceFunctions.hpp | 14 ++++++++- ...EBGeometry_MeshDistanceFunctionsImplem.hpp | 31 +++++++++++++++++++ Source/EBGeometry_Parser.hpp | 4 +-- Source/EBGeometry_ParserImplem.hpp | 4 +-- 5 files changed, 49 insertions(+), 5 deletions(-) diff --git a/Docs/Sphinx/source/Parsers.rst b/Docs/Sphinx/source/Parsers.rst index 0badd6c1..a77812e7 100644 --- a/Docs/Sphinx/source/Parsers.rst +++ b/Docs/Sphinx/source/Parsers.rst @@ -125,6 +125,7 @@ To turn this into a DCEL mesh, one should compress the triangle soup (get rid of .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp :language: c++ :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. diff --git a/Source/EBGeometry_MeshDistanceFunctions.hpp b/Source/EBGeometry_MeshDistanceFunctions.hpp index 30b4d38c..a6714bcc 100644 --- a/Source/EBGeometry_MeshDistanceFunctions.hpp +++ b/Source/EBGeometry_MeshDistanceFunctions.hpp @@ -264,6 +264,11 @@ template class FastTriMeshSDF : public SignedDistanceFunction { public: + /*! + @brief Alias for DCEL mesh type + */ + using Mesh = EBGeometry::DCEL::MeshT; + /*! @brief Alias for DCEL face type */ @@ -284,6 +289,13 @@ class FastTriMeshSDF : public SignedDistanceFunction */ FastTriMeshSDF() = delete; + /*! + @brief Full constructor. Takes a DCEL mesh and creates the input triangles. Then creates the BVH. + @param[in] a_mesh DCEL mesh + @param[in] a_build Specification of build method. Must be TopDown, Morton, or Nested. + */ + FastTriMeshSDF(const std::shared_ptr& a_mesh, const BVH::Build a_build = BVH::Build::TopDown) noexcept; + /*! @brief Full constructor. Takes the input triangles and creates the BVH. @param[in] a_triangles Input triangle soup @@ -313,7 +325,7 @@ class FastTriMeshSDF : public SignedDistanceFunction @return List of candidate triangles (potentially sorted) */ virtual std::vector, T>> - getClosestTriangles(const Vec3T& a_point, const bool a_sorted) const noexcept; + getClosestTriangles(const Vec3T& a_point, const bool a_sorted) const noexcept; /*! @brief Get the bounding volume hierarchy enclosing the mesh diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index 7405ec40..0cc1a22d 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -400,6 +400,37 @@ FastCompactMeshSDF::computeBoundingVolume() const noexcept return m_bvh->getBoundingVolume(); }; +template +FastTriMeshSDF::FastTriMeshSDF(const std::shared_ptr& a_mesh, const BVH::Build a_build) noexcept +{ + // Turn the input mesh into triangles + std::vector>> triangles; + + for (const auto& f : a_mesh->getFaces()) { + const auto normal = f->getNormal(); + const auto vertices = f->gatherVertices(); + const auto edges = f->gatherEdges(); + + if (vertices.size() != 3) { + std::cerr << "Parser::readIntoTriangles -- DCEL mesh not composed of only triangles!" << "\n"; + } + + // Create the triangle + auto tri = std::make_shared>(); + + tri->setNormal(normal); + tri->setVertexPositions({vertices[0]->getPosition(), vertices[1]->getPosition(), vertices[2]->getPosition()}); + tri->setVertexNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); + tri->setEdgeNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); + + triangles.emplace_back(tri); + } + + auto bvh = EBGeometry::buildTriMeshFullBVH(triangles, a_build); + + m_bvh = bvh->flattenTree(); +} + template FastTriMeshSDF::FastTriMeshSDF(const std::vector>& a_triangles, const BVH::Build a_build) noexcept diff --git a/Source/EBGeometry_Parser.hpp b/Source/EBGeometry_Parser.hpp index bed4ac9e..5ae8f109 100644 --- a/Source/EBGeometry_Parser.hpp +++ b/Source/EBGeometry_Parser.hpp @@ -144,7 +144,7 @@ namespace Parser { typename BV = EBGeometry::BoundingVolumes::AABBT, size_t K = 4> inline static std::vector>> - readIntoTriangleBVH(const std::vector a_files) noexcept; + readIntoTriangleBVH(const std::vector a_files) noexcept; /*! @brief Read a file containing a single watertight object and return it as an implicit function. @@ -160,7 +160,7 @@ namespace Parser { */ template inline static std::vector>>> - readIntoTriangles(const std::vector a_files) noexcept; + readIntoTriangles(const std::vector a_files) noexcept; /*! @brief Get file type diff --git a/Source/EBGeometry_ParserImplem.hpp b/Source/EBGeometry_ParserImplem.hpp index 7ebe6cb5..a4bd77bd 100644 --- a/Source/EBGeometry_ParserImplem.hpp +++ b/Source/EBGeometry_ParserImplem.hpp @@ -163,9 +163,9 @@ template inline std::shared_ptr> Parser::readIntoTriangleBVH(const std::string a_filename) noexcept { - auto triangles = EBGeometry::Parser::readIntoTriangles(a_filename); + const auto mesh = EBGeometry::Parser::readIntoDCEL(a_filename); - return std::make_shared>(triangles); + return std::make_shared>(mesh); } template From acdd08af0def07da0fab477784a65c0a5e9809a8 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 14:15:49 +0200 Subject: [PATCH 07/16] Update examples --- Examples/EBGeometry_DCEL/main.cpp | 2 -- Examples/EBGeometry_F18/main.cpp | 4 +--- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/Examples/EBGeometry_DCEL/main.cpp b/Examples/EBGeometry_DCEL/main.cpp index 240fcb87..0315f72d 100644 --- a/Examples/EBGeometry_DCEL/main.cpp +++ b/Examples/EBGeometry_DCEL/main.cpp @@ -10,8 +10,6 @@ #include "../../EBGeometry.hpp" -#warning "When the TriMesh class is done, all example files should be updated to use the new functionality." - using namespace EBGeometry; using namespace EBGeometry::DCEL; diff --git a/Examples/EBGeometry_F18/main.cpp b/Examples/EBGeometry_F18/main.cpp index cf8d3d57..f2884e0c 100644 --- a/Examples/EBGeometry_F18/main.cpp +++ b/Examples/EBGeometry_F18/main.cpp @@ -12,8 +12,6 @@ #include "../../EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - using namespace EBGeometry; using namespace EBGeometry::DCEL; @@ -23,7 +21,7 @@ using MetaData = std::map>; // Attach s using Vec3 = EBGeometry::Vec3T; using BV = EBGeometry::BoundingVolumes::AABBT; using SlowSDF = EBGeometry::MeshSDF; -using FastSDF = EBGeometry::FastCompactMeshSDF; +using FastSDF = EBGeometry::FastTriMeshSDF; int main(int argc, char* argv[]) From 568eff063916f26c226c5fe25b39f8b7e311b50a Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 14:24:31 +0200 Subject: [PATCH 08/16] Update -- lots of stuff to fix yet. But we're getting there --- Examples/AMReX_DCEL/main.cpp | 4 +--- Examples/AMReX_PaintEB/main.cpp | 3 ++- Source/EBGeometry_MeshDistanceFunctionsImplem.hpp | 2 ++ 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/Examples/AMReX_DCEL/main.cpp b/Examples/AMReX_DCEL/main.cpp index 18c76ef8..4bec038c 100644 --- a/Examples/AMReX_DCEL/main.cpp +++ b/Examples/AMReX_DCEL/main.cpp @@ -13,8 +13,6 @@ // Our include #include "../../EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - using namespace amrex; /*! @@ -33,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(a_filename); + m_sdf = EBGeometry::Parser::readIntoTriangleBVH(a_filename); } else { m_sdf = EBGeometry::Parser::readIntoMesh(a_filename); diff --git a/Examples/AMReX_PaintEB/main.cpp b/Examples/AMReX_PaintEB/main.cpp index 18dd09a0..d0a36eac 100644 --- a/Examples/AMReX_PaintEB/main.cpp +++ b/Examples/AMReX_PaintEB/main.cpp @@ -39,13 +39,14 @@ class AMReXSDF { // Read in the mesh into a DCEL mesh and partition it into a bounding volume hierarchy auto mesh = EBGeometry::Parser::readIntoDCEL(a_filename); - m_sdf = std::make_shared>(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>(mesh); } /*! diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index 0cc1a22d..dac21cab 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -423,6 +423,8 @@ FastTriMeshSDF::FastTriMeshSDF(const std::shared_ptr& a_me tri->setVertexNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); tri->setEdgeNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); +#warning "Must set metadata here" + triangles.emplace_back(tri); } From 32e3cdad0f6a8b744e4c9dfa4a4c30df56ee1144 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 14:52:41 +0200 Subject: [PATCH 09/16] Update --- Examples/AMReX_F18/main.cpp | 4 +--- Examples/AMReX_PaintEB/main.cpp | 4 +--- Examples/Chombo3_DCEL/main.cpp | 4 +--- Examples/Chombo3_F18/main.cpp | 4 +--- Source/EBGeometry_MeshDistanceFunctionsImplem.hpp | 7 +++---- 5 files changed, 7 insertions(+), 16 deletions(-) diff --git a/Examples/AMReX_F18/main.cpp b/Examples/AMReX_F18/main.cpp index 082eb2e1..972c9a56 100644 --- a/Examples/AMReX_F18/main.cpp +++ b/Examples/AMReX_F18/main.cpp @@ -16,8 +16,6 @@ // Our include #include "../../EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - using namespace amrex; constexpr size_t K = 4; @@ -27,7 +25,7 @@ using Vec3 = EBGeometry::Vec3T; using BV = EBGeometry::BoundingVolumes::AABBT; using Face = EBGeometry::DCEL::FaceT; using Mesh = EBGeometry::DCEL::MeshT; -using FastSDF = EBGeometry::FastCompactMeshSDF; +using FastSDF = EBGeometry::FastTriMeshSDF; // F18 geometry, using nifty EBGeometry bindings and accelerators. class F18 diff --git a/Examples/AMReX_PaintEB/main.cpp b/Examples/AMReX_PaintEB/main.cpp index d0a36eac..921e7231 100644 --- a/Examples/AMReX_PaintEB/main.cpp +++ b/Examples/AMReX_PaintEB/main.cpp @@ -13,8 +13,6 @@ // Our include #include "../../EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - using namespace amrex; /*! @@ -39,7 +37,7 @@ class AMReXSDF { // Read in the mesh into a DCEL mesh and partition it into a bounding volume hierarchy auto mesh = EBGeometry::Parser::readIntoDCEL(a_filename); - + // 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++) { diff --git a/Examples/Chombo3_DCEL/main.cpp b/Examples/Chombo3_DCEL/main.cpp index 7ef1c0f9..30b1c8d3 100644 --- a/Examples/Chombo3_DCEL/main.cpp +++ b/Examples/Chombo3_DCEL/main.cpp @@ -18,8 +18,6 @@ // Our includes #include "EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - // Binding for exposing EBGeometry's signed distance functions to Chombo template class ChomboSDF : public BaseIF @@ -29,7 +27,7 @@ class ChomboSDF : public BaseIF ChomboSDF(const std::string a_filename) { - m_implicitFunction = EBGeometry::Parser::readIntoLinearBVH(a_filename); + m_implicitFunction = EBGeometry::Parser::readIntoTriangleBVH(a_filename); m_implicitFunction = EBGeometry::Complement(m_implicitFunction); } diff --git a/Examples/Chombo3_F18/main.cpp b/Examples/Chombo3_F18/main.cpp index 18264906..99edc583 100644 --- a/Examples/Chombo3_F18/main.cpp +++ b/Examples/Chombo3_F18/main.cpp @@ -16,15 +16,13 @@ // Our includes #include "EBGeometry.hpp" -#warning "Example must be updated to use the trimesh class" - constexpr size_t K = 4; using T = Real; using MetaData = EBGeometry::DCEL::DefaultMetaData; using Face = EBGeometry::DCEL::FaceT; using Mesh = EBGeometry::DCEL::MeshT; using BV = EBGeometry::BoundingVolumes::AABBT; -using FastSDF = EBGeometry::FastCompactMeshSDF; +using FastSDF = EBGeometry::FastTriMeshBVH; using Vec3 = EBGeometry::Vec3T; class F18 : public BaseIF diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index dac21cab..0c72a10a 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -411,7 +411,7 @@ FastTriMeshSDF::FastTriMeshSDF(const std::shared_ptr& a_me const auto vertices = f->gatherVertices(); const auto edges = f->gatherEdges(); - if (vertices.size() != 3) { + if ((vertices.size() != 3) || (edges.size() != 3)) { std::cerr << "Parser::readIntoTriangles -- DCEL mesh not composed of only triangles!" << "\n"; } @@ -421,9 +421,8 @@ FastTriMeshSDF::FastTriMeshSDF(const std::shared_ptr& a_me tri->setNormal(normal); tri->setVertexPositions({vertices[0]->getPosition(), vertices[1]->getPosition(), vertices[2]->getPosition()}); tri->setVertexNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); - tri->setEdgeNormals({vertices[0]->getNormal(), vertices[1]->getNormal(), vertices[2]->getNormal()}); - -#warning "Must set metadata here" + tri->setEdgeNormals({edges[0]->getNormal(), edges[1]->getNormal(), edges[2]->getNormal()}); + tri->setMetaData(f->getMetaData()); triangles.emplace_back(tri); } From 323ab3d4e05b5541eaa5cfc356e96ada236947ba Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:02:19 +0200 Subject: [PATCH 10/16] Initializer for DCEL edge --- Source/EBGeometry_DCEL_EdgeImplem.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Source/EBGeometry_DCEL_EdgeImplem.hpp b/Source/EBGeometry_DCEL_EdgeImplem.hpp index 44c45aed..085f97e5 100644 --- a/Source/EBGeometry_DCEL_EdgeImplem.hpp +++ b/Source/EBGeometry_DCEL_EdgeImplem.hpp @@ -234,10 +234,11 @@ namespace DCEL { inline T EdgeT::signedDistance(const Vec3& a_x0) const noexcept { + T retval = std::numeric_limits::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); From 53b4d6069eaf9eb297b26c284bde8c41dfa08b43 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:05:54 +0200 Subject: [PATCH 11/16] commit before removing debug code --- Source/EBGeometry_Triangle2DImplem.hpp | 9 +++--- Source/EBGeometry_TriangleImplem.hpp | 44 ++++++++++++++++++-------- 2 files changed, 34 insertions(+), 19 deletions(-) diff --git a/Source/EBGeometry_Triangle2DImplem.hpp b/Source/EBGeometry_Triangle2DImplem.hpp index fea518d4..ab3b14b5 100644 --- a/Source/EBGeometry_Triangle2DImplem.hpp +++ b/Source/EBGeometry_Triangle2DImplem.hpp @@ -29,11 +29,11 @@ namespace EBGeometry { void Triangle2D::define(const Vec3T& a_normal, const std::array, 3>& a_vertices) noexcept { - int projectDir = 0; + int ignoreDir = 0; for (int dir = 1; dir < 3; dir++) { - if (std::abs(a_normal[dir]) > std::abs(a_normal[projectDir])) { - projectDir = dir; + if (std::abs(a_normal[dir]) > std::abs(a_normal[ignoreDir])) { + ignoreDir = dir; } } @@ -41,7 +41,7 @@ namespace EBGeometry { m_yDir = 0; for (int dir = 0; dir < 3; dir++) { - if (dir != projectDir) { + if (dir != ignoreDir) { m_xDir = std::min(m_xDir, dir); m_yDir = std::max(m_yDir, dir); } @@ -118,7 +118,6 @@ namespace EBGeometry { Vec2T Triangle2D::projectPoint(const Vec3T& a_point) const noexcept { - return Vec2T(a_point[m_xDir], a_point[m_yDir]); } diff --git a/Source/EBGeometry_TriangleImplem.hpp b/Source/EBGeometry_TriangleImplem.hpp index c5f5e875..b6c2194b 100644 --- a/Source/EBGeometry_TriangleImplem.hpp +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -14,6 +14,7 @@ // Our includes #include "EBGeometry_Triangle.hpp" +#include "EBGeometry_Polygon2D.hpp" namespace EBGeometry { @@ -34,9 +35,7 @@ namespace EBGeometry { void Triangle::setVertexPositions(const std::array, 3>& a_vertexPositions) noexcept { - for (int i = 0; i < 3; i++) { - this->m_vertexPositions[i] = a_vertexPositions[i]; - } + m_vertexPositions = a_vertexPositions; this->computeNormal(); @@ -47,18 +46,14 @@ namespace EBGeometry { void Triangle::setVertexNormals(const std::array, 3>& a_vertexNormals) noexcept { - for (int i = 0; i < 3; i++) { - this->m_vertexNormals[i] = a_vertexNormals[i]; - } + m_vertexNormals = a_vertexNormals; } template void Triangle::setEdgeNormals(const std::array, 3>& a_edgeNormals) noexcept { - for (int i = 0; i < 3; i++) { - this->m_edgeNormals[i] = a_edgeNormals[i]; - } + m_edgeNormals = a_edgeNormals; } template @@ -72,10 +67,10 @@ namespace EBGeometry { void Triangle::computeNormal() noexcept { - const Vec3T x1x0 = m_vertexPositions[1] - m_vertexPositions[0]; + const Vec3T x2x0 = m_vertexPositions[2] - m_vertexPositions[0]; const Vec3T x2x1 = m_vertexPositions[2] - m_vertexPositions[1]; - m_triangleNormal = x1x0.cross(x2x1); + m_triangleNormal = x2x1.cross(x2x0); m_triangleNormal = m_triangleNormal / m_triangleNormal.length(); } @@ -152,8 +147,30 @@ namespace EBGeometry { T Triangle::signedDistance(const Vec3T& a_point) const noexcept { - // Check if projection falls inside. - const bool isPointInside = m_triangle2D.isPointInsideWindingNumber(a_point); + // Check if projection falls inside. x is the projected point + const auto x = a_point - m_triangleNormal * (m_triangleNormal.dot(a_point - m_vertexPositions[0])); + +#if 1 + const bool isPointInside = m_triangle2D.isPointInsideCrossingNumber(x); +#else +#warning "Restore original code when debugging is done" + std::vector> vertices; + for (const auto& v : m_vertexPositions) { + vertices.emplace_back(v); + } + + auto tri2 = std::make_shared>(m_triangleNormal, m_vertexPositions); + auto poly2 = std::make_shared>(m_triangleNormal, vertices); + + const bool isPointInsideTri = tri2->isPointInsideCrossingNumber(x); + const bool isPointInsidePoly = poly2->isPointInside(x, Polygon2D::InsideOutsideAlgorithm::CrossingNumber); + + const bool isPointInside = isPointInsidePoly; + + if (isPointInsideTri != isPointInsidePoly) { + std::cout << "should not happen" << std::endl; + } +#endif T ret = std::numeric_limits::max(); @@ -171,7 +188,6 @@ namespace EBGeometry { ret = (delta.length() > std::abs(ret)) ? ret : delta.length() * sgn(m_vertexNormals[i].dot(delta)); } - // Check distances to edges for (int i = 0; i < 3; i++) { const Vec3T& a = m_vertexPositions[i]; const Vec3T& b = m_vertexPositions[(i + 1) % 3]; From d38ccadfe0f2b2997217c3160f81c1e694f590f1 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:09:03 +0200 Subject: [PATCH 12/16] Remove debug code --- Examples/AMReX_PaintEB/main.cpp | 2 +- Source/EBGeometry_TriangleImplem.hpp | 25 ++----------------------- 2 files changed, 3 insertions(+), 24 deletions(-) diff --git a/Examples/AMReX_PaintEB/main.cpp b/Examples/AMReX_PaintEB/main.cpp index 921e7231..38cdf22e 100644 --- a/Examples/AMReX_PaintEB/main.cpp +++ b/Examples/AMReX_PaintEB/main.cpp @@ -37,7 +37,7 @@ class AMReXSDF { // Read in the mesh into a DCEL mesh and partition it into a bounding volume hierarchy auto mesh = EBGeometry::Parser::readIntoDCEL(a_filename); - + // 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++) { diff --git a/Source/EBGeometry_TriangleImplem.hpp b/Source/EBGeometry_TriangleImplem.hpp index b6c2194b..a810c6ab 100644 --- a/Source/EBGeometry_TriangleImplem.hpp +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -14,7 +14,6 @@ // Our includes #include "EBGeometry_Triangle.hpp" -#include "EBGeometry_Polygon2D.hpp" namespace EBGeometry { @@ -147,32 +146,12 @@ namespace EBGeometry { T Triangle::signedDistance(const Vec3T& a_point) const noexcept { + T ret = std::numeric_limits::max(); + // Check if projection falls inside. x is the projected point const auto x = a_point - m_triangleNormal * (m_triangleNormal.dot(a_point - m_vertexPositions[0])); -#if 1 const bool isPointInside = m_triangle2D.isPointInsideCrossingNumber(x); -#else -#warning "Restore original code when debugging is done" - std::vector> vertices; - for (const auto& v : m_vertexPositions) { - vertices.emplace_back(v); - } - - auto tri2 = std::make_shared>(m_triangleNormal, m_vertexPositions); - auto poly2 = std::make_shared>(m_triangleNormal, vertices); - - const bool isPointInsideTri = tri2->isPointInsideCrossingNumber(x); - const bool isPointInsidePoly = poly2->isPointInside(x, Polygon2D::InsideOutsideAlgorithm::CrossingNumber); - - const bool isPointInside = isPointInsidePoly; - - if (isPointInsideTri != isPointInsidePoly) { - std::cout << "should not happen" << std::endl; - } -#endif - - T ret = std::numeric_limits::max(); if (isPointInside) { ret = m_triangleNormal.dot(a_point - m_vertexPositions[0]); From 9ff2e169f0aaaacd3ace6213202f4dc3a3e3ef50 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:13:43 +0200 Subject: [PATCH 13/16] update ci file --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 9ac7a6c3..61cfe527 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -85,7 +85,7 @@ jobs: working-directory: ${{ matrix.directory }} run: | source /opt/intel/oneapi/setvars.sh - icpc -std=c++17 -O0 \ + icx -std=c++17 -O0 \ -Wall -Warray-bounds -Wchar-subscripts -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 \ From 2329866354ce83acb7be3af654e526d742315518 Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:18:24 +0200 Subject: [PATCH 14/16] Update with the goal of choking Intel --- .github/workflows/CI.yml | 2 +- Source/EBGeometry_ParserImplem.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 61cfe527..13103644 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -85,7 +85,7 @@ jobs: working-directory: ${{ matrix.directory }} run: | source /opt/intel/oneapi/setvars.sh - icx -std=c++17 -O0 \ + icpx -std=c++17 -O0 \ -Wall -Warray-bounds -Wchar-subscripts -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 \ diff --git a/Source/EBGeometry_ParserImplem.hpp b/Source/EBGeometry_ParserImplem.hpp index a4bd77bd..da42dc7d 100644 --- a/Source/EBGeometry_ParserImplem.hpp +++ b/Source/EBGeometry_ParserImplem.hpp @@ -169,7 +169,7 @@ Parser::readIntoTriangleBVH(const std::string a_filename) noexcept } template -inline std::vector>> +inline std::vector>> Parser::readIntoTriangleBVH(const std::vector a_files) noexcept { std::vector>> implicitFunctions; From 104bf65c9f266c0c402c21031bba520cb93bdaef Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:31:19 +0200 Subject: [PATCH 15/16] update --- .github/workflows/CI.yml | 2 +- ...EBGeometry_MeshDistanceFunctionsImplem.hpp | 72 +++++++++---------- 2 files changed, 34 insertions(+), 40 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 13103644..c4ba3647 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -86,7 +86,7 @@ jobs: run: | source /opt/intel/oneapi/setvars.sh icpx -std=c++17 -O0 \ - -Wall -Warray-bounds -Wchar-subscripts -Wcomment -Wenum-compare -Wformat -Wuninitialized -Wmaybe-uninitialized -Wmain -Wnarrowing -Wnonnull \ + -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 diff --git a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp index 0c72a10a..151d53e6 100644 --- a/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp +++ b/Source/EBGeometry_MeshDistanceFunctionsImplem.hpp @@ -146,7 +146,7 @@ template T FastMeshSDF::signedDistance(const Vec3T& a_point) const noexcept { - T minDist = std::numeric_limits::infinity(); + T minDist = std::numeric_limits::max(); BVH::Updater updater = [&minDist, &a_point](const std::vector>& faces) noexcept -> void { @@ -157,16 +157,15 @@ FastMeshSDF::signedDistance(const Vec3T& a_point) const noexc } }; - BVH::Visiter visiter = [&minDist, &a_point](const Node& a_node, const T& a_bvDist) noexcept -> bool { + BVH::Visiter visiter = [&minDist](const Node& a_node, const T& a_bvDist) noexcept -> bool { return a_bvDist <= std::abs(minDist); }; BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; @@ -193,7 +192,7 @@ FastMeshSDF::getClosestFaces(const Vec3T& a_point, const bool using BVHMeta = T; // Shortest distance so far. - BVHMeta shortestDistanceSoFar = std::numeric_limits::infinity(); + BVHMeta shortestDistanceSoFar = std::numeric_limits::max(); // Visitation pattern - go into the node if the point is inside or the distance to the BV is shorter than // the shortest distance so far. @@ -204,11 +203,10 @@ FastMeshSDF::getClosestFaces(const Vec3T& a_point, const bool // Sorter for BVH nodes, visit closest nodes first EBGeometry::BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; @@ -280,7 +278,7 @@ template T FastCompactMeshSDF::signedDistance(const Vec3T& a_point) const noexcept { - T minDist = std::numeric_limits::infinity(); + T minDist = std::numeric_limits::max(); BVH::Updater updater = [&minDist, &a_point](const std::vector>& faces) noexcept -> void { @@ -291,16 +289,15 @@ FastCompactMeshSDF::signedDistance(const Vec3T& a_point) cons } }; - BVH::Visiter visiter = [&minDist, &a_point](const Node& a_node, const T& a_bvDist) noexcept -> bool { + BVH::Visiter visiter = [&minDist](const Node& a_node, const T& a_bvDist) noexcept -> bool { return a_bvDist <= std::abs(minDist); }; BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; @@ -326,7 +323,7 @@ FastCompactMeshSDF::getClosestFaces(const Vec3T& a_point, con using BVHMeta = T; // Shortest distance so far. - BVHMeta shortestDistanceSoFar = std::numeric_limits::infinity(); + BVHMeta shortestDistanceSoFar = std::numeric_limits::max(); // Visitation pattern - go into the node if the point is inside or the distance to the BV is shorter than // the shortest distance so far. @@ -337,11 +334,10 @@ FastCompactMeshSDF::getClosestFaces(const Vec3T& a_point, con // Sorter for BVH nodes, visit closest nodes first EBGeometry::BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; @@ -445,7 +441,7 @@ template T FastTriMeshSDF::signedDistance(const Vec3T& a_point) const noexcept { - T minDist = std::numeric_limits::infinity(); + T minDist = std::numeric_limits::max(); BVH::Updater updater = [&minDist, &a_point](const std::vector>& triangles) noexcept -> void { @@ -456,16 +452,15 @@ FastTriMeshSDF::signedDistance(const Vec3T& a_point) const no } }; - BVH::Visiter visiter = [&minDist, &a_point](const Node& a_node, const T& a_bvDist) noexcept -> bool { + BVH::Visiter visiter = [&minDist](const Node& a_node, const T& a_bvDist) noexcept -> bool { return a_bvDist <= std::abs(minDist); }; BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; @@ -492,7 +487,7 @@ FastTriMeshSDF::getClosestTriangles(const Vec3T& a_point, con using BVHMeta = T; // Shortest distance so far. - BVHMeta shortestDistanceSoFar = std::numeric_limits::infinity(); + BVHMeta shortestDistanceSoFar = std::numeric_limits::max(); // Visitation pattern - go into the node if the point is inside or the distance to the BV is shorter than // the shortest distance so far. @@ -503,11 +498,10 @@ FastTriMeshSDF::getClosestTriangles(const Vec3T& a_point, con // Sorter for BVH nodes, visit closest nodes first EBGeometry::BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); }; From 2b4e2a47c39a47ad013b3a1fa13315a52c36e0ee Mon Sep 17 00:00:00 2001 From: Robert Marskar Date: Sun, 24 Aug 2025 18:36:00 +0200 Subject: [PATCH 16/16] Yet again I update --- Examples/EBGeometry_F18/main.cpp | 3 +++ Source/EBGeometry_CSGImplem.hpp | 9 ++++----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/Examples/EBGeometry_F18/main.cpp b/Examples/EBGeometry_F18/main.cpp index f2884e0c..c3f69374 100644 --- a/Examples/EBGeometry_F18/main.cpp +++ b/Examples/EBGeometry_F18/main.cpp @@ -137,6 +137,9 @@ main(int argc, char* argv[]) if (std::abs(sumSlowSlow) - std::abs(sumFastFast) > std::numeric_limits::min()) { std::cerr << "Got wrong distance! Diff = " << std::abs(sumSlowFast) - std::abs(sumFastFast) << "\n"; } + if (std::abs(sumSlowSlow) - std::abs(sumFastSlow) > std::numeric_limits::min()) { + std::cerr << "Got wrong distance! Diff = " << std::abs(sumSlowFast) - std::abs(sumFastSlow) << "\n"; + } const std::chrono::duration slowSlowTime = (t1 - t0); const std::chrono::duration slowFastTime = (t2 - t1); diff --git a/Source/EBGeometry_CSGImplem.hpp b/Source/EBGeometry_CSGImplem.hpp index 1f0a679c..1d61e26c 100644 --- a/Source/EBGeometry_CSGImplem.hpp +++ b/Source/EBGeometry_CSGImplem.hpp @@ -330,11 +330,10 @@ FastUnionIF::value(const Vec3T& a_point) const noexcept }; BVH::Sorter sorter = - [&a_point](std::array, T>, K>& a_leaves) noexcept -> void { - std::sort( - a_leaves.begin(), - a_leaves.end(), - [&a_point](const std::pair, T>& n1, + [](std::array, T>, K>& a_leaves) noexcept -> void { + std::sort(a_leaves.begin(), + a_leaves.end(), + [](const std::pair, T>& n1, const std::pair, T>& n2) -> bool { return n1.second > n2.second; }); };