Change Cartesian to AxisAligned, style fixes, update CHANGELOG, vario…

…us minor comment fixes

CHANGELOG.md

@@ -16,7 +16,9 @@ The format of this changelog is based on
   - Added Apptainer/Singularity container build definition for Palace.
   - Added support for non axis aligned lumped ports and current sources. Key
     words `"X"`, `"Y"`, `"Z"` and `"R"`, with optional prefix `"+"` or `"-"` still work,
-    but now alternative directions can be specified as vectors with 3 components.
+    but now alternative directions can be specified as vectors with 3
+    components. Users will be warned, an ultimately errored, if the specified
+    directions do not agree with axes directions discovered from the geometry.
 
 ## [0.11.2] - 2023-07-14
 

examples/rings/mesh/mesh.jl

@@ -4,7 +4,7 @@
 using Gmsh: gmsh
 using LinearAlgebra
 
-function generate_ring_mesh(gui::Bool = false)
+function generate_ring_mesh(gui::Bool=false)
     kernel = gmsh.model.occ
 
     gmsh.initialize()
@@ -113,8 +113,8 @@ function generate_ring_mesh(gui::Bool = false)
     )
 
     # Apply a rotation transformation to all entities in the model
-    rc = [0.0,0.0,0.0]
-    ra = [1.0,1.0,1.0]
+    rc = [0.0, 0.0, 0.0]
+    ra = [1.0, 1.0, 1.0]
     θ = pi / 3
     ra ./= norm(ra)
     kernel.rotate(kernel.getEntities(), rc[1], rc[2], rc[3], ra[1], ra[2], ra[3], θ)
@@ -137,7 +137,6 @@ function generate_ring_mesh(gui::Bool = false)
     inner_gap_group = gmsh.model.addPhysicalGroup(2, [inner_gap], -1, "hole_inner")
     outer_gap_group = gmsh.model.addPhysicalGroup(2, [outer_gap], -1, "hole_outer")
 
-
     # Generate mesh
     gmsh.option.setNumber("Mesh.MeshSizeMin", l_ring)
     gmsh.option.setNumber("Mesh.MeshSizeMax", l_farfield)
@@ -212,6 +211,5 @@ function generate_ring_mesh(gui::Bool = false)
         gmsh.fltk.run()
     end
 
-    gmsh.finalize()
-
-end
+    return gmsh.finalize()
+end

examples/rings/rings.json

@@ -44,12 +44,12 @@
       {
         "Index": 1,
         "Attributes": [4],  // Inner ring
-        "Direction": [ -0.33333333333333337, 0.6666666666666667, 0.6666666666666666] // "Y" rotated pi/3 around (1,1,1)
+        "Direction": [ -0.33333333333333337, 0.6666666666666667, 0.6666666666666666] // "+Y" rotated pi/3 around (1,1,1)
       },
       {
         "Index": 2,
         "Attributes": [5],  // Outer ring
-        "Direction": [ -0.33333333333333337, 0.6666666666666667, 0.6666666666666666] // "Y" rotated pi/3 around (1,1,1)
+        "Direction": [ -0.33333333333333337, 0.6666666666666667, 0.6666666666666666] // "+Y" rotated pi/3 around (1,1,1)
       }
     ],
     "Postprocessing":  // Inductance from flux instead of energy

palace/fem/lumpedelement.hpp

@@ -8,8 +8,8 @@
 #include <string>
 #include <mfem.hpp>
 #include "fem/integrator.hpp"
-#include "utils/geodata.hpp"
 #include "utils/communication.hpp"
