Fix bounding ball calculation for coaxial lumped ports

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
 

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<double, 3> &input_dir,
   const mfem::ParMesh &mesh = *fespace.GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> 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);
 
   // 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((!bounding_box.planar || (std::abs(A - bounding_box.Area()) / A < rel_tol)),
+  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!");
+                  << " 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<double, 3> &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,34 +107,37 @@ UniformElementData::GetModeCoefficient(double coef) const
       attr_list, source);
 }
 
-CoaxialElementData::CoaxialElementData(const std::array<double, 3> &direction,
+CoaxialElementData::CoaxialElementData(const std::array<double, 3> &input_dir,
                                        const mfem::Array<int> &attr_list,
                                        mfem::ParFiniteElementSpace &fespace)
   : LumpedElementData(attr_list)
 {
   const mfem::ParMesh &mesh = *fespace.GetParMesh();
   int bdr_attr_max = mesh.bdr_attributes.Size() ? mesh.bdr_attributes.Max() : 0;
   mfem::Array<int> attr_marker = mesh::AttrToMarker(bdr_attr_max, attr_list);
-  bounding_ball = mesh::GetBoundingBall(mesh, attr_marker, true);
+  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!");
-  sign = (direction[0] > 0);
+
+  // Direction of the excitation as +/-r̂.
+  direction = (input_dir[0] > 0);
+  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.
+  r_outer = 0.5 * bounding_ball.Lengths()[0];
   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);
+  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<mfem::VectorCoefficient>
 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;

palace/fem/lumpedelement.hpp

@@ -6,7 +6,6 @@
 
 #include <memory>
 #include <mfem.hpp>
-#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<double, 3> &direction,
+  CoaxialElementData(const std::array<double, 3> &input_dir,
                      const mfem::Array<int> &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<mfem::VectorCoefficient>