From b2c5266e68fbf3acf9836832f72b4e134bffe7db Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Wed, 23 Aug 2023 17:43:11 -0400 Subject: [PATCH 1/5] Fix bug in bounding box failing to account for origin shift, and switch to more robust detection of out of plane --- palace/utils/geodata.cpp | 73 ++++++++++++++++++---------------------- 1 file changed, 32 insertions(+), 41 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index ab7cbc72d..b6619c551 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -764,39 +764,31 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // pick the greatest deviation. const Eigen::Vector3d n_3 = n_1.cross(n_2).normalized(); - auto OutOfPlaneComp = [&n_1, &n_2, &n_3, &origin](const auto &v_1, const auto &v_2) + auto OutOfPlaneDistance = [&n_1, &n_2, &origin](const Eigen::Vector3d &v) { - // Precedence of directions is in reverse order of discovery of the directions. The - // most important deciding feature is the distance in the out of plane direction, - // then in the first off-diagonal, then finally in the diagonal direction. - const std::array dist_1{(v_1 - origin).dot(n_3), (v_1 - origin).dot(n_2), - (v_1 - origin).dot(n_1)}; - const std::array dist_2{(v_2 - origin).dot(n_3), (v_2 - origin).dot(n_2), - (v_2 - origin).dot(n_1)}; - return dist_1 < dist_2; + return ((v - origin) - (v - origin).dot(n_1) * n_1 - (v - origin).dot(n_2) * n_2) + .norm(); }; - // There is a degeneration if the final point is within the plane defined by (v_000, - // v_001, v_111). - const auto &t_1 = *std::max_element(vertices.begin(), vertices.end(), OutOfPlaneComp); - constexpr double planar_tolerance = 1.0e-9; - box.planar = std::abs(n_3.dot(t_0)) < planar_tolerance && - std::abs(n_3.dot(t_1)) < planar_tolerance; - if (!box.planar) - { - MFEM_VERIFY(&t_0 != &t_1, "Degenerate coordinates!"); - MFEM_VERIFY(&t_0 != &v_000, "Degenerate coordinates!"); - MFEM_VERIFY(&t_0 != &v_111, "Degenerate coordinates!"); - MFEM_VERIFY(&t_1 != &v_000, "Degenerate coordinates!"); - MFEM_VERIFY(&t_1 != &v_111, "Degenerate coordinates!"); - } - - // If t_1 points to v_000, t_0 or v_111, then the data is coplanar. Establish if t_0 is - // a diagonal or not (using Pythagoras). Only pick t_1 for v_001 if the points are non- - // planar, and t_0 is longer. - bool t_0_gt_t_1 = t_0.norm() >= t_1.norm(); - const auto &v_001 = !box.planar && t_0_gt_t_1 ? t_1 : t_0; - const auto &v_011 = !box.planar && t_0_gt_t_1 ? t_0 : t_1; + // Filter the list of vertices to remove those that are within the plane. Planar + // detection is done using a tolerance given floating point precision. + std::vector vertices_out_of_plane; + const double planar_tolerance = 1.0e-3 * (v_111 - v_000).norm(); + std::copy_if(vertices.begin(), vertices.end(), + std::back_inserter(vertices_out_of_plane), + [OutOfPlaneDistance, planar_tolerance](const auto &v) + { return OutOfPlaneDistance(v) > planar_tolerance; }); + box.planar = vertices_out_of_plane.empty(); + + // Given candidates t_0 and t_1, the closer to origin defines v_001. + const auto &t_1 = box.planar + ? t_0 + : *std::max_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 + 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); // Compute the center as halfway along the main diagonal. Vector3dMap(box.center.data()) = 0.5 * (v_000 + v_111); @@ -804,10 +796,8 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // The length in each direction is then given by traversing the edges of the cuboid in // turn. Vector3dMap(box.normals[0].data()) = 0.5 * (v_001 - v_000); - Vector3dMap(box.normals[1].data()) = - box.planar ? (0.5 * (v_111 - v_001)).eval() : (0.5 * (v_011 - v_001)).eval(); - Vector3dMap(box.normals[2].data()) = - box.planar ? Eigen::Vector3d(0, 0, 0) : (0.5 * (v_111 - v_011)).eval(); + Vector3dMap(box.normals[1].data()) = 0.5 * (v_011 - v_001); + Vector3dMap(box.normals[2].data()) = 0.5 * (v_111 - v_011); // Make sure the longest dimension comes first. std::sort(box.normals.begin(), box.normals.end(), @@ -875,10 +865,10 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, vertices.begin(), vertices.end(), [&delta, PerpendicularDistance](const auto &x, const auto &y) { return PerpendicularDistance({delta}, x) < PerpendicularDistance({delta}, y); }); - constexpr double radius_tol = 1.0e-6; + constexpr double radius_tol = 1.0e-3; MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <= radius_tol * ball.radius, - "A point perpendicular must be contained in the ball: " + "Furthest point perpendicular must be on the exterior of the ball: " << PerpendicularDistance({delta}, perp) << " vs. " << ball.radius << "!"); @@ -887,7 +877,7 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, Vector3dMap(ball.planar_normal.data()) = delta.cross(n_radial).normalized(); // Compute the point furthest out of the plane discovered. If below tolerance, this - // means the circle is 2D. + // 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) @@ -896,13 +886,14 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, PerpendicularDistance({delta, n_radial}, y); }); - constexpr double planar_tolerance = 1.0e-9; - ball.planar = PerpendicularDistance({delta, n_radial}, out_of_plane) < planar_tolerance; + ball.planar = + PerpendicularDistance({delta, n_radial}, out_of_plane) < radius_tol * ball.radius; if (!ball.planar) { // The points are not functionally coplanar, zero out the normal. - MFEM_VERIFY(PerpendicularDistance({delta, n_radial}, out_of_plane) <= ball.radius, - "A point perpendicular must be contained in the ball!"); + MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <= + radius_tol * ball.radius, + "Furthest point perpendicular must be on the exterior of the sphere!"); Vector3dMap(ball.planar_normal.data()) *= 0; } } From c81777eb631706cfdf0a0cb8e50b467aa2c5ee2e Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Fri, 25 Aug 2023 11:11:54 -0400 Subject: [PATCH 2/5] Add more diagonistic information --- palace/fem/lumpedelement.hpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index 807f5b389..a5d6ebe9f 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -69,7 +69,16 @@ class UniformElementData : public LumpedElementData bounding_box(mesh::GetBoundingBox(*fespace.GetParMesh(), marker, true)), direction(3) { MFEM_VERIFY(bounding_box.planar, - "Boundary elements must be coplanar to define a lumped element!"); + [this]() + { + std::stringstream ss; + ss << "Boundary elements must be coplanar to define a lumped element:\n"; + auto l = bounding_box.Lengths(); + ss << "Dimensions (" << l[0] << " x " << l[1] << " x " << l[2] << ")\n"; + ss << "Area " << bounding_box.Area() << "\n"; + ss << "Volume " << bounding_box.Volume() << "\n"; + return ss.str(); + }()); // Check that the bounding box discovered matches the area. This validates that the // boundary elements form a right angled quadrilateral port. From 27377e6443c3114044a931d070f2eb534854d45a Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Fri, 25 Aug 2023 17:07:31 -0400 Subject: [PATCH 3/5] Fix for better treatment of topologically two dimensional objects embedded in three dimensions - Additionally introduce GetDirectionalExtent which computes the maximum separation between points within a geometric object, aligned with a given direction. --- palace/fem/lumpedelement.hpp | 53 +++++++++++++++-------------- palace/utils/geodata.cpp | 66 ++++++++++++++++++++++++++++++------ palace/utils/geodata.hpp | 6 ++++ 3 files changed, 89 insertions(+), 36 deletions(-) diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index a5d6ebe9f..779d6b5d8 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -68,34 +68,22 @@ class UniformElementData : public LumpedElementData : LumpedElementData(fespace.GetParMesh()->SpaceDimension(), marker), bounding_box(mesh::GetBoundingBox(*fespace.GetParMesh(), marker, true)), direction(3) { - MFEM_VERIFY(bounding_box.planar, - [this]() - { - std::stringstream ss; - ss << "Boundary elements must be coplanar to define a lumped element:\n"; - auto l = bounding_box.Lengths(); - ss << "Dimensions (" << l[0] << " x " << l[1] << " x " << l[2] << ")\n"; - ss << "Area " << bounding_box.Area() << "\n"; - ss << "Volume " << bounding_box.Volume() << "\n"; - return ss.str(); - }()); - // 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); - MFEM_VERIFY(std::abs(A - bounding_box.Area()) / A < rel_tol, + MFEM_VERIFY((!bounding_box.planar || std::abs(A - bounding_box.Area()) / A < rel_tol), "Assumed bounding box area " << bounding_box.Area() << " and integrated area " << A - << " do not match: Port geometry is not rectangular!"); + << " 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; constexpr double angle_error_deg = 1.0; auto lengths = bounding_box.Lengths(); auto deviation_deg = bounding_box.Deviation(input_dir); - if ((deviation_deg[0] > angle_warning_deg && deviation_deg[1] > angle_warning_deg) || - std::isnan(deviation_deg[0]) || std::isnan(deviation_deg[1])) + if (!std::any_of(deviation_deg.begin(), deviation_deg.end(), + [](double x) { return x < angle_warning_deg; })) { auto normal_0 = bounding_box.normals[0]; for (auto &x : normal_0) @@ -107,21 +95,36 @@ class UniformElementData : public LumpedElementData { x /= lengths[1]; } - Mpi::Warning("User specified direction {} does not align with either bounding box " - "axis up to {} degrees.\n" - "Axis 1: {} ({} degrees)\nAxis 2: {} ({} degrees)!\n", - input_dir, angle_warning_deg, normal_0, deviation_deg[0], normal_1, - deviation_deg[1]); + auto normal_2 = bounding_box.normals[2]; + for (auto &x : normal_2) + { + x /= lengths[2]; + } + Mpi::Warning( + "User specified direction {:.3e} does not align with either bounding box " + "axis up to {:.3e} degrees.\n" + "Axis 1: {:.3e} ({:.3e} degrees)\nAxis 2: {:.3e} ({:.3e} degrees)\nAxis 3: {:.3e} ({:.3e} degrees)!\n", + input_dir, angle_warning_deg, normal_0, deviation_deg[0], normal_1, + deviation_deg[1], normal_2, deviation_deg[2]); } - MFEM_VERIFY(deviation_deg[0] <= angle_error_deg || deviation_deg[1] <= angle_error_deg, + 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] << ' ' << angle_error_deg - << '!'); + << deviation_deg[0] << ' ' << deviation_deg[1] << ' ' + << deviation_deg[2] << " tolerance " << angle_error_deg << '!'); std::copy(input_dir.begin(), input_dir.end(), direction.begin()); direction /= direction.Norml2(); // Compute the length from the most aligned normal direction. - l = lengths[deviation_deg[0] < deviation_deg[1] ? 0 : 1]; + l = lengths[std::distance( + deviation_deg.begin(), + std::min_element(deviation_deg.begin(), deviation_deg.end()))]; + + MFEM_ASSERT( + (l - mesh::GetDirectionalExtent(*fespace.GetParMesh(), marker, true, input_dir)) / + l < + 1e-3, + "Bounding box discovered length should match projected length"); w = A / l; } diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index b6619c551..019de7fe9 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -736,6 +736,7 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // useful for checking pointers later. 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 n_1 = (v_111 - v_000).normalized(); @@ -770,20 +771,27 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, .norm(); }; - // Filter the list of vertices to remove those that are within the plane. Planar - // detection is done using a tolerance given floating point precision. + // 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); })); + + constexpr double rel_tol = 1e-3; + box.planar = max_distance < (rel_tol * (v_111 - v_000).norm()); + + // Given numerical tolerance, collect other points with an almost matching distance. std::vector vertices_out_of_plane; - const double planar_tolerance = 1.0e-3 * (v_111 - v_000).norm(); - std::copy_if(vertices.begin(), vertices.end(), - std::back_inserter(vertices_out_of_plane), - [OutOfPlaneDistance, planar_tolerance](const auto &v) - { return OutOfPlaneDistance(v) > planar_tolerance; }); - box.planar = vertices_out_of_plane.empty(); + 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; }); // Given candidates t_0 and t_1, the closer to origin defines v_001. const auto &t_1 = box.planar ? t_0 - : *std::max_element(vertices_out_of_plane.begin(), + : *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 @@ -793,8 +801,7 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // Compute the center as halfway along the main diagonal. Vector3dMap(box.center.data()) = 0.5 * (v_000 + v_111); - // The length in each direction is then given by traversing the edges of the cuboid in - // turn. + // The length in each direction is given by traversing the edges of the cuboid in turn. Vector3dMap(box.normals[0].data()) = 0.5 * (v_001 - v_000); Vector3dMap(box.normals[1].data()) = 0.5 * (v_011 - v_001); Vector3dMap(box.normals[2].data()) = 0.5 * (v_111 - v_011); @@ -907,8 +914,45 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, return ball; } +double LengthFromPointCloud(MPI_Comm comm, const std::vector &vertices, + int dominant_rank, const std::array &dir) +{ + double length; + if (dominant_rank == Mpi::Rank(comm)) + { + MFEM_VERIFY(dir.size() == 3, "Direction must be a 3 vector!\n"); + 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); + return length; +} + } // namespace +double GetDirectionalExtent(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 GetDirectionalExtent(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 GetDirectionalExtent(mesh, marker, bdr, dir); +} + BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr) { std::vector vertices; diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 41977deaf..1b64951ef 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -112,6 +112,12 @@ struct BoundingBall BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr); BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr); +// Helper function for computing the direction aligned length of a marked group. +double GetDirectionalExtent(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, + const std::array &dir); +double GetDirectionalExtent(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(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr); From 8006ee82eca0944485deb12062e9e3da247bbc6b Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Tue, 29 Aug 2023 09:54:24 -0400 Subject: [PATCH 4/5] Renaming and style fixes --- palace/fem/lumpedelement.hpp | 18 +++++++++--------- palace/utils/geodata.cpp | 13 +++++-------- palace/utils/geodata.hpp | 8 ++++---- 3 files changed, 18 insertions(+), 21 deletions(-) diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index 779d6b5d8..aac2b2146 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -70,10 +70,10 @@ class UniformElementData : public LumpedElementData { // 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; + constexpr double rel_tol = 1.0e-3; double A = GetArea(fespace); - MFEM_VERIFY((!bounding_box.planar || std::abs(A - bounding_box.Area()) / A < rel_tol), - "Assumed bounding box area " + MFEM_VERIFY((!bounding_box.planar || (std::abs(A - bounding_box.Area()) / A < rel_tol)), + "Discovered bounding box area " << bounding_box.Area() << " and integrated area " << A << " do not match: Planar port geometry is not a quadrilateral!"); @@ -82,7 +82,7 @@ class UniformElementData : public LumpedElementData constexpr double angle_error_deg = 1.0; auto lengths = bounding_box.Lengths(); auto deviation_deg = bounding_box.Deviation(input_dir); - if (!std::any_of(deviation_deg.begin(), deviation_deg.end(), + if (std::none_of(deviation_deg.begin(), deviation_deg.end(), [](double x) { return x < angle_warning_deg; })) { auto normal_0 = bounding_box.normals[0]; @@ -103,12 +103,13 @@ class UniformElementData : public LumpedElementData Mpi::Warning( "User specified direction {:.3e} does not align with either bounding box " "axis up to {:.3e} degrees.\n" - "Axis 1: {:.3e} ({:.3e} degrees)\nAxis 2: {:.3e} ({:.3e} degrees)\nAxis 3: {:.3e} ({:.3e} degrees)!\n", + "Axis 1: {:.3e} ({:.3e} degrees)\nAxis 2: {:.3e} ({:.3e} degrees)\nAxis 3: " + "{:.3e} ({:.3e} degrees)!\n", input_dir, angle_warning_deg, normal_0, deviation_deg[0], normal_1, deviation_deg[1], normal_2, deviation_deg[2]); } MFEM_VERIFY(std::any_of(deviation_deg.begin(), deviation_deg.end(), - [](double x) { return x < angle_error_deg; }), + [angle_error_deg](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 << '!'); @@ -121,9 +122,8 @@ class UniformElementData : public LumpedElementData std::min_element(deviation_deg.begin(), deviation_deg.end()))]; MFEM_ASSERT( - (l - mesh::GetDirectionalExtent(*fespace.GetParMesh(), marker, true, input_dir)) / - l < - 1e-3, + (l - mesh::GetProjectedLength(*fespace.GetParMesh(), marker, true, input_dir)) / l < + rel_tol, "Bounding box discovered length should match projected length"); w = A / l; } diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 019de7fe9..711dc5e62 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -763,8 +763,6 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // 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_3 = n_1.cross(n_2).normalized(); - 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) @@ -920,7 +918,6 @@ double LengthFromPointCloud(MPI_Comm comm, const std::vector &v double length; if (dominant_rank == Mpi::Rank(comm)) { - MFEM_VERIFY(dir.size() == 3, "Direction must be a 3 vector!\n"); CVector3dMap direction(dir.data()); auto Dot = [&](const auto &x, const auto &y) @@ -936,21 +933,21 @@ double LengthFromPointCloud(MPI_Comm comm, const std::vector &v } // namespace -double GetDirectionalExtent(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, - const std::array &dir) +double GetProjectedLength(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 GetDirectionalExtent(mfem::ParMesh &mesh, int attr, bool bdr, - const std::array &dir) +double GetProjectedLength(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 GetDirectionalExtent(mesh, marker, bdr, dir); + return GetProjectedLength(mesh, marker, bdr, dir); } BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr) diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 1b64951ef..099e10391 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -113,10 +113,10 @@ BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array &marker, BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr); // Helper function for computing the direction aligned length of a marked group. -double GetDirectionalExtent(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, - const std::array &dir); -double GetDirectionalExtent(mfem::ParMesh &mesh, int attr, bool bdr, - const std::array &dir); +double GetProjectedLength(mfem::ParMesh &mesh, const mfem::Array &marker, bool bdr, + const std::array &dir); +double GetProjectedLength(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. From eca90d25b4f996495f0e7a1c0a0797d71f9d404c Mon Sep 17 00:00:00 2001 From: Hugh Carson Date: Tue, 29 Aug 2023 19:45:09 -0400 Subject: [PATCH 5/5] Tighten relative tolerance, fix warning messages --- palace/fem/lumpedelement.hpp | 6 +++--- palace/utils/geodata.cpp | 11 +++++------ 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/palace/fem/lumpedelement.hpp b/palace/fem/lumpedelement.hpp index aac2b2146..64101cb17 100644 --- a/palace/fem/lumpedelement.hpp +++ b/palace/fem/lumpedelement.hpp @@ -70,7 +70,7 @@ class UniformElementData : public LumpedElementData { // 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-3; + constexpr double rel_tol = 1.0e-6; double A = GetArea(fespace); MFEM_VERIFY((!bounding_box.planar || (std::abs(A - bounding_box.Area()) / A < rel_tol)), "Discovered bounding box area " @@ -102,7 +102,7 @@ class UniformElementData : public LumpedElementData } Mpi::Warning( "User specified direction {:.3e} does not align with either bounding box " - "axis up to {:.3e} degrees.\n" + "axis up to {:.3e} degrees!\n" "Axis 1: {:.3e} ({:.3e} degrees)\nAxis 2: {:.3e} ({:.3e} degrees)\nAxis 3: " "{:.3e} ({:.3e} degrees)!\n", input_dir, angle_warning_deg, normal_0, deviation_deg[0], normal_1, @@ -124,7 +124,7 @@ class UniformElementData : public LumpedElementData MFEM_ASSERT( (l - mesh::GetProjectedLength(*fespace.GetParMesh(), marker, true, input_dir)) / l < rel_tol, - "Bounding box discovered length should match projected length"); + "Bounding box discovered length should match projected length!"); w = A / l; } diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 711dc5e62..716304af5 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -736,7 +736,6 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, // useful for checking pointers later. 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 n_1 = (v_111 - v_000).normalized(); @@ -775,7 +774,7 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, [OutOfPlaneDistance](const auto &x, const auto &y) { return OutOfPlaneDistance(x) < OutOfPlaneDistance(y); })); - constexpr double rel_tol = 1e-3; + constexpr double rel_tol = 1e-6; box.planar = max_distance < (rel_tol * (v_111 - v_000).norm()); // Given numerical tolerance, collect other points with an almost matching distance. @@ -870,9 +869,9 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, vertices.begin(), vertices.end(), [&delta, PerpendicularDistance](const auto &x, const auto &y) { return PerpendicularDistance({delta}, x) < PerpendicularDistance({delta}, y); }); - constexpr double radius_tol = 1.0e-3; + constexpr double rel_tol = 1.0e-6; MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <= - radius_tol * ball.radius, + rel_tol * ball.radius, "Furthest point perpendicular must be on the exterior of the ball: " << PerpendicularDistance({delta}, perp) << " vs. " << ball.radius << "!"); @@ -892,12 +891,12 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm, }); ball.planar = - PerpendicularDistance({delta, n_radial}, out_of_plane) < radius_tol * ball.radius; + 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) <= - radius_tol * ball.radius, + rel_tol * ball.radius, "Furthest point perpendicular must be on the exterior of the sphere!"); Vector3dMap(ball.planar_normal.data()) *= 0; }