Fix for non-planar lumped ports

palace/fem/lumpedelement.hpp

@@ -68,25 +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,
-                "Boundary elements must be coplanar to define a 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);
-    MFEM_VERIFY(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: 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::none_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)
@@ -98,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(),
+                            [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] << ' ' << 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::GetProjectedLength(*fespace.GetParMesh(), marker, true, input_dir)) / l <
+            rel_tol,
+        "Bounding box discovered length should match projected length!");
     w = A / l;
   }
 

palace/utils/geodata.cpp

@@ -762,52 +762,46 @@ 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 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<double, 3> dist_1{(v_1 - origin).dot(n_3), (v_1 - origin).dot(n_2),
-                                         (v_1 - origin).dot(n_1)};
-      const std::array<double, 3> 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;
+    // 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-6;
+    box.planar = max_distance < (rel_tol * (v_111 - v_000).norm());
+
+    // Given numerical tolerance, collect other points with an almost matching distance.
+    std::vector<Eigen::Vector3d> 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; });
+
+    // 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 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);
 
-    // 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()) =
-        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 +869,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 rel_tol = 1.0e-6;
     MFEM_VERIFY(std::abs(PerpendicularDistance({delta}, perp) - ball.radius) <=
-                    radius_tol * ball.radius,
-                "A point perpendicular must be contained in the ball: "
+                    rel_tol * ball.radius,
+                "Furthest point perpendicular must be on the exterior of the ball: "
                     << PerpendicularDistance({delta}, perp) << " vs. " << ball.radius
                     << "!");
 
@@ -887,7 +881,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 +890,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) / ball.radius < rel_tol;
     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) <=
+                      rel_tol * ball.radius,
+                  "Furthest point perpendicular must be on the exterior of the sphere!");
       Vector3dMap(ball.planar_normal.data()) *= 0;
     }
   }
@@ -916,8 +911,44 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
   return ball;
 }
 
+double LengthFromPointCloud(MPI_Comm comm, const std::vector<Eigen::Vector3d> &vertices,
+                            int dominant_rank, const std::array<double, 3> &dir)
+{
+  double length;
+  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);
+  return length;
+}
+
 }  // namespace
 
+double GetProjectedLength(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+                          const std::array<double, 3> &dir)
+{
+  std::vector<Eigen::Vector3d> vertices;
+  int dominant_rank = CollectPointCloudOnRoot(mesh, marker, bdr, vertices);
+  return LengthFromPointCloud(mesh.GetComm(), vertices, dominant_rank, dir);
+}
+
+double GetProjectedLength(mfem::ParMesh &mesh, int attr, bool bdr,
+                          const std::array<double, 3> &dir)
+{
+  mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
+  marker = 0;
+  marker[attr - 1] = 1;
+  return GetProjectedLength(mesh, marker, bdr, dir);
+}
+
 BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr)
 {
   std::vector<Eigen::Vector3d> vertices;

palace/utils/geodata.hpp

@@ -112,6 +112,12 @@ struct BoundingBall
 BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &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 GetProjectedLength(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
+                          const std::array<double, 3> &dir);
+double GetProjectedLength(mfem::ParMesh &mesh, int attr, bool bdr,
+                          const std::array<double, 3> &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<int> &marker, bool bdr);