Further style fixes prior to merge, including column-width for docs a…

…nd comments and allowing GlobalMinLoc/MaxLoc to use arbitrary type for the second argument as supported by the available MPI types

docs/src/config/boundaries.md

@@ -278,20 +278,19 @@ with
 boundary. If this port is to be a multielement lumped port with more than a single lumped
 element, use the `"Elements"` array described below.
 
-`"Direction" [None]` :  Direction to define the polarization direction of the
-port field mode on this lumped port boundary. Axis aligned lumped ports can be
-specified using keywords: `"+X"`, `"-X"`, `"+Y"`, `"-Y"`, `"+Z"`, `"-Z"`, while
-coaxial lumped ports can be specified using `"+R"`, `"-R"`. The direction can
-alternatively be specified as a normalized array of three values, for example
-`[0, 1, 0]`. If a vector direction is specified, the `"CoordinateSystem"` value
-specifies the coordinate system it is expressed in. If this port is to be a
-multielement lumped port with more than a single lumped element, use the
-`"Elements"` array described below.
+`"Direction" [None]` :  Direction to define the polarization direction of the port field
+mode on this lumped port boundary. Axis aligned lumped ports can be specified using
+keywords: `"+X"`, `"-X"`, `"+Y"`, `"-Y"`, `"+Z"`, `"-Z"`, while coaxial lumped ports can be
+specified using `"+R"`, `"-R"`. The direction can alternatively be specified as a
+normalized array of three values, for example `[0, 1, 0]`. If a vector direction is
+specified, the `"CoordinateSystem"` value specifies the coordinate system it is expressed
+in. If this port is to be a multielement lumped port with more than a single lumped
+element, use the `"Elements"` array described below.
 
-`"CoordinateSystem" ["Cartesian"]` : Coordinate system used to express the
-`"Direction"` vector, the options are `"Cartesian"` and `"Cylindrical"`. If a
-keyword argument is used for `"Direction"` this value is ignored, and the
-appropriate coordinate system is used instead.
+`"CoordinateSystem" ["Cartesian"]` : Coordinate system used to express the `"Direction"`
+vector, the options are `"Cartesian"` and `"Cylindrical"`. If a keyword argument is used
+for `"Direction"` this value is ignored, and the appropriate coordinate system is used
+instead.
 
 `"Excitation" [false]` :  Turns on or off port excitation for this lumped port boundary for
 driven or transient simulation types.
@@ -330,10 +329,10 @@ should not be combined with the `"Direction"` field described above. Each elemen
 multielement lumped port can be described by its own unique direction, which is specified
 here. The elements of a multielement port add in parallel.
 
-`"Elements"[]["CoordinateSystem"] ["Cartesian"]` :  This option is for multielement
-lumped ports and should not be combined with the `"CoordinateSystem"` field
-described above. Each element of a multielement lumped port can be described by
-its own unique direction, and corresponding coordinate system.
+`"Elements"[]["CoordinateSystem"] ["Cartesian"]` :  This option is for multielement lumped
+ports and should not be combined with the `"CoordinateSystem"` field described above. Each
+element of a multielement lumped port can be described by its own unique direction, and
+corresponding coordinate system.
 
 ## `boundaries["WavePort"]`
 
@@ -423,9 +422,9 @@ single lumped element, use the `"Elements"` array described below.
 
 `"CoordinateSystem" ["Cartesian"]` :  Defines the coordinate system for the source current
 direction for this surface current boundary. The available options are the same as under
