diff --git a/EBGeometry.hpp b/EBGeometry.hpp index 32e0494d..f0b984ef 100644 --- a/EBGeometry.hpp +++ b/EBGeometry.hpp @@ -15,7 +15,6 @@ #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/Examples/AMReX_Shapes/main.cpp b/Examples/AMReX_Shapes/main.cpp index 63b1182e..ea971374 100644 --- a/Examples/AMReX_Shapes/main.cpp +++ b/Examples/AMReX_Shapes/main.cpp @@ -169,6 +169,11 @@ main(int argc, char* argv[]) func = std::make_shared>(0.5, 2.0 * Vec3::one(), 0.5, 4); } + else if (whichGeom == 14) { // Rounded cylinder + rb = RealBox({-1, -1, -1}, {1, 1, 1}); + + func = std::make_shared>(0.5, 0.1, 1.0); + } // AMReX uses the opposite sign. func = EBGeometry::Complement(func); diff --git a/Source/EBGeometry_AnalyticDistanceFunctions.hpp b/Source/EBGeometry_AnalyticDistanceFunctions.hpp index 3c24e2d0..d8fe7739 100644 --- a/Source/EBGeometry_AnalyticDistanceFunctions.hpp +++ b/Source/EBGeometry_AnalyticDistanceFunctions.hpp @@ -1091,6 +1091,67 @@ class PerlinSDF : public SignedDistanceFunction }; }; +/*! + @brief Rounded cylinder signed-distance function +*/ +template +class RoundedCylinderSDF : public SignedDistanceFunction +{ +public: + /*! + @brief Disallowed default constructor + */ + RoundedCylinderSDF() = delete; + + /*! + @brief Full constructor. User inputs the radius, curvature, and height. Adjust parameters as necessary. + @param[in] a_radius Cylinder outer radius + @param[in] a_curvature Corner curvature + @param[in] a_height Cylinder height + */ + RoundedCylinderSDF(const T a_radius, const T a_curvature, const T a_height) noexcept + { + m_majorRadius = a_radius - a_curvature; + m_minorRadius = a_curvature; + m_height = 0.5 * a_height - a_curvature; + } + + /*! + @brief Destructor (does nothing). + */ + virtual ~RoundedCylinderSDF() noexcept + {} + + /*! + @brief Signed distance function for the rounded box + */ + virtual T + signedDistance(const Vec3T& a_point) const noexcept override + { + const auto xz = sqrt(a_point[2] * a_point[2] + a_point[0] * a_point[0]); + const auto d1 = Vec2T(xz - 2.0 * m_majorRadius + m_minorRadius, std::abs(a_point[1]) - m_height); + const auto d2 = Vec2T(std::max(d1.x, 0.0), std::max(d1.y, 0.0)); + + return std::min(std::max(d1.x, d1.y), 0.0) + sqrt(d2.x * d2.x + d2.y * d2.y) - m_minorRadius; + } + +protected: + /*! + @brief Major radius + */ + T m_majorRadius; + + /*! + @brief Minor radius + */ + T m_minorRadius; + + /*! + @brief Minor radius + */ + T m_height; +}; + #include "EBGeometry_NamespaceFooter.hpp" #endif diff --git a/Source/EBGeometry_Triangle.hpp b/Source/EBGeometry_Triangle.hpp index efb71f0b..36103da4 100644 --- a/Source/EBGeometry_Triangle.hpp +++ b/Source/EBGeometry_Triangle.hpp @@ -17,7 +17,6 @@ // Our includes #include "EBGeometry_Vec.hpp" -#include "EBGeometry_Triangle2D.hpp" namespace EBGeometry { @@ -228,32 +227,16 @@ namespace EBGeometry { @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 diff --git a/Source/EBGeometry_Triangle2D.hpp b/Source/EBGeometry_Triangle2D.hpp deleted file mode 100644 index 83c75501..00000000 --- a/Source/EBGeometry_Triangle2D.hpp +++ /dev/null @@ -1,204 +0,0 @@ -/* 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 deleted file mode 100644 index ab3b14b5..00000000 --- a/Source/EBGeometry_Triangle2DImplem.hpp +++ /dev/null @@ -1,214 +0,0 @@ -/* 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 ignoreDir = 0; - - for (int dir = 1; dir < 3; dir++) { - if (std::abs(a_normal[dir]) > std::abs(a_normal[ignoreDir])) { - ignoreDir = dir; - } - } - - m_xDir = 3; - m_yDir = 0; - - for (int dir = 0; dir < 3; dir++) { - if (dir != ignoreDir) { - 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 index a810c6ab..f8dd2e82 100644 --- a/Source/EBGeometry_TriangleImplem.hpp +++ b/Source/EBGeometry_TriangleImplem.hpp @@ -37,8 +37,6 @@ namespace EBGeometry { m_vertexPositions = a_vertexPositions; this->computeNormal(); - - m_triangle2D.define(m_triangleNormal, m_vertexPositions); } template @@ -146,54 +144,45 @@ namespace EBGeometry { T Triangle::signedDistance(const Vec3T& a_point) const noexcept { + // Here is a message from the past: If one wants, one can precompute v21, v32, v13 + // as well as many other quantities (e.g., v21.cross(m_triangleNormal)). This might + // be helpful in order to speed things up a little bit. 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])); - - const bool isPointInside = m_triangle2D.isPointInsideCrossingNumber(x); - - 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]; + auto sgn = [](const T x) -> int { return (x > 0.0) ? 1 : -1; }; - ret = (delta.length() > std::abs(ret)) ? ret : delta.length() * sgn(m_vertexNormals[i].dot(delta)); - } + const Vec3 v21 = m_vertexPositions[1] - m_vertexPositions[0]; + const Vec3 v32 = m_vertexPositions[2] - m_vertexPositions[1]; + const Vec3 v13 = m_vertexPositions[0] - m_vertexPositions[2]; - for (int i = 0; i < 3; i++) { - const Vec3T& a = m_vertexPositions[i]; - const Vec3T& b = m_vertexPositions[(i + 1) % 3]; + const Vec3 p1 = a_point - m_vertexPositions[0]; + const Vec3 p2 = a_point - m_vertexPositions[1]; + const Vec3 p3 = a_point - m_vertexPositions[2]; - const T t = this->projectPointToEdge(a_point, a, b); + const T s0 = sgn(dot(v21.cross(m_triangleNormal), p1)); + const T s1 = sgn(dot(v32.cross(m_triangleNormal), p2)); + const T s2 = sgn(dot(v13.cross(m_triangleNormal), p3)); - if (t > 0.0 && t < 1.0) { - const Vec3T delta = a_point - (a + t * (b - a)); + const T t1 = dot(p1, v21) / dot(v21, v21); + const T t2 = dot(p2, v32) / dot(v32, v32); + const T t3 = dot(p3, v13) / dot(v13, v13); - ret = (delta.length() > std::abs(ret)) ? ret : delta.length() * sgn(m_edgeNormals[i].dot(delta)); - } - } - } + const Vec3 y1 = p1 - t1 * v21; + const Vec3 y2 = p2 - t2 * v32; + const Vec3 y3 = p3 - t3 * v13; - return ret; - } + // Distance to vertices + ret = (p1.length() > std::abs(ret)) ? ret : p1.length() * sgn(m_vertexNormals[0].dot(p1)); + ret = (p2.length() > std::abs(ret)) ? ret : p2.length() * sgn(m_vertexNormals[1].dot(p2)); + ret = (p3.length() > std::abs(ret)) ? ret : p3.length() * sgn(m_vertexNormals[2].dot(p3)); - 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; + // Distance to edges + ret = (t1 > 0.0 && t1 < 1.0 && y1.length() < std::abs(ret)) ? y1.length() * sgn(m_edgeNormals[0].dot(y1)) : ret; + ret = (t2 > 0.0 && t2 < 1.0 && y2.length() < std::abs(ret)) ? y2.length() * sgn(m_edgeNormals[1].dot(y2)) : ret; + ret = (t3 > 0.0 && t3 < 1.0 && y3.length() < std::abs(ret)) ? y3.length() * sgn(m_edgeNormals[2].dot(y3)) : ret; - return a.dot(b) / (b.dot(b)); + // Note that s0 + s1 + s2 >= 2.0 is a point-in-polygon test. + return (s0 + s1 + s2 >= 2.0) ? m_triangleNormal.dot(p1) : ret; } } // namespace EBGeometry