From b2da4ca1adbf6ed35d35e939b0661449827c773d Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Mon, 6 May 2024 15:49:16 -0700 Subject: [PATCH 01/11] Fix bounding ball calculation for coaxial lumped ports --- palace/utils/geodata.cpp | 196 +++++++++++++++++++-------------------- 1 file changed, 94 insertions(+), 102 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index f96064728..66d62a54c 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -744,10 +744,10 @@ bool EigenLE(const Eigen::Vector3d &x, const Eigen::Vector3d &y) int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, std::vector &vertices) { - std::set vertex_indices; if (!mesh.GetNodes()) { // Linear mesh, work with element vertices directly. + std::set vertex_indices; mfem::Array v; if (bdr) { @@ -884,6 +884,19 @@ int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array &m return dominant_rank; } +// Compute the distance from a point orthogonal to the list of normal axes, relative to +// the given origin. +auto PerpendicularDistance(const std::initializer_list &normals, + const Eigen::Vector3d &origin, const Eigen::Vector3d &v) +{ + Eigen::Vector3d v0 = v - origin; + for (const auto &n : normals) + { + v0 -= n.dot(v0) * n; + } + return v0.norm(); +}; + // Calculates a bounding box from a point cloud, result is broadcast across all processes. BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, const std::vector &vertices, @@ -910,57 +923,53 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, MFEM_VERIFY(std::max_element(vertices.begin(), vertices.end(), DistFromP_000) == p_111, "p_000 and p_111 must be mutually opposing points!"); - // Define a diagonal of the ASSUMED cuboid bounding box. Store references as this is - // useful for checking pointers later. + // Define a diagonal of the ASSUMED cuboid bounding box. const auto &v_000 = *p_000; const auto &v_111 = *p_111; MFEM_VERIFY(&v_000 != &v_111, "Minimum and maximum extents cannot be identical!"); - const auto origin = v_000; + const Eigen::Vector3d origin = v_000; const Eigen::Vector3d n_1 = (v_111 - v_000).normalized(); - // Compute the distance from the normal axis. Note: everything has been oriented - // relative to v_000 == (0,0,0). - auto PerpendicularDistance = [&n_1, &origin](const Eigen::Vector3d &v) - { return ((v - origin) - (v - origin).dot(n_1) * n_1).norm(); }; - // Find the vertex furthest from the diagonal axis. We cannot know yet if this defines // (001) or (011). - const auto &t_0 = - *std::max_element(vertices.begin(), vertices.end(), - [PerpendicularDistance](const auto &x, const auto &y) - { return PerpendicularDistance(x) < PerpendicularDistance(y); }); - MFEM_VERIFY(&t_0 != &v_000, "Vertices are degenerate!"); - MFEM_VERIFY(&t_0 != &v_111, "Vertices are degenerate!"); - - // Use the discovered vertex to define a second direction and thus a plane. + const auto &t_0 = *std::max_element(vertices.begin(), vertices.end(), + [&](const auto &x, const auto &y) + { + return PerpendicularDistance({n_1}, origin, x) < + PerpendicularDistance({n_1}, origin, y); + }); + MFEM_VERIFY(&t_0 != &v_000 && &t_0 != &v_111, "Vertices are degenerate!"); + + // Use the discovered vertex to define a second direction and thus a plane. n_1 and n_2 + // now define a planar coordinate system intersecting the main diagonal, and two + // opposite edges of the cuboid. Now look for a component that maximizes distance from + // the planar system: complete the axes with a cross, then use a dot product to pick + // the greatest deviation. const Eigen::Vector3d n_2 = ((t_0 - origin) - (t_0 - origin).dot(n_1) * n_1).normalized(); - // n_1 and n_2 now define a planar coordinate system intersecting the main diagonal, and - // two opposite edges of the cuboid. Now look for a component that maximizes distance - // from the planar system: complete the axes with a cross, then use a dot product to - // pick the greatest deviation. - auto OutOfPlaneDistance = [&n_1, &n_2, &origin](const Eigen::Vector3d &v) - { - return ((v - origin) - (v - origin).dot(n_1) * n_1 - (v - origin).dot(n_2) * n_2) - .norm(); - }; - // Collect the furthest point from the plane. - auto max_distance = OutOfPlaneDistance(*std::max_element( - vertices.begin(), vertices.end(), [OutOfPlaneDistance](const auto &x, const auto &y) - { return OutOfPlaneDistance(x) < OutOfPlaneDistance(y); })); - + auto max_distance = PerpendicularDistance( + {n_1, n_2}, origin, + *std::max_element(vertices.begin(), vertices.end(), + [&](const auto &x, const auto &y) + { + return PerpendicularDistance({n_1, n_2}, origin, x) < + PerpendicularDistance({n_1, n_2}, origin, y); + })); constexpr double rel_tol = 1e-6; - box.planar = max_distance < (rel_tol * (v_111 - v_000).norm()); + box.planar = (max_distance < (rel_tol * (v_111 - v_000).norm())); // Given numerical tolerance, collect other points with an almost matching distance. + const double cooincident_tol = rel_tol * max_distance; std::vector vertices_out_of_plane; - const double cooincident_tolerance = rel_tol * max_distance; - std::copy_if( - vertices.begin(), vertices.end(), std::back_inserter(vertices_out_of_plane), - [OutOfPlaneDistance, cooincident_tolerance, max_distance](const auto &v) - { return std::abs(OutOfPlaneDistance(v) - max_distance) < cooincident_tolerance; }); + std::copy_if(vertices.begin(), vertices.end(), + std::back_inserter(vertices_out_of_plane), + [&](const auto &v) + { + return std::abs(PerpendicularDistance({n_1, n_2}, origin, v) - + max_distance) < cooincident_tol; + }); // Given candidates t_0 and t_1, the closer to origin defines v_001. const auto &t_1 = box.planar @@ -1001,80 +1010,63 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, BoundingBall ball; if (dominant_rank == Mpi::Rank(comm)) { - // Pick a candidate 000 vertex using lexicographic sort. This can be vulnerable to - // floating point precision if there is no directly opposed vertex. - // Pick candidate 111 as the furthest from this candidate, then reassign 000 as the - // furthest from 111. Such a pair has to form the diagonal for a point cloud defining a - // ball. Verify that p_111 is also the maximum distance from p_000 -> a diagonal is - // found. - MFEM_VERIFY(vertices.size() >= 3, - "A bounding ball requires a minimum of three vertices for this algorithm!"); - auto p_000 = std::min_element(vertices.begin(), vertices.end(), EigenLE); - auto DistFromP_000 = [&p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - *p_000).norm() < (y - *p_000).norm(); }; - auto p_111 = std::max_element(vertices.begin(), vertices.end(), DistFromP_000); - auto DistFromP_111 = [&p_111](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - *p_111).norm() < (y - *p_111).norm(); }; - p_000 = std::max_element(vertices.begin(), vertices.end(), DistFromP_111); - MFEM_VERIFY(std::max_element(vertices.begin(), vertices.end(), DistFromP_000) == p_111, - "p_000 and p_111 must be mutually opposing points!"); - - const auto &min = *p_000; - const auto &max = *p_111; - Eigen::Vector3d delta = max - min; - ball.radius = 0.5 * delta.norm(); - Vector3dMap(ball.center.data()) = 0.5 * (min + max); - - // Project onto this candidate diameter, and pick a vertex furthest away. Check that - // this resulting distance is less than or equal to the radius, and use the resulting - // direction to compute another in plane vector. Assumes all delta are normalized, and - // applies a common origin as part of the projection. - auto PerpendicularDistance = [min](const std::initializer_list &deltas, - const Eigen::Vector3d &vin) - { - Eigen::Vector3d v = vin - min; - for (const auto &d : deltas) + // Pick the ball center as the centroid of all of the points in the cloud. Then, pick + // the radial direction by finding the point furthest away from the center. + const auto origin = [&]() + { + Eigen::Vector3d v = Eigen::Vector3d::Zero(); + for (const auto &p : vertices) { - v -= d.dot(v) * d; + v += p; } - return v.norm(); - }; + v /= vertices.size(); + return v; + }(); + auto DistFromOrigin = [&origin](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - origin).norm() < (y - origin).norm(); }; + const auto &t_0 = *std::max_element(vertices.begin(), vertices.end(), DistFromOrigin); + const Eigen::Vector3d n_1 = (t_0 - origin).normalized(); + + // Initialize the bounding ball data. + Vector3dMap(ball.center.data()) = origin; + ball.radius = (t_0 - origin).norm(); - delta.normalize(); - const auto &perp = *std::max_element( - vertices.begin(), vertices.end(), - [&delta, PerpendicularDistance](const auto &x, const auto &y) - { return PerpendicularDistance({delta}, x) < PerpendicularDistance({delta}, y); }); + // Project onto this candidate diameter, and pick a vertex furthest away. Check that + // this resulting distance is less than or equal to the radius. + const auto &t_1 = *std::max_element(vertices.begin(), vertices.end(), + [&](const auto &x, const auto &y) + { + return PerpendicularDistance({n_1}, origin, x) < + PerpendicularDistance({n_1}, origin, y); + }); + MFEM_VERIFY(&t_1 != &t_0, "Vertices are degenerate!"); constexpr double rel_tol = 1.0e-6; - MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <= - rel_tol * ball.radius, + MFEM_VERIFY(PerpendicularDistance({n_1}, origin, t_1) < (1.0 + rel_tol) * ball.radius, "Furthest point perpendicular must be on the exterior of the ball: " - << PerpendicularDistance({delta}, perp) << " vs. " << ball.radius + << PerpendicularDistance({n_1}, origin, t_1) << " vs. " << ball.radius << "!"); - // Compute a perpendicular to the circle using the cross product. - const Eigen::Vector3d n_radial = (perp - CVector3dMap(ball.center.data())).normalized(); - Vector3dMap(ball.planar_normal.data()) = delta.cross(n_radial).normalized(); + // Use the resulting point to compute another in-plane vector. This vector may not be + // perfectly orthogonal to n_1, so we explicitly orthonormalize it even though it does + // not have to be. + const Eigen::Vector3d n_2 = + ((t_1 - origin) - (t_1 - origin).dot(n_1) * n_1).normalized(); // Compute the point furthest out of the plane discovered. If below tolerance, this - // means the ball is 2D. - const auto &out_of_plane = *std::max_element( - vertices.begin(), vertices.end(), - [&delta, &n_radial, PerpendicularDistance](const auto &x, const auto &y) - { - return PerpendicularDistance({delta, n_radial}, x) < - PerpendicularDistance({delta, n_radial}, y); - }); - - ball.planar = - PerpendicularDistance({delta, n_radial}, out_of_plane) / ball.radius < rel_tol; - if (!ball.planar) - { - // The points are not functionally coplanar, zero out the normal. - MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <= - rel_tol * ball.radius, - "Furthest point perpendicular must be on the exterior of the sphere!"); - Vector3dMap(ball.planar_normal.data()) *= 0; + // means the ball is 2D. For the planar case, compute a perpendicular to the circle + // using the cross product. + auto max_distance = PerpendicularDistance( + {n_1, n_2}, origin, + *std::max_element(vertices.begin(), vertices.end(), + [&](const auto &x, const auto &y) + { + return PerpendicularDistance({n_1, n_2}, origin, x) < + PerpendicularDistance({n_1, n_2}, origin, y); + })); + ball.planar = (max_distance < rel_tol * ball.radius); + if (ball.planar) + { + Vector3dMap(ball.planar_normal.data()) = n_1.cross(n_2).normalized(); } } From 7ba4f14589c4321e798a778340ec0a223bf67726 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 7 May 2024 10:51:20 -0700 Subject: [PATCH 02/11] Remove BoundingBall and replace use in coaxial lumped ports with BoundingBox --- palace/fem/lumpedelement.cpp | 52 ++++++++++------- palace/fem/lumpedelement.hpp | 20 +++---- palace/utils/geodata.cpp | 105 ++--------------------------------- palace/utils/geodata.hpp | 36 ++---------- 4 files changed, 51 insertions(+), 162 deletions(-) diff --git a/palace/fem/lumpedelement.cpp b/palace/fem/lumpedelement.cpp index fd4b3b550..52a708571 100644 --- a/palace/fem/lumpedelement.cpp +++ b/palace/fem/lumpedelement.cpp @@ -7,6 +7,7 @@ #include "fem/integrator.hpp" #include "linalg/vector.hpp" #include "utils/communication.hpp" +#include "utils/geodata.hpp" namespace palace { @@ -36,16 +37,18 @@ UniformElementData::UniformElementData(const std::array &input_dir, const mfem::ParMesh &mesh = *fespace.GetParMesh(); int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); - bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true); + auto bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true); + MFEM_VERIFY(bounding_box.planar, + "Boundary elements must be coplanar to define a rectangular lumped element!"); // Check that the bounding box discovered matches the area. This validates that the // boundary elements form a right angled quadrilateral port. constexpr double rel_tol = 1.0e-6; double A = GetArea(fespace, attr_marker); - MFEM_VERIFY((!bounding_box.planar || (std::abs(A - bounding_box.Area()) / A < rel_tol)), + MFEM_VERIFY(std::abs(A - bounding_box.Area()) < rel_tol * A, "Discovered bounding box area " << bounding_box.Area() << " and integrated area " << A - << " do not match: Planar port geometry is not a quadrilateral!"); + << " do not match, planar port geometry is not a quadrilateral!"); // Check the user specified direction aligns with an axis direction. constexpr double angle_warning_deg = 0.1; @@ -79,9 +82,9 @@ UniformElementData::UniformElementData(const std::array &input_dir, } MFEM_VERIFY(std::any_of(deviation_deg.begin(), deviation_deg.end(), [](double x) { return x < angle_error_deg; }), - "Specified direction does not align sufficiently with bounding box axes: " - << deviation_deg[0] << ' ' << deviation_deg[1] << ' ' << deviation_deg[2] - << " tolerance " << angle_error_deg << '!'); + "Specified direction does not align sufficiently with bounding box axes (" + << deviation_deg[0] << ", " << deviation_deg[1] << ", " + << deviation_deg[2] << " vs. tolerance " << angle_error_deg << ")!"); direction.SetSize(input_dir.size()); std::copy(input_dir.begin(), input_dir.end(), direction.begin()); direction /= direction.Norml2(); @@ -104,7 +107,7 @@ UniformElementData::GetModeCoefficient(double coef) const attr_list, source); } -CoaxialElementData::CoaxialElementData(const std::array &direction, +CoaxialElementData::CoaxialElementData(const std::array &input_dir, const mfem::Array &attr_list, mfem::ParFiniteElementSpace &fespace) : LumpedElementData(attr_list) @@ -112,26 +115,37 @@ CoaxialElementData::CoaxialElementData(const std::array &direction, const mfem::ParMesh &mesh = *fespace.GetParMesh(); int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); - bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true); - MFEM_VERIFY(bounding_ball.planar, + auto bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true); + MFEM_VERIFY(bounding_box.planar, "Boundary elements must be coplanar to define a coaxial lumped element!"); - sign = (direction[0] > 0); - // Get inner radius of annulus assuming full 2π circumference. + // Direction of the excitation as +/-r̂. + direction = (input_dir[0] > 0); + origin.SetSize(bounding_box.center.size()); + std::copy(bounding_box.center.begin(), bounding_box.center.end(), origin.begin()); + + // Get inner radius of annulus assuming full 2π circumference. Coarse tolerance to allow + // for approximately-meshed circles. + constexpr double rel_tol = 1.0e-3; double A = GetArea(fespace, attr_marker); - MFEM_VERIFY(bounding_ball.radius > 0.0 && - std::pow(bounding_ball.radius, 2) - A / M_PI > 0.0, - "Coaxial element boundary is not defined correctly: Radius " - << bounding_ball.radius << ", area " << A << "!"); - ra = std::sqrt(std::pow(bounding_ball.radius, 2) - A / M_PI); + + auto lengths = bounding_box.Lengths(); + MFEM_VERIFY((1.0 - rel_tol) * lengths[1] < lengths[0] && + lengths[0] < (1.0 + rel_tol) * lengths[1], + "Bounding box for coaxial lumped port elements is not approximately square (" + << "sides = " << lengths[0] << ", " << lengths[1] << ")!"); + r_outer = 0.25 * (lengths[0] + lengths[1]); + MFEM_VERIFY(std::pow(r_outer, 2) - A / M_PI > 0.0, + "Coaxial element boundary is not defined correctly (radius " + << r_outer << ", area " << A << ")!"); + r_inner = std::sqrt(std::pow(r_outer, 2) - A / M_PI); } std::unique_ptr CoaxialElementData::GetModeCoefficient(double coef) const { - coef *= (sign ? 1.0 : -1.0); - mfem::Vector x0(bounding_ball.center.size()); - std::copy(bounding_ball.center.begin(), bounding_ball.center.end(), x0.begin()); + coef *= direction; + mfem::Vector x0(origin); auto Source = [coef, x0](const mfem::Vector &x, mfem::Vector &f) -> void { f = x; diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index f154e058c..f6e3b18ef 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -6,7 +6,6 @@ #include #include -#include "utils/geodata.hpp" namespace palace { @@ -37,9 +36,6 @@ class LumpedElementData class UniformElementData : public LumpedElementData { private: - // Bounding box defining the rectangular lumped port. - mesh::BoundingBox bounding_box; - // Cartesian vector specifying signed direction of incident field. mfem::Vector direction; @@ -61,21 +57,21 @@ class UniformElementData : public LumpedElementData class CoaxialElementData : public LumpedElementData { private: - // Bounding ball defined by boundary element. - mesh::BoundingBall bounding_ball; + // Sign of incident field, +1 for +r̂, -1 for -r̂. + double direction; - // Sign of incident field, +r̂ if true. - bool sign; + // Origin of the coaxial annulus. + mfem::Vector origin; - // Inner radius of coaxial annulus. - double ra; + // Outer and inner radii of coaxial annulus. + double r_outer, r_inner; public: - CoaxialElementData(const std::array &direction, + CoaxialElementData(const std::array &input_dir, const mfem::Array &attr_list, mfem::ParFiniteElementSpace &fespace); - double GetGeometryLength() const override { return std::log(bounding_ball.radius / ra); } + double GetGeometryLength() const override { return std::log(r_outer / r_inner); } double GetGeometryWidth() const override { return 2.0 * M_PI; } std::unique_ptr diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 66d62a54c..90baa7f2c 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1002,83 +1002,6 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, return box; } -// Calculates a bounding ball from a point cloud, result is broadcast across all processes. -BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, - const std::vector &vertices, - int dominant_rank) -{ - BoundingBall ball; - if (dominant_rank == Mpi::Rank(comm)) - { - // Pick the ball center as the centroid of all of the points in the cloud. Then, pick - // the radial direction by finding the point furthest away from the center. - const auto origin = [&]() - { - Eigen::Vector3d v = Eigen::Vector3d::Zero(); - for (const auto &p : vertices) - { - v += p; - } - v /= vertices.size(); - return v; - }(); - auto DistFromOrigin = [&origin](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - origin).norm() < (y - origin).norm(); }; - const auto &t_0 = *std::max_element(vertices.begin(), vertices.end(), DistFromOrigin); - const Eigen::Vector3d n_1 = (t_0 - origin).normalized(); - - // Initialize the bounding ball data. - Vector3dMap(ball.center.data()) = origin; - ball.radius = (t_0 - origin).norm(); - - // Project onto this candidate diameter, and pick a vertex furthest away. Check that - // this resulting distance is less than or equal to the radius. - const auto &t_1 = *std::max_element(vertices.begin(), vertices.end(), - [&](const auto &x, const auto &y) - { - return PerpendicularDistance({n_1}, origin, x) < - PerpendicularDistance({n_1}, origin, y); - }); - MFEM_VERIFY(&t_1 != &t_0, "Vertices are degenerate!"); - constexpr double rel_tol = 1.0e-6; - MFEM_VERIFY(PerpendicularDistance({n_1}, origin, t_1) < (1.0 + rel_tol) * ball.radius, - "Furthest point perpendicular must be on the exterior of the ball: " - << PerpendicularDistance({n_1}, origin, t_1) << " vs. " << ball.radius - << "!"); - - // Use the resulting point to compute another in-plane vector. This vector may not be - // perfectly orthogonal to n_1, so we explicitly orthonormalize it even though it does - // not have to be. - const Eigen::Vector3d n_2 = - ((t_1 - origin) - (t_1 - origin).dot(n_1) * n_1).normalized(); - - // Compute the point furthest out of the plane discovered. If below tolerance, this - // means the ball is 2D. For the planar case, compute a perpendicular to the circle - // using the cross product. - auto max_distance = PerpendicularDistance( - {n_1, n_2}, origin, - *std::max_element(vertices.begin(), vertices.end(), - [&](const auto &x, const auto &y) - { - return PerpendicularDistance({n_1, n_2}, origin, x) < - PerpendicularDistance({n_1, n_2}, origin, y); - })); - ball.planar = (max_distance < rel_tol * ball.radius); - if (ball.planar) - { - Vector3dMap(ball.planar_normal.data()) = n_1.cross(n_2).normalized(); - } - } - - // Broadcast result to all processors. - Mpi::Broadcast(3, ball.center.data(), dominant_rank, comm); - Mpi::Broadcast(3, ball.planar_normal.data(), dominant_rank, comm); - Mpi::Broadcast(1, &ball.radius, dominant_rank, comm); - Mpi::Broadcast(1, &ball.planar, dominant_rank, comm); - - return ball; -} - double LengthFromPointCloud(MPI_Comm comm, const std::vector &vertices, int dominant_rank, const std::array &dir) { @@ -1100,23 +1023,6 @@ double LengthFromPointCloud(MPI_Comm comm, const std::vector &v } // namespace -double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &marker, - bool bdr, const std::array &dir) -{ - std::vector vertices; - int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices); - return LengthFromPointCloud(mesh.GetComm(), vertices, dominant_rank, dir); -} - -double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, - const std::array &dir) -{ - mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); - marker = 0; - marker[attr - 1] = 1; - return GetProjectedLength(mesh, marker, bdr, dir); -} - BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr) { @@ -1133,20 +1039,21 @@ BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr) return GetBoundingBox(mesh, marker, bdr); } -BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &marker, - bool bdr) +double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &marker, + bool bdr, const std::array &dir) { std::vector vertices; int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices); - return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank); + return LengthFromPointCloud(mesh.GetComm(), vertices, dominant_rank, dir); } -BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr) +double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, + const std::array &dir) { mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); marker = 0; marker[attr - 1] = 1; - return GetBoundingBall(mesh, marker, bdr); + return GetProjectedLength(mesh, marker, bdr, dir); } mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array &marker, diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 28723cff4..2aced85e2 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -67,7 +67,8 @@ mfem::Array AttrToMarker(int max_attr, const T &attr_list, bool skip_invali return marker; } -// Helper function to construct the bounding box for all elements with the given attribute. +// Helper function to construct the axis-aligned bounding box for all elements with the +// given attribute. void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, mfem::Vector &min, mfem::Vector &max); void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr, @@ -101,31 +102,8 @@ struct BoundingBox std::array Deviation(const std::array &direction) const; }; -// Struct describing a bounding ball in terms of a center and radius. If a ball is two -// dimensional, additionally provides a normal to the plane. -struct BoundingBall -{ - // The centroid of the ball. - std::array center; - - // The radius of the ball from the center. - double radius; - - // If the ball is two dimensional, the normal defining the planar surface. Zero magnitude - // if a sphere. - std::array planar_normal; - - // Whether or not this bounding ball is two dimensional. - bool planar; - - // Compute the area of the bounding box spanned by the first two normals. - double Area() const { return M_PI * std::pow(radius, 2.0); } - - // Compute the volume of a 3D bounding box. Returns zero if planar. - double Volume() const { return planar ? 0.0 : (4 * M_PI / 3) * std::pow(radius, 3.0); } -}; - -// Helper functions for computing bounding boxes from a mesh and markers. +// Helper functions for computing bounding boxes from a mesh and markers. These do not need +// to be axis-aligned. BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr); BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr); @@ -136,12 +114,6 @@ double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &mar double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, const std::array &dir); -// Given a mesh and a marker, compute the diameter of a bounding circle/sphere, assuming -// that the extrema points are in the marked group. -BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &marker, - bool bdr); -BoundingBall GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr); - // Helper function to compute the average surface normal for all elements with the given // attribute. mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array &marker, From c22bf04f2197c6f3058a1d0ccb1bd48cb7a53bc9 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 7 May 2024 15:55:05 -0700 Subject: [PATCH 03/11] Fix bounding box for non-rectangular boundaries --- palace/fem/lumpedelement.cpp | 16 ++++++++-------- palace/utils/geodata.cpp | 2 +- palace/utils/geodata.hpp | 13 +++++++++---- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/palace/fem/lumpedelement.cpp b/palace/fem/lumpedelement.cpp index 52a708571..5b5111d30 100644 --- a/palace/fem/lumpedelement.cpp +++ b/palace/fem/lumpedelement.cpp @@ -38,14 +38,14 @@ UniformElementData::UniformElementData(const std::array &input_dir, int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); auto bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true); - MFEM_VERIFY(bounding_box.planar, - "Boundary elements must be coplanar to define a rectangular lumped element!"); // Check that the bounding box discovered matches the area. This validates that the - // boundary elements form a right angled quadrilateral port. + // boundary elements form a right angled quadrilateral port. Rectangular elements are + // allowed to be non-planar, for example a union of 2D quadrilaterals in 3D space (which + // can be "unfolded" to be a planar rectangle). constexpr double rel_tol = 1.0e-6; double A = GetArea(fespace, attr_marker); - MFEM_VERIFY(std::abs(A - bounding_box.Area()) < rel_tol * A, + MFEM_VERIFY(!bounding_box.planar || (std::abs(A - bounding_box.Area()) < rel_tol * A), "Discovered bounding box area " << bounding_box.Area() << " and integrated area " << A << " do not match, planar port geometry is not a quadrilateral!"); @@ -125,16 +125,16 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, std::copy(bounding_box.center.begin(), bounding_box.center.end(), origin.begin()); // Get inner radius of annulus assuming full 2π circumference. Coarse tolerance to allow - // for approximately-meshed circles. - constexpr double rel_tol = 1.0e-3; + // for approximately-meshed circles. The returned bounding box is not a bounding box, it + // is inscribed in the outer circle. + constexpr double rel_tol = 1.0e-2; double A = GetArea(fespace, attr_marker); - auto lengths = bounding_box.Lengths(); MFEM_VERIFY((1.0 - rel_tol) * lengths[1] < lengths[0] && lengths[0] < (1.0 + rel_tol) * lengths[1], "Bounding box for coaxial lumped port elements is not approximately square (" << "sides = " << lengths[0] << ", " << lengths[1] << ")!"); - r_outer = 0.25 * (lengths[0] + lengths[1]); + r_outer = 0.25 * std::sqrt(2.0) * (lengths[0] + lengths[1]); MFEM_VERIFY(std::pow(r_outer, 2) - A / M_PI > 0.0, "Coaxial element boundary is not defined correctly (radius " << r_outer << ", area " << A << ")!"); diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 90baa7f2c..fcfa2192a 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -781,7 +781,7 @@ int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array &m } else { - // Nonlinear mesh, need to process point matrices. + // Curved mesh, need to process point matrices. const int ref = mesh.GetNodes()->FESpace()->GetMaxElementOrder(); mfem::DenseMatrix pointmat; // 3 x N mfem::IsoparametricTransformation T; diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 2aced85e2..3b3bba66a 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -92,18 +92,23 @@ struct BoundingBox // Compute the area of the bounding box spanned by the first two normals. double Area() const; - // Compute the volume of a 3D bounding box. Returns zero if planar. + // Compute the volume of the 3D bounding box. Returns zero if planar. double Volume() const; - // Compute the lengths of each axis. + // Compute the lengths along each axis. std::array Lengths() const; - // Compute the deviation in degrees of a vector from each of the normal directions. + // Compute the deviation in degrees of a vector from each of the axis directions. std::array Deviation(const std::array &direction) const; }; // Helper functions for computing bounding boxes from a mesh and markers. These do not need -// to be axis-aligned. +// to be axis-aligned. Note: This function only returns an oriented bounding box for points +// which exactly form a rectangle or rectangular prism. For other shapes, the result is less +// predictable. For shapes where the convex hull is a circle or sphere, the convex hull +// is the circumscribed sphere of the resulting box, or equivalently the box is inscribed by +// the convex hull. It is thus clearly not a bounding box but can be useful to compute +// nonetheless. BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr); BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr); From d6a869fbb50dec0fc965be9258011300c92eb311 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 13:29:45 -0700 Subject: [PATCH 04/11] Restore GetBoundingBall calculation which uses Welzl's algorithm --- palace/fem/lumpedelement.cpp | 20 +--- palace/utils/geodata.cpp | 205 ++++++++++++++++++++++++++--------- palace/utils/geodata.hpp | 89 ++++++++++++--- 3 files changed, 232 insertions(+), 82 deletions(-) diff --git a/palace/fem/lumpedelement.cpp b/palace/fem/lumpedelement.cpp index 5b5111d30..30770ecca 100644 --- a/palace/fem/lumpedelement.cpp +++ b/palace/fem/lumpedelement.cpp @@ -115,26 +115,18 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, const mfem::ParMesh &mesh = *fespace.GetParMesh(); int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); - auto bounding_box = mesh::GetBoundingBox(mesh, attr_marker, true); - MFEM_VERIFY(bounding_box.planar, + auto bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true); + MFEM_VERIFY(bounding_ball.planar, "Boundary elements must be coplanar to define a coaxial lumped element!"); // Direction of the excitation as +/-r̂. direction = (input_dir[0] > 0); - origin.SetSize(bounding_box.center.size()); - std::copy(bounding_box.center.begin(), bounding_box.center.end(), origin.begin()); + origin.SetSize(bounding_ball.center.size()); + std::copy(bounding_ball.center.begin(), bounding_ball.center.end(), origin.begin()); - // Get inner radius of annulus assuming full 2π circumference. Coarse tolerance to allow - // for approximately-meshed circles. The returned bounding box is not a bounding box, it - // is inscribed in the outer circle. - constexpr double rel_tol = 1.0e-2; + // Get inner radius of annulus assuming full 2π circumference. + r_outer = 0.5 * bounding_ball.Lengths()[0]; double A = GetArea(fespace, attr_marker); - auto lengths = bounding_box.Lengths(); - MFEM_VERIFY((1.0 - rel_tol) * lengths[1] < lengths[0] && - lengths[0] < (1.0 + rel_tol) * lengths[1], - "Bounding box for coaxial lumped port elements is not approximately square (" - << "sides = " << lengths[0] << ", " << lengths[1] << ")!"); - r_outer = 0.25 * std::sqrt(2.0) * (lengths[0] + lengths[1]); MFEM_VERIFY(std::pow(r_outer, 2) - A / M_PI > 0.0, "Coaxial element boundary is not defined correctly (radius " << r_outer << ", area " << A << ")!"); diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index fcfa2192a..996cb49d6 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -6,8 +6,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -679,23 +681,6 @@ void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array Mpi::GlobalMax(dim, max.HostReadWrite(), mesh.GetComm()); } -void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr, - mfem::Vector &min, mfem::Vector &max) -{ - mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); - marker = 0; - marker[attr - 1] = 1; - GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max); -} - -void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, mfem::Vector &min, - mfem::Vector &max) -{ - mfem::Array marker(mesh.attributes.Max()); - marker = 1; - GetAxisAlignedBoundingBox(mesh, marker, false, min, max); -} - double BoundingBox::Area() const { return 4.0 * @@ -735,7 +720,7 @@ namespace bool EigenLE(const Eigen::Vector3d &x, const Eigen::Vector3d &y) { return std::lexicographical_compare(x.begin(), x.end(), y.begin(), y.end()); -}; +} // Helper for collecting a point cloud from a mesh, used in calculating bounding boxes and // bounding balls. Returns the dominant rank, for which the vertices argument will be @@ -895,7 +880,7 @@ auto PerpendicularDistance(const std::initializer_list &normals v0 -= n.dot(v0) * n; } return v0.norm(); -}; +} // Calculates a bounding box from a point cloud, result is broadcast across all processes. BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, @@ -977,7 +962,7 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, : *std::min_element(vertices_out_of_plane.begin(), vertices_out_of_plane.end(), DistFromP_000); const bool t_0_gt_t_1 = - (t_0 - origin).norm() > (t_1 - origin).norm(); // If planar t_1 == t_0 + (t_0 - origin).norm() > (t_1 - origin).norm(); // If planar, t_1 == t_0 const auto &v_001 = t_0_gt_t_1 ? t_1 : t_0; const auto &v_011 = box.planar ? v_111 : (t_0_gt_t_1 ? t_0 : t_1); @@ -1002,6 +987,149 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, return box; } +// Use 4 points to define a sphere in 3D. If the points are coplanar, 3 of them are used to +// define a circle which is interpreted as the equator of the sphere. We assume the points +// are unique and not collinear. +void SphereFromPoints(const std::vector &indices, + const std::vector &vertices, Eigen::Vector3d &origin, + double &radius, bool &coplanar) +{ + // Given 0, 1, or 2 points, just return a radius of 0. + MFEM_VERIFY( + indices.size() < 4, + "Determining a sphere in 3D requires 4 points (and a circle requires 3 points)!"); + if (indices.size() < 3) + { + origin = Eigen::Vector3d::Zero(); + radius = 0.0; + coplanar = true; + return; + } + + // Check for coplanarity. + constexpr double rel_tol = 1.0e-6; + const Eigen::Vector3d AB = vertices[indices[1]] - vertices[indices[0]]; + const Eigen::Vector3d AC = vertices[indices[2]] - vertices[indices[0]]; + const Eigen::Vector3d ABAC = AB.cross(AC); + Eigen::Vector3d AD; + coplanar = (indices.size() < 4); + if (!coplanar) + { + AD = vertices[indices[3]] - vertices[indices[0]]; + coplanar = (std::abs(AD.dot(ABAC)) < rel_tol * AD.norm() * ABAC.norm()); + } + + // Construct a circle passing through 3 points. + // See: https://en.wikipedia.org/wiki/Circumcircle#Higher_dimensions. + if (coplanar) + { + origin = (0.5 / ABAC.squaredNorm()) * + ((AB.squaredNorm() * AC) - (AC.squaredNorm() * AB)).cross(ABAC); + radius = origin.norm(); + origin += vertices[indices[0]]; +#if defined(MFEM_DEBUG) + const auto r1 = (vertices[indices[1]] - origin).norm(); + const auto r2 = (vertices[indices[2]] - origin).norm(); + MFEM_VERIFY((1.0 - rel_tol) * radius < r1 && r1 < (1.0 + rel_tol) * radius && + (1.0 - rel_tol) * radius < r2 && r2 < (1.0 + rel_tol) * radius, + "Invalid circle calculated from 3 points!"); +#endif + return; + } + + // Construct a sphere passing through 4 points. + // See: https://steve.hollasch.net/cgindex/geometry/sphere4pts.html. + Eigen::Matrix3d C; + Eigen::Vector3d d; + const auto s = vertices[indices[0]].squaredNorm(); + C.row(0) = AB.transpose(); + C.row(1) = AC.transpose(); + C.row(2) = AD.transpose(); + d(0) = 0.5 * (vertices[indices[1]].squaredNorm() - s); + d(0) = 0.5 * (vertices[indices[2]].squaredNorm() - s); + d(0) = 0.5 * (vertices[indices[3]].squaredNorm() - s); + origin = C.inverse() * d; // 3x3 matrix inverse might be faster than general LU + radius = (vertices[indices[0]] - origin).norm(); +#if defined(MFEM_DEBUG) + const auto r1 = (vertices[indices[1]] - origin).norm(); + const auto r2 = (vertices[indices[2]] - origin).norm(); + const auto r3 = (vertices[indices[3]] - origin).norm(); + MFEM_VERIFY((1.0 - rel_tol) * radius < r1 && r1 < (1.0 + rel_tol) * radius && + (1.0 - rel_tol) * radius < r2 && r2 < (1.0 + rel_tol) * radius && + (1.0 - rel_tol) * radius < r3 && r3 < (1.0 + rel_tol) * radius, + "Invalid sphere calculated from 3 points!"); +#endif +} + +void Welzl(std::vector P, std::vector R, + const std::vector &vertices, Eigen::Vector3d &origin, + double &radius, bool &planar) +{ + // Base case. + if (R.size() == 4 || P.empty()) + { + return SphereFromPoints(R, vertices, origin, radius, planar); + } + + // Choose a p ∈ P randomly, and recurse for (P \ {p}, R). The set P has already been + // randomized on input. + const std::size_t p = P.back(); + P.pop_back(); + Welzl(P, R, vertices, origin, radius, planar); + + // If p is outside the sphere, recurse for (P \ {p}, R U {p}). + constexpr double rel_tol = 1.0e-6; + if ((vertices[p] - origin).norm() >= (1.0 + rel_tol) * radius) + { + R.push_back(p); + Welzl(P, R, vertices, origin, radius, planar); + } +} + +// Calculates a bounding ball from a point cloud using Welzl's algorithm, result is +// broadcast across all processes. We don't operate on the convex hull, since the number of +// points should be small enough that operating on the full set should be OK. If only three +// points are provided, the bounding circle is computed (likewise for if the points are +// coplanar). +BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, + const std::vector &vertices, + int dominant_rank) +{ + BoundingBox ball; + if (dominant_rank == Mpi::Rank(comm)) + { + // Randomly permute the point set. + MFEM_VERIFY(vertices.size() >= 3, + "A bounding ball requires a minimum of three vertices for this algorithm!"); + std::vector indices(vertices.size()); + for (std::size_t i = 0; i < vertices.size(); i++) + { + indices[i] = i; + } + { + std::random_device rd; + std::mt19937 g(rd()); + std::shuffle(indices.begin(), indices.end(), g); + } + + // Compute the bounding ball. + Eigen::Vector3d origin; + double radius; + Welzl(indices, {}, vertices, origin, radius, ball.planar); + Vector3dMap(ball.center.data()) = origin; + Vector3dMap(ball.normals[0].data()) = Eigen::Vector3d(radius, 0.0, 0.0); + Vector3dMap(ball.normals[1].data()) = Eigen::Vector3d(0.0, radius, 0.0); + Vector3dMap(ball.normals[2].data()) = Eigen::Vector3d(0.0, 0.0, radius); + } + + // Broadcast result to all processors. + Mpi::Broadcast(3, ball.center.data(), dominant_rank, comm); + Mpi::Broadcast(3 * 3, ball.normals.data()->data(), dominant_rank, comm); + Mpi::Broadcast(1, &ball.planar, dominant_rank, comm); + + return ball; +} + double LengthFromPointCloud(MPI_Comm comm, const std::vector &vertices, int dominant_rank, const std::array &dir) { @@ -1009,12 +1137,10 @@ double LengthFromPointCloud(MPI_Comm comm, const std::vector &v if (dominant_rank == Mpi::Rank(comm)) { CVector3dMap direction(dir.data()); - auto Dot = [&](const auto &x, const auto &y) { return direction.dot(x) < direction.dot(y); }; auto p_min = std::min_element(vertices.begin(), vertices.end(), Dot); auto p_max = std::max_element(vertices.begin(), vertices.end(), Dot); - length = (*p_max - *p_min).dot(direction.normalized()); } Mpi::Broadcast(1, &length, dominant_rank, comm); @@ -1031,12 +1157,12 @@ BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &ma return BoundingBoxFromPointCloud(mesh.GetComm(), vertices, dominant_rank); } -BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr) +BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &marker, + bool bdr) { - mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); - marker = 0; - marker[attr - 1] = 1; - return GetBoundingBox(mesh, marker, bdr); + std::vector vertices; + int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices); + return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank); } double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &marker, @@ -1047,15 +1173,6 @@ double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &mar return LengthFromPointCloud(mesh.GetComm(), vertices, dominant_rank, dir); } -double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, - const std::array &dir) -{ - mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); - marker = 0; - marker[attr - 1] = 1; - return GetProjectedLength(mesh, marker, bdr, dir); -} - mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array &marker, bool average) { @@ -1174,22 +1291,6 @@ mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array return normal; } -mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr, bool average) -{ - const bool bdr = (mesh.Dimension() == mesh.SpaceDimension()); - mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); - marker = 0; - marker[attr - 1] = 1; - return GetSurfaceNormal(mesh, marker, average); -} - -mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average) -{ - const bool bdr = (mesh.Dimension() == mesh.SpaceDimension()); - const auto &attributes = bdr ? mesh.bdr_attributes : mesh.attributes; - return GetSurfaceNormal(mesh, AttrToMarker(attributes.Max(), attributes), average); -} - double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) { BlockTimer bt0(Timer::REBALANCE); diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 3b3bba66a..df29ea990 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -59,8 +59,10 @@ ElementTypeInfo CheckElements(const mfem::Mesh &mesh); template void AttrToMarker(int max_attr, const T &attr_list, mfem::Array &marker, bool skip_invalid = false); + template -mfem::Array AttrToMarker(int max_attr, const T &attr_list, bool skip_invalid = false) +inline mfem::Array AttrToMarker(int max_attr, const T &attr_list, + bool skip_invalid = false) { mfem::Array marker; AttrToMarker(max_attr, attr_list, marker, skip_invalid); @@ -71,10 +73,23 @@ mfem::Array AttrToMarker(int max_attr, const T &attr_list, bool skip_invali // given attribute. void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, mfem::Vector &min, mfem::Vector &max); -void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr, - mfem::Vector &min, mfem::Vector &max); -void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, mfem::Vector &min, - mfem::Vector &max); + +inline void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr, + mfem::Vector &min, mfem::Vector &max) +{ + mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); + marker = 0; + marker[attr - 1] = 1; + GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max); +} + +inline void GetAxisAlignedBoundingBox(const mfem::ParMesh &mesh, mfem::Vector &min, + mfem::Vector &max) +{ + mfem::Array marker(mesh.attributes.Max()); + marker = 1; + GetAxisAlignedBoundingBox(mesh, marker, false, min, max); +} // Struct describing a bounding box in terms of the center and face normals. The normals // specify the direction from the center of the box. @@ -103,28 +118,70 @@ struct BoundingBox }; // Helper functions for computing bounding boxes from a mesh and markers. These do not need -// to be axis-aligned. Note: This function only returns an oriented bounding box for points -// which exactly form a rectangle or rectangular prism. For other shapes, the result is less -// predictable. For shapes where the convex hull is a circle or sphere, the convex hull -// is the circumscribed sphere of the resulting box, or equivalently the box is inscribed by -// the convex hull. It is thus clearly not a bounding box but can be useful to compute -// nonetheless. +// to be axis-aligned. Note: This function only returns a minimum oriented bounding box for +// points whose convex hull exactly forms a rectangle or rectangular prism, implementing a +// vastly simplified version of QuickHull for this case. For other shapes, the result is +// less predictable, and may not even form a bounding box of the sampled point cloud. BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr); -BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr); + +inline BoundingBox GetBoundingBox(const mfem::ParMesh &mesh, int attr, bool bdr) +{ + mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); + marker = 0; + marker[attr - 1] = 1; + return GetBoundingBox(mesh, marker, bdr); +} + +// Given a mesh and a marker, compute the bounding circle/sphere of the marked elements. In +// this case the normals of the bounding box object are arbitrary, and the Area and Volume +// members should not be used, but the Lengths function returns the ball diameter. This +// function implements Welzl's algorithm. +BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &marker, + bool bdr); + +inline BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, int attr, bool bdr) +{ + mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); + marker = 0; + marker[attr - 1] = 1; + return GetBoundingBall(mesh, marker, bdr); +} // Helper function for computing the direction aligned length of a marked group. double GetProjectedLength(const mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, const std::array &dir); -double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, - const std::array &dir); + +inline double GetProjectedLength(const mfem::ParMesh &mesh, int attr, bool bdr, + const std::array &dir) +{ + mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); + marker = 0; + marker[attr - 1] = 1; + return GetProjectedLength(mesh, marker, bdr, dir); +} // Helper function to compute the average surface normal for all elements with the given // attribute. mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, const mfem::Array &marker, bool average = true); -mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr, bool average = true); -mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average = true); + +inline mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, int attr, + bool average = true) +{ + const bool bdr = (mesh.Dimension() == mesh.SpaceDimension()); + mfem::Array marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max()); + marker = 0; + marker[attr - 1] = 1; + return GetSurfaceNormal(mesh, marker, average); +} + +inline mfem::Vector GetSurfaceNormal(const mfem::ParMesh &mesh, bool average = true) +{ + const bool bdr = (mesh.Dimension() == mesh.SpaceDimension()); + const auto &attributes = bdr ? mesh.bdr_attributes : mesh.attributes; + return GetSurfaceNormal(mesh, AttrToMarker(attributes.Max(), attributes), average); +} // Helper function responsible for rebalancing the mesh, and optionally writing meshes from // the intermediate stages to disk. Returns the imbalance ratio before rebalancing. From 70de55e54b16202d42548aaa5869e475ebe0658d Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 13:30:24 -0700 Subject: [PATCH 05/11] Improve implementation for minumum bounding ball, test 3D case --- palace/fem/lumpedelement.cpp | 20 +++- palace/utils/geodata.cpp | 202 ++++++++++++++++++++++++----------- 2 files changed, 160 insertions(+), 62 deletions(-) diff --git a/palace/fem/lumpedelement.cpp b/palace/fem/lumpedelement.cpp index 30770ecca..291d31179 100644 --- a/palace/fem/lumpedelement.cpp +++ b/palace/fem/lumpedelement.cpp @@ -116,8 +116,8 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); auto bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true); - MFEM_VERIFY(bounding_ball.planar, - "Boundary elements must be coplanar to define a coaxial lumped element!"); + // MFEM_VERIFY(bounding_ball.planar, + // "Boundary elements must be coplanar to define a coaxial lumped element!"); // Direction of the excitation as +/-r̂. direction = (input_dir[0] > 0); @@ -131,6 +131,22 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, "Coaxial element boundary is not defined correctly (radius " << r_outer << ", area " << A << ")!"); r_inner = std::sqrt(std::pow(r_outer, 2) - A / M_PI); + + + + + // XX TODO WIP... + const auto ¢er = bounding_ball.center; + const auto &lengths = bounding_ball.Lengths(); + std::cout << "\norigin: " << center[0] << ", " << center[1] << ", " << center[2] + << "\n\n"; + std::cout << "\nlengths: " << lengths[0] << ", " << lengths[1] << ", " << lengths[2] + << "\n\n"; + std::cout << "\nr_outer: " << r_outer << "\n\n"; + std::cout << "\nr_inner: " << r_inner << "\n\n"; + + + } std::unique_ptr diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 996cb49d6..dcfe717cc 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -899,13 +899,17 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, MFEM_VERIFY(vertices.size() >= 4, "A bounding box requires a minimum of four vertices for this algorithm!"); auto p_000 = std::min_element(vertices.begin(), vertices.end(), EigenLE); - auto DistFromP_000 = [&p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - *p_000).norm() < (y - *p_000).norm(); }; - auto p_111 = std::max_element(vertices.begin(), vertices.end(), DistFromP_000); - auto DistFromP_111 = [&p_111](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - *p_111).norm() < (y - *p_111).norm(); }; - p_000 = std::max_element(vertices.begin(), vertices.end(), DistFromP_111); - MFEM_VERIFY(std::max_element(vertices.begin(), vertices.end(), DistFromP_000) == p_111, + auto p_111 = + std::max_element(vertices.begin(), vertices.end(), + [p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_000).norm() < (y - *p_000).norm(); }); + p_000 = std::max_element(vertices.begin(), vertices.end(), + [p_111](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_111).norm() < (y - *p_111).norm(); }); + MFEM_ASSERT(std::max_element(vertices.begin(), vertices.end(), + [p_000](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_000).norm() < (y - *p_000).norm(); }) == + p_111, "p_000 and p_111 must be mutually opposing points!"); // Define a diagonal of the ASSUMED cuboid bounding box. @@ -957,10 +961,12 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, }); // Given candidates t_0 and t_1, the closer to origin defines v_001. - const auto &t_1 = box.planar - ? t_0 - : *std::min_element(vertices_out_of_plane.begin(), - vertices_out_of_plane.end(), DistFromP_000); + const auto &t_1 = + box.planar + ? t_0 + : *std::min_element(vertices_out_of_plane.begin(), vertices_out_of_plane.end(), + [&](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - origin).norm() < (y - origin).norm(); }); const bool t_0_gt_t_1 = (t_0 - origin).norm() > (t_1 - origin).norm(); // If planar, t_1 == t_0 const auto &v_001 = t_0_gt_t_1 ? t_1 : t_0; @@ -987,23 +993,41 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, return box; } +// For the public interface, we can use a BoundingBox as a generalization of a BoundingBall. +// Internally, however, it's nice to work with a specific ball data type. +struct BoundingBall +{ + Eigen::Vector3d origin; + double radius; + bool planar; +}; + // Use 4 points to define a sphere in 3D. If the points are coplanar, 3 of them are used to // define a circle which is interpreted as the equator of the sphere. We assume the points // are unique and not collinear. -void SphereFromPoints(const std::vector &indices, - const std::vector &vertices, Eigen::Vector3d &origin, - double &radius, bool &coplanar) +BoundingBall SphereFromPoints(const std::vector &indices, + const std::vector &vertices) { - // Given 0, 1, or 2 points, just return a radius of 0. + // Given 0 or 1 points, just return a radius of 0. MFEM_VERIFY( - indices.size() < 4, + indices.size() <= 4, "Determining a sphere in 3D requires 4 points (and a circle requires 3 points)!"); - if (indices.size() < 3) + BoundingBall ball; + ball.planar = (indices.size() < 4); + if (indices.size() < 2) + { + ball.origin = Eigen::Vector3d::Zero(); + ball.radius = 0.0; + return ball; + } + + // For two points, construct a circle with the segment as its diameter. This could also + // handle the collinear case for more than 2 points. + if (indices.size() == 2) { - origin = Eigen::Vector3d::Zero(); - radius = 0.0; - coplanar = true; - return; + ball.origin = 0.5 * (vertices[indices[0]] + vertices[indices[1]]); + ball.radius = (vertices[indices[0]] - ball.origin).norm(); + return ball; } // Check for coplanarity. @@ -1012,29 +1036,29 @@ void SphereFromPoints(const std::vector &indices, const Eigen::Vector3d AC = vertices[indices[2]] - vertices[indices[0]]; const Eigen::Vector3d ABAC = AB.cross(AC); Eigen::Vector3d AD; - coplanar = (indices.size() < 4); - if (!coplanar) + if (!ball.planar) { AD = vertices[indices[3]] - vertices[indices[0]]; - coplanar = (std::abs(AD.dot(ABAC)) < rel_tol * AD.norm() * ABAC.norm()); + ball.planar = (std::abs(AD.dot(ABAC)) < rel_tol * AD.norm() * ABAC.norm()); } // Construct a circle passing through 3 points. // See: https://en.wikipedia.org/wiki/Circumcircle#Higher_dimensions. - if (coplanar) + if (ball.planar) { - origin = (0.5 / ABAC.squaredNorm()) * - ((AB.squaredNorm() * AC) - (AC.squaredNorm() * AB)).cross(ABAC); - radius = origin.norm(); - origin += vertices[indices[0]]; + ball.origin = (0.5 / ABAC.squaredNorm()) * + ((AB.squaredNorm() * AC) - (AC.squaredNorm() * AB)).cross(ABAC); + ball.radius = ball.origin.norm(); + ball.origin += vertices[indices[0]]; #if defined(MFEM_DEBUG) - const auto r1 = (vertices[indices[1]] - origin).norm(); - const auto r2 = (vertices[indices[2]] - origin).norm(); - MFEM_VERIFY((1.0 - rel_tol) * radius < r1 && r1 < (1.0 + rel_tol) * radius && - (1.0 - rel_tol) * radius < r2 && r2 < (1.0 + rel_tol) * radius, + const auto r1 = (vertices[indices[1]] - ball.origin).norm(); + const auto r2 = (vertices[indices[2]] - ball.origin).norm(); + MFEM_VERIFY((1.0 - rel_tol) * ball.radius < r1 && r1 < (1.0 + rel_tol) * ball.radius && + (1.0 - rel_tol) * ball.radius < r2 && + r2 < (1.0 + rel_tol) * ball.radius, "Invalid circle calculated from 3 points!"); #endif - return; + return ball; } // Construct a sphere passing through 4 points. @@ -1046,44 +1070,48 @@ void SphereFromPoints(const std::vector &indices, C.row(1) = AC.transpose(); C.row(2) = AD.transpose(); d(0) = 0.5 * (vertices[indices[1]].squaredNorm() - s); - d(0) = 0.5 * (vertices[indices[2]].squaredNorm() - s); - d(0) = 0.5 * (vertices[indices[3]].squaredNorm() - s); - origin = C.inverse() * d; // 3x3 matrix inverse might be faster than general LU - radius = (vertices[indices[0]] - origin).norm(); + d(1) = 0.5 * (vertices[indices[2]].squaredNorm() - s); + d(2) = 0.5 * (vertices[indices[3]].squaredNorm() - s); + ball.origin = C.inverse() * d; // 3x3 matrix inverse might be faster than general LU + // if Eigen uses the explicit closed-form solution + ball.radius = (vertices[indices[0]] - ball.origin).norm(); #if defined(MFEM_DEBUG) - const auto r1 = (vertices[indices[1]] - origin).norm(); - const auto r2 = (vertices[indices[2]] - origin).norm(); - const auto r3 = (vertices[indices[3]] - origin).norm(); - MFEM_VERIFY((1.0 - rel_tol) * radius < r1 && r1 < (1.0 + rel_tol) * radius && - (1.0 - rel_tol) * radius < r2 && r2 < (1.0 + rel_tol) * radius && - (1.0 - rel_tol) * radius < r3 && r3 < (1.0 + rel_tol) * radius, + const auto r1 = (vertices[indices[1]] - ball.origin).norm(); + const auto r2 = (vertices[indices[2]] - ball.origin).norm(); + const auto r3 = (vertices[indices[3]] - ball.origin).norm(); + MFEM_VERIFY((1.0 - rel_tol) * ball.radius < r1 && r1 < (1.0 + rel_tol) * ball.radius && + (1.0 - rel_tol) * ball.radius < r2 && + r2 < (1.0 + rel_tol) * ball.radius && + (1.0 - rel_tol) * ball.radius < r3 && r3 < (1.0 + rel_tol) * ball.radius, "Invalid sphere calculated from 3 points!"); #endif + return ball; } -void Welzl(std::vector P, std::vector R, - const std::vector &vertices, Eigen::Vector3d &origin, - double &radius, bool &planar) +BoundingBall Welzl(std::vector P, std::vector R, + const std::vector &vertices) { // Base case. if (R.size() == 4 || P.empty()) { - return SphereFromPoints(R, vertices, origin, radius, planar); + return SphereFromPoints(R, vertices); } // Choose a p ∈ P randomly, and recurse for (P \ {p}, R). The set P has already been // randomized on input. const std::size_t p = P.back(); P.pop_back(); - Welzl(P, R, vertices, origin, radius, planar); + BoundingBall D = Welzl(P, R, vertices); // If p is outside the sphere, recurse for (P \ {p}, R U {p}). constexpr double rel_tol = 1.0e-6; - if ((vertices[p] - origin).norm() >= (1.0 + rel_tol) * radius) + if ((vertices[p] - D.origin).norm() >= (1.0 + rel_tol) * D.radius) { R.push_back(p); - Welzl(P, R, vertices, origin, radius, planar); + D = Welzl(P, R, vertices); } + + return D; } // Calculates a bounding ball from a point cloud using Welzl's algorithm, result is @@ -1101,7 +1129,7 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, // Randomly permute the point set. MFEM_VERIFY(vertices.size() >= 3, "A bounding ball requires a minimum of three vertices for this algorithm!"); - std::vector indices(vertices.size()); + std::vector indices(vertices.size() + 4); for (std::size_t i = 0; i < vertices.size(); i++) { indices[i] = i; @@ -1109,17 +1137,58 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, { std::random_device rd; std::mt19937 g(rd()); - std::shuffle(indices.begin(), indices.end(), g); + std::shuffle(indices.begin(), indices.end() - 4, g); + } + + // Acceleration from https://informatica.vu.lt/journal/INFORMATICA/article/1251. Allow + // for duplicate points and just add the 4 points to the end of the indicies list to be + // considered first. The two points are not necessarily the maximizer of the distance + // between all pairs, but they should be a good estimate. + { + auto p_1 = std::min_element(vertices.begin(), vertices.end(), EigenLE); + auto p_2 = std::max_element(vertices.begin(), vertices.end(), + [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_1).norm() < (y - *p_1).norm(); }); + p_1 = std::max_element(vertices.begin(), vertices.end(), + [p_2](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_2).norm() < (y - *p_2).norm(); }); + + + // //XX TODO + // auto p_test = std::max_element(vertices.begin(), vertices.end(), + // [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + // { return (x - *p_1).norm() < (y - *p_1).norm(); }); + // MFEM_ASSERT(p_test == p_2, "p_1 and p_2 must be mutually opposing points!"); + + + // MFEM_ASSERT(std::max_element(vertices.begin(), vertices.end(), + // [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + // { return (x - *p_1).norm() < (y - *p_1).norm(); }) == + // p_2, + // "p_1 and p_2 must be mutually opposing points!"); + + + auto p_12 = 0.5 * (*p_1 + *p_2); + auto p_3 = + std::max_element(vertices.begin(), vertices.end(), + [&p_12](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - p_12).norm() < (y - p_12).norm(); }); + auto p_4 = std::max_element(vertices.begin(), vertices.end(), + [p_3](const Eigen::Vector3d &x, const Eigen::Vector3d &y) + { return (x - *p_3).norm() < (y - *p_3).norm(); }); + indices[indices.size() - 1] = p_1 - vertices.begin(); + indices[indices.size() - 2] = p_2 - vertices.begin(); + indices[indices.size() - 3] = p_3 - vertices.begin(); + indices[indices.size() - 4] = p_4 - vertices.begin(); } // Compute the bounding ball. - Eigen::Vector3d origin; - double radius; - Welzl(indices, {}, vertices, origin, radius, ball.planar); - Vector3dMap(ball.center.data()) = origin; - Vector3dMap(ball.normals[0].data()) = Eigen::Vector3d(radius, 0.0, 0.0); - Vector3dMap(ball.normals[1].data()) = Eigen::Vector3d(0.0, radius, 0.0); - Vector3dMap(ball.normals[2].data()) = Eigen::Vector3d(0.0, 0.0, radius); + BoundingBall min_ball = Welzl(indices, {}, vertices); + Vector3dMap(ball.center.data()) = min_ball.origin; + Vector3dMap(ball.normals[0].data()) = Eigen::Vector3d(min_ball.radius, 0.0, 0.0); + Vector3dMap(ball.normals[1].data()) = Eigen::Vector3d(0.0, min_ball.radius, 0.0); + Vector3dMap(ball.normals[2].data()) = Eigen::Vector3d(0.0, 0.0, min_ball.radius); + ball.planar = min_ball.planar; } // Broadcast result to all processors. @@ -1162,6 +1231,19 @@ BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &m { std::vector vertices; int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices); + + + + //XX TODO WIP FOR TEST... + bool neg = false; + for (auto &v : vertices) + { + v[0] += 0.1 * (neg ? -1.0 : 1.0); + neg = !neg; + } + + + return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank); } From f98872d26fe930a5fe936027780fb6dad8e180c9 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 13:31:11 -0700 Subject: [PATCH 06/11] Remove debugging code --- palace/fem/lumpedelement.cpp | 20 ++------------------ palace/utils/geodata.cpp | 29 ----------------------------- 2 files changed, 2 insertions(+), 47 deletions(-) diff --git a/palace/fem/lumpedelement.cpp b/palace/fem/lumpedelement.cpp index 291d31179..30770ecca 100644 --- a/palace/fem/lumpedelement.cpp +++ b/palace/fem/lumpedelement.cpp @@ -116,8 +116,8 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0; mfem::Array attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list); auto bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true); - // MFEM_VERIFY(bounding_ball.planar, - // "Boundary elements must be coplanar to define a coaxial lumped element!"); + MFEM_VERIFY(bounding_ball.planar, + "Boundary elements must be coplanar to define a coaxial lumped element!"); // Direction of the excitation as +/-r̂. direction = (input_dir[0] > 0); @@ -131,22 +131,6 @@ CoaxialElementData::CoaxialElementData(const std::array &input_dir, "Coaxial element boundary is not defined correctly (radius " << r_outer << ", area " << A << ")!"); r_inner = std::sqrt(std::pow(r_outer, 2) - A / M_PI); - - - - - // XX TODO WIP... - const auto ¢er = bounding_ball.center; - const auto &lengths = bounding_ball.Lengths(); - std::cout << "\norigin: " << center[0] << ", " << center[1] << ", " << center[2] - << "\n\n"; - std::cout << "\nlengths: " << lengths[0] << ", " << lengths[1] << ", " << lengths[2] - << "\n\n"; - std::cout << "\nr_outer: " << r_outer << "\n\n"; - std::cout << "\nr_inner: " << r_inner << "\n\n"; - - - } std::unique_ptr diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index dcfe717cc..4462640c8 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1152,22 +1152,6 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, p_1 = std::max_element(vertices.begin(), vertices.end(), [p_2](const Eigen::Vector3d &x, const Eigen::Vector3d &y) { return (x - *p_2).norm() < (y - *p_2).norm(); }); - - - // //XX TODO - // auto p_test = std::max_element(vertices.begin(), vertices.end(), - // [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - // { return (x - *p_1).norm() < (y - *p_1).norm(); }); - // MFEM_ASSERT(p_test == p_2, "p_1 and p_2 must be mutually opposing points!"); - - - // MFEM_ASSERT(std::max_element(vertices.begin(), vertices.end(), - // [p_1](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - // { return (x - *p_1).norm() < (y - *p_1).norm(); }) == - // p_2, - // "p_1 and p_2 must be mutually opposing points!"); - - auto p_12 = 0.5 * (*p_1 + *p_2); auto p_3 = std::max_element(vertices.begin(), vertices.end(), @@ -1231,19 +1215,6 @@ BoundingBox GetBoundingBall(const mfem::ParMesh &mesh, const mfem::Array &m { std::vector vertices; int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices); - - - - //XX TODO WIP FOR TEST... - bool neg = false; - for (auto &v : vertices) - { - v[0] += 0.1 * (neg ? -1.0 : 1.0); - neg = !neg; - } - - - return BoundingBallFromPointCloud(mesh.GetComm(), vertices, dominant_rank); } From 78e2de96765918d291df979b2773231f9da8e456 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 14:04:50 -0700 Subject: [PATCH 07/11] Remove unneeded include --- palace/utils/geodata.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 4462640c8..8814f89a4 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -6,7 +6,6 @@ #include #include #include -#include #include #include #include From 32ebd0c207e95e6a463007517689f75af0240dd1 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 15:50:08 -0700 Subject: [PATCH 08/11] Address PR feedback and improve initialization for Welzl --- palace/utils/geodata.cpp | 47 +++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 8814f89a4..9de9d1e23 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1125,19 +1125,10 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, BoundingBox ball; if (dominant_rank == Mpi::Rank(comm)) { - // Randomly permute the point set. MFEM_VERIFY(vertices.size() >= 3, "A bounding ball requires a minimum of three vertices for this algorithm!"); - std::vector indices(vertices.size() + 4); - for (std::size_t i = 0; i < vertices.size(); i++) - { - indices[i] = i; - } - { - std::random_device rd; - std::mt19937 g(rd()); - std::shuffle(indices.begin(), indices.end() - 4, g); - } + std::vector indices(vertices.size()); + std::iota(indices.begin(), indices.end(), 0); // Acceleration from https://informatica.vu.lt/journal/INFORMATICA/article/1251. Allow // for duplicate points and just add the 4 points to the end of the indicies list to be @@ -1151,18 +1142,34 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, p_1 = std::max_element(vertices.begin(), vertices.end(), [p_2](const Eigen::Vector3d &x, const Eigen::Vector3d &y) { return (x - *p_2).norm() < (y - *p_2).norm(); }); - auto p_12 = 0.5 * (*p_1 + *p_2); - auto p_3 = - std::max_element(vertices.begin(), vertices.end(), - [&p_12](const Eigen::Vector3d &x, const Eigen::Vector3d &y) - { return (x - p_12).norm() < (y - p_12).norm(); }); + + // Find the p_3 as the vertex furthest from the initial axis. + const Eigen::Vector3d n_1 = (*p_2 - *p_1).normalized(); + auto p_3 = std::max_element(vertices.begin(), vertices.end(), + [&](const auto &x, const auto &y) { + return PerpendicularDistance({n_1}, *p_1, x) < + PerpendicularDistance({n_1}, *p_1, y); + }); auto p_4 = std::max_element(vertices.begin(), vertices.end(), [p_3](const Eigen::Vector3d &x, const Eigen::Vector3d &y) { return (x - *p_3).norm() < (y - *p_3).norm(); }); - indices[indices.size() - 1] = p_1 - vertices.begin(); - indices[indices.size() - 2] = p_2 - vertices.begin(); - indices[indices.size() - 3] = p_3 - vertices.begin(); - indices[indices.size() - 4] = p_4 - vertices.begin(); + MFEM_VERIFY(p_3 != p_1 && p_3 != p_2 && p_4 != p_1 && p_4 != p_2, + "Vertices are degenerate!"); + + // Start search with these points, which should be roughly extremal. With the search + // for p_3 done in an orthogonal direction, p_1, p_2, p_3, and p_4 should all be + // unique. + std::swap(indices[indices.size() - 1], indices[p_1 - vertices.begin()]); + std::swap(indices[indices.size() - 2], indices[p_2 - vertices.begin()]); + std::swap(indices[indices.size() - 3], indices[p_3 - vertices.begin()]); + std::swap(indices[indices.size() - 4], indices[p_4 - vertices.begin()]); + } + + // Randomly permute the point set. + { + std::random_device rd; + std::mt19937 g(rd()); + std::shuffle(indices.begin(), indices.end() - 4, g); } // Compute the bounding ball. From c82393ec119860599586ee98b2cab5c328226ee6 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 15:59:54 -0700 Subject: [PATCH 09/11] Silence uninitialized variable warning --- palace/utils/geodata.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 9de9d1e23..67e25ab88 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1034,7 +1034,7 @@ BoundingBall SphereFromPoints(const std::vector &indices, const Eigen::Vector3d AB = vertices[indices[1]] - vertices[indices[0]]; const Eigen::Vector3d AC = vertices[indices[2]] - vertices[indices[0]]; const Eigen::Vector3d ABAC = AB.cross(AC); - Eigen::Vector3d AD; + Eigen::Vector3d AD = Eigen::Vector3d::Zero(); if (!ball.planar) { AD = vertices[indices[3]] - vertices[indices[0]]; From cfc317e65c491a1c3f14f375df6809b30380b370 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 16:03:23 -0700 Subject: [PATCH 10/11] Revise typo in comment --- palace/utils/geodata.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 67e25ab88..a8ff0da65 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -1143,7 +1143,7 @@ BoundingBox BoundingBallFromPointCloud(MPI_Comm comm, [p_2](const Eigen::Vector3d &x, const Eigen::Vector3d &y) { return (x - *p_2).norm() < (y - *p_2).norm(); }); - // Find the p_3 as the vertex furthest from the initial axis. + // Find the next point as the vertex furthest from the initial axis. const Eigen::Vector3d n_1 = (*p_2 - *p_1).normalized(); auto p_3 = std::max_element(vertices.begin(), vertices.end(), [&](const auto &x, const auto &y) { From d4c52bdceac4e8b90c293fba91b419fa6ec54f4a Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 8 May 2024 16:04:50 -0700 Subject: [PATCH 11/11] Update changelog --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15118fd7d..32b3ca85e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -43,6 +43,8 @@ The format of this changelog is based on - Fixed a bug when computing the energy associated with lumped elements with more than one nonzero R, L, or C. This also affects the inductive EPR for lumped inductors with and associated parallel capacitance. + - Fixed a bug for coaxial lumped ports which led to incorrect extraction of the geometric + parameters, especially when coarsely-meshed or non-axis-aligned. ## [0.12.0] - 2023-12-21