-[`config["Boundaries"]["LumpedPort"]["CoordinateSystem"]`](#boundaries%5B%22LumpedPort%22%5D). If
-this source is to be a multielement source which distributes the source across more than a
-single lumped element, use the `"Elements"` array described below.
+[`config["Boundaries"]["LumpedPort"]["CoordinateSystem"]`](#boundaries%5B%22LumpedPort%22%5D).
+If this source is to be a multielement source which distributes the source across more than
+a single lumped element, use the `"Elements"` array described below.
 
 `"Elements"[]["Attributes"] [None]` :  This option is for multielement surface current
 boundaries should not be combined with the `"Attributes"` field described above. Each
@@ -544,10 +543,10 @@ postprocessing boundary.
 
 `"Direction" [None]` :  Defines the global direction with which to orient the surface
 normals with computing the magnetic flux for this inductance postprocessing boundary. The
-available options are: `"+X"`, `"-X"`, `"+Y"`, `"-Y"`, `"+Z"`, `"-Z"`. The
-direction can alternatively be specified as a normalized array of three values,
-for example `[0, 1, 0]`. The true surface normal is used in the calculation,
-`"Direction"` is only used to ensure the correct choice of orientation of the normal.
+available options are: `"+X"`, `"-X"`, `"+Y"`, `"-Y"`, `"+Z"`, `"-Z"`. The direction can
+alternatively be specified as a normalized array of three values, for example `[0, 1, 0]`.
+The true surface normal is used in the calculation, `"Direction"` is only used to ensure
+the correct choice of orientation of the normal.
 
 ## `boundaries["Postprocessing"]["Dielectric"]`
 
@@ -592,12 +591,13 @@ use the `"Elements"` array described below.
 `"Side" [None]` :  Defines the postprocessing side when this dielectric interface is an
 internal boundary surface (and thus the electric field on the boundary is in general
 double-valued). The available options are: `"+X"`, `"-X"`, `"+Y"`, `"-Y"`, `"+Z"`, `"-Z"`.
-If the boundary is not axis-aligned, the field value is taken from the side which is
-oriented along the specified direction. If no `"Side"` is specified, the field solution is
-taken from the neighboring element with the smaller electrical permittivity, which is an
-attempt to get the field in the domain corresponding to vacuum. If the interface consists
-of multiple elements with different `"Side"` values, use the `"Elements"` array described
-below.
+The direction can alternatively be specified as a normalized array of three values, for
+example `[0, 1, 0]`. If the boundary is not axis-aligned, the field value is taken from the
+side which is oriented along the specified direction. If no `"Side"` is specified, the
+field solution is taken from the neighboring element with the smaller electrical
+permittivity, which is an attempt to get the field in the domain corresponding to vacuum.
+If the interface consists of multiple elements with different `"Side"` values, use the
+`"Elements"` array described below.
 
 `"Thickness" [None]` :  Thickness of this dielectric interface, specified in mesh length
 units.

palace/fem/lumpedelement.hpp

@@ -71,8 +71,8 @@ class UniformElementData : public LumpedElementData
     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.
+    // 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,

palace/utils/communication.hpp

@@ -127,49 +127,45 @@ inline MPI_Datatype DataType<bool>()
   return MPI_C_BOOL;
 }
 
-// Helper for computing a reduction along with an identifier, for example rank
-// of the owning process.
-template <typename T, typename U = signed int>
-struct ValueAndId
+template <typename T, typename U>
+struct ValueAndLoc
 {
-  // The value used in a reduction.
   T val;
-  // The id that corresponds to the reduced value.
-  U id;
+  U loc;
 };
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<float, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<float, signed int>>()
 {
   return MPI_FLOAT_INT;
 }
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<double, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<double, signed int>>()
 {
   return MPI_DOUBLE_INT;
 }
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<long double, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<long double, signed int>>()
 {
   return MPI_LONG_DOUBLE_INT;
 }
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<signed short, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<signed short, signed int>>()
 {
   return MPI_SHORT_INT;
 }
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<signed int, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<signed int, signed int>>()
 {
   return MPI_2INT;
 }
 
 template <>
-inline MPI_Datatype DataType<ValueAndId<signed long int, signed int>>()
+inline MPI_Datatype DataType<ValueAndLoc<signed long int, signed int>>()
 {
   return MPI_LONG_INT;
 }
@@ -266,41 +262,39 @@ class Mpi
     GlobalOp(len, buff, MPI_SUM, comm);
   }
 