+#include "utils/geodata.hpp"
 
 namespace palace
 {
@@ -81,11 +81,8 @@ class UniformElementData : public LumpedElementData
 
     // Check the user specified direction aligns with an axis direction
     auto Dot = [](const auto &x_1, const auto &x_2)
-    {
-      return x_1[0] * x_2[0] + x_1[1] * x_2[1] + x_1[2] * x_2[2];
-    };
-    auto AbsAngle = [Dot](const auto &x_1, const auto &x_2)
-    {
+    { return x_1[0] * x_2[0] + x_1[1] * x_2[1] + x_1[2] * x_2[2]; };
+    auto AbsAngle = [Dot](const auto &x_1, const auto &x_2) {
       return std::acos(std::abs(Dot(x_1, x_2)) / std::sqrt(Dot(x_1, x_1) * Dot(x_2, x_2)));
     };
 
@@ -97,21 +94,29 @@ class UniformElementData : public LumpedElementData
     if (deviation_0_deg > angle_warning_deg && deviation_1_deg > angle_warning_deg)
     {
       auto normalized_0 = bounding_box.normals[0];
-      for (auto &x : normalized_0) { x /= std::sqrt(Dot(bounding_box.normals[0], bounding_box.normals[0])); }
+      for (auto &x : normalized_0)
+      {
+        x /= std::sqrt(Dot(bounding_box.normals[0], bounding_box.normals[0]));
+      }
       auto normalized_1 = bounding_box.normals[1];
-      for (auto &x : normalized_1) { x /= std::sqrt(Dot(bounding_box.normals[1], bounding_box.normals[1])); }
-      Mpi::Warning("User specified direction {} does not align with either bounding box axis up to {} degrees.\n"
-      "Axis 1: {} ({} degrees)\nAxis 2: {} ({} degrees)",
-        input_dir, angle_warning_deg, normalized_0, deviation_0_deg, normalized_1, deviation_1_deg);
+      for (auto &x : normalized_1)
+      {
+        x /= std::sqrt(Dot(bounding_box.normals[1], bounding_box.normals[1]));
+      }
+      Mpi::Warning("User specified direction {} does not align with either bounding box "
+                   "axis up to {} degrees.\n"
+                   "Axis 1: {} ({} degrees)\nAxis 2: {} ({} degrees)",
+                   input_dir, angle_warning_deg, normalized_0, deviation_0_deg,
+                   normalized_1, deviation_1_deg);
     }
 
     MFEM_VERIFY(deviation_0_deg <= angle_error_deg || deviation_1_deg <= angle_error_deg,
-      "Specified direction does not align sufficiently with bounding box axes");
+                "Specified direction does not align sufficiently with bounding box axes");
 
     // Compute the length from the most aligned normal direction.
     l = 2 * std::sqrt(deviation_0_deg < deviation_1_deg
-                        ? Dot(bounding_box.normals[0], bounding_box.normals[0])
-                        : Dot(bounding_box.normals[1], bounding_box.normals[1]));
+                          ? Dot(bounding_box.normals[0], bounding_box.normals[0])
+                          : Dot(bounding_box.normals[1], bounding_box.normals[1]));
     w = bounding_box.Area() / l;
   }
 

palace/utils/communication.hpp

@@ -133,7 +133,7 @@ struct ValueAndRank
 {
   // The reduced value from the MINLOC or MAXLOC operation.
   T value;
-  // The processor rank that has the value/
+  // The processor rank that has the value.
   int rank;
   // Convenience decay operators to allow use as the underlying value.
   inline operator const T &() const { return value; }

palace/utils/configfile.cpp

@@ -815,7 +815,6 @@ auto ParseDataNode(json &j, const std::string &key_word, bool is_port = true)
     }
   }
 
-
   return node;
 }
 }  // namespace

palace/utils/configfile.hpp

@@ -360,15 +360,10 @@ enum class CoordinateSystem
 // Exactly one of either the direction or the normal direction must be defined.
 struct DataNode
 {
-  // Optional string defining source excitation field direction. Options are "X", "Y", "Z",
-  // or "R", preceeded with a "+" or "-" for the direction. The first three options are for
-  // uniform lumped ports and the last is for a coaxial lumped port (radial excitation).
-  // std::string direction = "";
-
   // Vector defining the normal direction for this port. In a Cartesian system
   // with  "X", "Y", and "Z" mapping to (1,0,0), (0,1,0), and (0,0,1)
   // respectively. Any user specified value is normalized to a tolerance of
-  // 1e-6, and if normalization is needed, printed to the terminal.
+  // 1e-5, and if normalization is needed, printed to the terminal.
   std::array<double, 3> normal{{0.0, 0.0, 0.0}};
 
   // Coordinate system that the normal vector is expressed in
@@ -463,12 +458,6 @@ struct CapacitancePostData : public internal::DataMap<CapacitanceData>
   void SetUp(json &postpro);
 };
 