-  // Global minimum with id (in-place, result is broadcast to all processes).
-  template <typename T>
-  static void GlobalMinLoc(int len, T *val, int *id, MPI_Comm comm)
+  // Global minimum with index (in-place, result is broadcast to all processes).
+  template <typename T, typename U>
+  static void GlobalMinLoc(int len, T *val, U *loc, MPI_Comm comm)
   {
-    std::vector<mpi::ValueAndId<T, int>> buffer(len);
+    std::vector<mpi::ValueAndLoc<T, U>> buffer(len);
     for (int i = 0; i < len; i++)
     {
       buffer[i].val = val[i];
-      buffer[i].id = id[i];
+      buffer[i].loc = loc[i];
     }
-
     GlobalOp(len, buffer.data(), MPI_MINLOC, comm);
     for (int i = 0; i < len; i++)
     {
       val[i] = buffer[i].val;
-      id[i] = buffer[i].id;
+      loc[i] = buffer[i].loc;
     }
   }
 
-  // Global minimum with id (in-place, result is broadcast to all processes).
-  template <typename T>
-  static void GlobalMaxLoc(int len, T *val, int *id, MPI_Comm comm)
+  // Global maximum with index (in-place, result is broadcast to all processes).
+  template <typename T, typename U>
+  static void GlobalMaxLoc(int len, T *val, U *loc, MPI_Comm comm)
   {
-    std::vector<mpi::ValueAndId<T, int>> buffer(len);
+    std::vector<mpi::ValueAndLoc<T, U>> buffer(len);
     for (int i = 0; i < len; i++)
     {
       buffer[i].val = val[i];
-      buffer[i].id = id[i];
+      buffer[i].loc = loc[i];
     }
-
     GlobalOp(len, buffer.data(), MPI_MAXLOC, comm);
     for (int i = 0; i < len; i++)
     {
       val[i] = buffer[i].val;
-      id[i] = buffer[i].id;
+      loc[i] = buffer[i].loc;
     }
   }
 

palace/utils/configfile.hpp

@@ -61,7 +61,7 @@ struct DataMap
 };
 
 // An ElementData consists of a list of attributes making up a single element of a
-// potentially multielement node, and a direction and/or a normal defining the incident
+// potentially multielement boundary, and a direction and/or a normal defining the incident
 // field. These are used for lumped ports, terminals, surface currents, and other boundary
 // postprocessing objects.
 struct ElementData

palace/utils/geodata.cpp

@@ -708,21 +708,21 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, std::vector<Eigen::Vector3d
     MFEM_VERIFY(vertices.size() >= 4,
                 "A bounding box requires a minimum of four vertices for this algorithm!");
 
-    // 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. Store references as this is
+    // useful for checking pointers later.
     const auto &v_000 = *std::min_element(vertices.begin(), vertices.end(), EigenLE);
     const auto &v_111 = *std::max_element(vertices.begin(), vertices.end(), EigenLE);
     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();
 
-    // Compute the distance from the normal axis. Note: everything has been
-    // oriented relative to v_000 == (0,0,0).
+    // 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).
+    // 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)
@@ -734,27 +734,26 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, std::vector<Eigen::Vector3d
     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.
+    // 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_3 = n_1.cross(n_2).normalized();
 
     auto OutOfPlaneComp = [&n_1, &n_2, &n_3, &origin](const auto &v_1, const auto &v_2)
     {
-      // 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.
+      // 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;
     };
 
-    // There is a degeneration if the final point is within the plane defined by
-    // (v_000, v_001, v_111).
+    // 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 &&
@@ -768,18 +767,18 @@ BoundingBox BoundingBoxFromPointCloud(MPI_Comm comm, std::vector<Eigen::Vector3d
       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.
+    // 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;
 
     // 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 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();
@@ -858,10 +857,8 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
 
     // 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.
+    // 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<Eigen::Vector3d> &deltas,
                                        const Eigen::Vector3d &vin)
     {
@@ -889,8 +886,8 @@ BoundingBall BoundingBallFromPointCloud(MPI_Comm comm,
     const Eigen::Vector3d n_radial = (perp - CVector3dMap(ball.center.data())).normalized();
     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.
+    // 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)

palace/utils/geodata.hpp

@@ -80,8 +80,7 @@ struct BoundingBox
   // Compute the lengths of each axis.
   std::array<double, 3> 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 normal directions.
   std::array<double, 3> Deviation(const std::array<double, 3> &direction) const;
 };