-// using InductanceData = DataNode;
-// struct InductanceData : public DataNode
-// {
-//   using DataNode::DataNode;
-// };
-
 struct InductancePostData : public internal::DataMap<DataNode>
 {
 public:

palace/utils/geodata.cpp

@@ -849,14 +849,11 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
     n_radial[1] /= radius_2;
     n_radial[2] /= radius_2;
 
-    // Compute a perpendicular to the circle using the cross product
-    ball.planar_normal[0] = delta[1] * n_radial[2] - delta[2] * n_radial[1];
-    ball.planar_normal[1] = delta[2] * n_radial[0] - delta[0] * n_radial[2];
-    ball.planar_normal[2] = delta[0] * n_radial[1] - delta[1] * n_radial[0];
+    // Compute a perpendicular to the circle using the cross product.
+    ball.planar_normal = Cross(delta, n_radial);
 
     // Compute the point furthest out of the plane discovered. If below
     // tolerance, this means the circle is 2D.
-
     const auto &out_of_plane = *std::max_element(
         vertices.begin(), vertices.end(),
         [&delta, &n_radial, PerpendicularDistance](const auto &x, const auto &y)
@@ -866,15 +863,12 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
         });
 
     constexpr double planar_tolerance = 1e-9;
-    if (PerpendicularDistance({delta, n_radial}, out_of_plane) < planar_tolerance)
+    if (PerpendicularDistance({delta, n_radial}, out_of_plane) > planar_tolerance)
     {
-      // The points are functionally coplanar, zero out the normal
+      // The points are not functionally coplanar, zero out the normal
       ball.planar_normal[0] *= 0;
       ball.planar_normal[1] *= 0;
       ball.planar_normal[2] *= 0;
-    }
-    else
-    {
       MFEM_VERIFY(PerpendicularDistance({delta, n_radial}, out_of_plane) <= ball.radius,
                   "A point perpendicular must be contained in the ball");
     }
@@ -939,17 +933,17 @@ BoundingBall GetBoundingBall(mfem::ParMesh &mesh, int attr, bool bdr)
   return GetBoundingBall(mesh, marker, bdr);
 }
 
-void GetCartesianBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                             mfem::Vector &max)
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
+                               mfem::Vector &max)
 {
   mfem::Array<int> marker(bdr ? mesh.bdr_attributes.Max() : mesh.attributes.Max());
   marker = 0;
   marker[attr - 1] = 1;
-  GetCartesianBoundingBox(mesh, marker, bdr, min, max);
+  GetAxisAlignedBoundingBox(mesh, marker, bdr, min, max);
 }
 
-void GetCartesianBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
-                             mfem::Vector &min, mfem::Vector &max)
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                               bool bdr, mfem::Vector &min, mfem::Vector &max)
 {
   int dim = mesh.SpaceDimension();
   min.SetSize(dim);

palace/utils/geodata.hpp

@@ -52,13 +52,13 @@ void AttrToMarker(int max_attr, const mfem::Array<int> &attrs, mfem::Array<int>
 void AttrToMarker(int max_attr, const std::vector<int> &attrs, mfem::Array<int> &marker);
 
 // Helper function to construct the bounding box for all elements with the given attribute.
-void GetCartesianBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
-                             mfem::Vector &max);
-void GetCartesianBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr,
-                             mfem::Vector &min, mfem::Vector &max);
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr, mfem::Vector &min,
+                               mfem::Vector &max);
+void GetAxisAlignedBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker,
+                               bool bdr, mfem::Vector &min, mfem::Vector &max);
 
 // Struct describing a bounding box in terms of the center and face normals. The
-// normals specify the distance from the center of the box.
+// normals specify the direction from the center of the box.
 struct BoundingBox
 {
   // The central point of the bounding box.
@@ -81,19 +81,21 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm,
 BoundingBox GetBoundingBox(mfem::ParMesh &mesh, const mfem::Array<int> &marker, bool bdr);
 BoundingBox GetBoundingBox(mfem::ParMesh &mesh, int attr, bool bdr);
 
+// 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
+  // The centroid of the ball.
   std::array<double, 3> center;
-  // The radius of the ball from the 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<double, 3> planar_normal;
 };
 
-// Helper for computing a diameter from a point cloud, assumes that contains a
-// set of points on a circle, and uses lexicographic min and max.
+// Helper for computing a diameter from a point cloud. The cloud is assumed to
+// contain a set of points on a circle.
 BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
                                         std::vector<std::array<double, 3>> vertices);
 // Given a mesh and a marker, compute the diameter of a bounding circle/sphere,