Merge branch 'develop' into feature/han12/cuda_tutorial_image

RELEASE-NOTES.md

@@ -42,6 +42,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 - Improves support for `axom::Array` allocated in unified and pinned memory on GPU platforms.
   Use of GPU-based operations for Arrays allocated in a unified memory space is controlled with
   a new method, `Array::setDevicePreference()`.
+- Adds `svg2contours` script to convert paths in an SVG file to an MFEM NURBS mesh
+- Quest: Adds an example to query winding numbers on an MFEM NURBS mesh
 
 ### Changed
 - `MarchingCubes` masking now uses the mask field's integer values instead of
@@ -63,6 +65,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
   external `fmt` and axom's vendored copy.
 - Turn off CMake finding dependencies on system paths.
 - `axom::Array`: trivially-copyable types with a non-trivial constructor are now initialized on the GPU.
+- SLIC no longer outputs the rank count in the `RANK` format string in parallel loggers. You can access
+  the rank count via new format option `RANK_COUNT`.
 
 ### Removed
 - Removes config option `AXOM_ENABLE_ANNOTATIONS`. Annotations are now provided by `caliper` 

src/axom/primal/geometry/BezierCurve.hpp

@@ -26,6 +26,7 @@
 
 #include <vector>
 #include <ostream>
+#include "axom/fmt.hpp"
 
 namespace axom
 {
@@ -919,4 +920,10 @@ std::ostream& operator<<(std::ostream& os, const BezierCurve<T, NDIMS>& bCurve)
 }  // namespace primal
 }  // namespace axom
 
+/// Overload to format a primal::BezierCurve using fmt
+template <typename T, int NDIMS>
+struct axom::fmt::formatter<axom::primal::BezierCurve<T, NDIMS>>
+  : ostream_formatter
+{ };
+
 #endif  // AXOM_PRIMAL_BEZIERCURVE_HPP_

src/axom/quest/examples/CMakeLists.txt

@@ -532,3 +532,15 @@ if (ENABLE_FORTRAN)
         endif()
     endif()
 endif()
+
+# Quest winding number example ------------------------------------------------
+
+if( MFEM_FOUND)
+    axom_add_executable(
+        NAME        quest_winding_number_ex
+        SOURCES     quest_winding_number.cpp
+        OUTPUT_DIR  ${EXAMPLE_OUTPUT_DIRECTORY}
+        DEPENDS_ON  axom mfem cli11 fmt
+        FOLDER      axom/quest/examples
+        )
+endif()

src/axom/quest/examples/quest_winding_number.cpp

@@ -0,0 +1,306 @@
+// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+/*!
+ * \file quest_winding_number.cpp
+ * \brief Example that computes the winding number of a grid of points
+ * against a collection of 2D parametric rational curves.
+ * Supports MFEM meshes in the cubic positive Bernstein basis or the (rational)
+ * NURBS basis.
+ */
+
+#include "axom/config.hpp"
+#include "axom/core.hpp"
+#include "axom/slic.hpp"
+#include "axom/primal.hpp"
+#include "axom/quest.hpp"
+
+#include "axom/CLI11.hpp"
+#include "axom/fmt.hpp"
+
+#include "mfem.hpp"
+
+namespace primal = axom::primal;
+using Point2D = primal::Point<double, 2>;
+using BezierCurve2D = primal::BezierCurve<double, 2>;
+using BoundingBox2D = primal::BoundingBox<double, 2>;
+
+/*!
+ * Given an mfem mesh, convert element with id \a elem_id to a (rational) BezierCurve
+ * \pre Assumes the elements of the mfem mesh are in the positive (Bernstein)
+ * basis, or in the NURBS basis
+ */
+BezierCurve2D segment_to_curve(const mfem::Mesh* mesh, int elem_id)
+{
+  const auto* fes = mesh->GetNodes()->FESpace();
+  const auto* fec = fes->FEColl();
+
+  const bool isBernstein =
+    dynamic_cast<const mfem::H1Pos_FECollection*>(fec) != nullptr;
+  const bool isNURBS =
+    dynamic_cast<const mfem::NURBSFECollection*>(fec) != nullptr;
+
+  SLIC_ERROR_IF(
+    !(isBernstein || isNURBS),
+    "MFEM mesh elements must be in either the Bernstein or NURBS basis");
+
+  const int NE = isBernstein ? mesh->GetNE() : fes->GetNURBSext()->GetNP();
+  SLIC_ERROR_IF(NE < elem_id,
+                axom::fmt::format("Mesh does not have {} elements", elem_id));
+
+  const int order =
+    isBernstein ? fes->GetOrder(elem_id) : mesh->NURBSext->GetOrders()[elem_id];
+  SLIC_ERROR_IF(order != 3,
+                axom::fmt::format(
+                  "This example currently requires the input mfem mesh to "
+                  "contain cubic elements, but the order of element {} is {}",
+                  elem_id,
+                  order));
+
+  mfem::Array<int> dofs;
+  mfem::Array<int> vdofs;
+
+  mfem::Vector dvec;
+  mfem::Vector v;
+
+  fes->GetElementDofs(elem_id, dofs);
+  fes->GetElementVDofs(elem_id, vdofs);
+  mesh->GetNodes()->GetSubVector(vdofs, v);
+
+  // Currently hard-coded for 3rd order. This can easily be extended to arbitrary order
+  axom::Array<Point2D> points(4, 4);
+  if(isBernstein)
+  {
+    points[0] = Point2D {v[0], v[0 + 4]};
+    points[1] = Point2D {v[2], v[2 + 4]};
+    points[2] = Point2D {v[3], v[3 + 4]};
+    points[3] = Point2D {v[1], v[1 + 4]};
+
+    return BezierCurve2D(points, fec->GetOrder());
+  }
+  else  // isNURBS
+  {
+    // temporary assumption is that there are no interior knots
+    // i.e. the NURBS curve is essentially a rational Bezier curve
+
+    points[0] = Point2D {v[0], v[0 + 4]};
+    points[1] = Point2D {v[1], v[1 + 4]};
+    points[2] = Point2D {v[2], v[2 + 4]};
+    points[3] = Point2D {v[3], v[3 + 4]};
+
+    fes->GetNURBSext()->GetWeights().GetSubVector(dofs, dvec);
+    axom::Array<double> weights {dvec[0], dvec[1], dvec[2], dvec[3]};
+
+    return BezierCurve2D(points, weights, fec->GetOrder());
+  }
+}
+
+bool check_mesh_valid(const mfem::Mesh* mesh)
+{
+  const auto* fes = mesh->GetNodes()->FESpace();
+  if(fes == nullptr)
+  {
+    SLIC_WARNING("MFEM mesh finite element space was null");
+    return false;
+  }
+
+  const auto* fec = fes->FEColl();
+  if(fec == nullptr)
+  {
+    SLIC_WARNING("MFEM mesh finite element collection was null");
+    return false;
+  }
+
+  const bool isBernstein =
+    dynamic_cast<const mfem::H1Pos_FECollection*>(fec) != nullptr;
+  const bool isNURBS =
+    dynamic_cast<const mfem::NURBSFECollection*>(fec) != nullptr;
+  const bool isValidFEC = isBernstein || isNURBS;
+
+  // TODO: Convert from Lagrange to Bernstein, if/when necessary
+  if(!isValidFEC)
+  {
+    SLIC_WARNING(
+      "Example only currently supports 1D NURBS meshes "
+      "or meshes with nodes in the Bernstein basis");
+    return false;
+  }
+
+  if(fes->GetVDim() != 2)
+  {
+    SLIC_WARNING("Example only currently supports 2D meshes");
+    return false;
+  }
+
+  const int NE = isBernstein ? mesh->GetNE() : fes->GetNURBSext()->GetNP();
+  int order = -1;
+  if(isBernstein)
+  {
+    order = NE > 0 ? fes->GetOrder(0) : 3;
+  }
+  else  // isNURBS
+  {
+    //SLIC_INFO("nurbsext order :" << mesh->NURBSext->GetOrder());
+    order = NE > 0 ? mesh->NURBSext->GetOrders()[0] : 3;
+  }
+
+  if(order != 3)
+  {
+    SLIC_WARNING(axom::fmt::format(
+      "This example currently requires the input mfem mesh to contain cubic "
+      "elements, but the provided mesh has order {}",
+      order));
+    return false;
+  }
+
+  return true;
+}
+
+int main(int argc, char** argv)
+{
+  axom::slic::SimpleLogger raii_logger;
+
+  axom::CLI::App app {
+    "Load mesh containing collection of curves"
+    " and optionally generate a query mesh of winding numbers."};
+
+  std::string inputFile;
+  std::string outputPrefix = {"winding"};
+
+  bool verbose {false};
+
+  // Query mesh parameters
+  std::vector<double> boxMins;
+  std::vector<double> boxMaxs;
+  std::vector<int> boxResolution;
+  int queryOrder {1};
+
+  app.add_option("-i,--input", inputFile)
+    ->description("MFEM mesh containing contours (1D segments)")
+    ->required()
+    ->check(axom::CLI::ExistingFile);
+
+  app.add_option("-o,--output-prefix", outputPrefix)
+    ->description(
+      "Prefix for output 2D query mesh (in MFEM format) mesh containing "
+      "winding number calculations")
+    ->capture_default_str();
+
+  app.add_flag("-v,--verbose", verbose, "verbose output")->capture_default_str();
+
+  auto* query_mesh_subcommand =
+    app.add_subcommand("query_mesh")
+      ->description("Options for setting up a query mesh")
+      ->fallthrough();
+  query_mesh_subcommand->add_option("--min", boxMins)
+    ->description("Min bounds for box mesh (x,y)")
+    ->expected(2)
+    ->required();
+  query_mesh_subcommand->add_option("--max", boxMaxs)
+    ->description("Max bounds for box mesh (x,y)")
+    ->expected(2)
+    ->required();
+  query_mesh_subcommand->add_option("--res", boxResolution)
+    ->description("Resolution of the box mesh (i,j)")
+    ->expected(2)
+    ->required();
+  query_mesh_subcommand->add_option("--order", queryOrder)
+    ->description("polynomial order of the query mesh")
+    ->check(axom::CLI::PositiveNumber);
+
+  CLI11_PARSE(app, argc, argv);
+
+  mfem::Mesh mesh(inputFile);
+  SLIC_INFO(
+    axom::fmt::format("Curve mesh has a topological dimension of {}d, "
+                      "has {} vertices and {} elements",
+                      mesh.Dimension(),
+                      mesh.GetNV(),
+                      mesh.GetNE()));
+
+  if(!check_mesh_valid(&mesh))
+  {
+    return 1;
+  }
+
+  axom::Array<int> segments;
+  axom::Array<BezierCurve2D> curves;
+
+  // Loop through mesh elements, retaining the (curved) 1D segments
+  for(int i = 0; i < mesh.GetNE(); ++i)
+  {
+    auto* el = mesh.GetElement(i);
+    if(el->GetGeometryType() == mfem::Geometry::SEGMENT)
+    {
+      segments.push_back(i);
+    }
+  }
+
+  // Extract the curves and compute their bounding boxes along the way
+  BoundingBox2D bbox;
+  for(int i = 0; i < segments.size(); ++i)
+  {
+    auto curve = segment_to_curve(&mesh, i);
+    SLIC_INFO_IF(verbose, axom::fmt::format("Element {}: {}", i, curve));
+
+    bbox.addBox(curve.boundingBox());
+
+    curves.emplace_back(std::move(curve));
+  }
+
+  SLIC_INFO(axom::fmt::format("Curve mesh bounding box: {}", bbox));
+
+  // Early return if user didn't set up a query mesh
+  if(boxMins.empty())
+  {
+    return 0;
+  }
+
+  // Generate a Cartesian (high order) mesh for the query points
+  const auto query_res = primal::NumericArray<int, 2>(boxResolution.data());
+  const auto query_box =
+    BoundingBox2D(Point2D(boxMins.data()), Point2D(boxMaxs.data()));
+
+  auto query_mesh = std::unique_ptr<mfem::Mesh>(
+    axom::quest::util::make_cartesian_mfem_mesh_2D(query_box,
+                                                   query_res,
+                                                   queryOrder));
+  auto fec = mfem::H1_FECollection(queryOrder, 2);
+  auto fes = mfem::FiniteElementSpace(query_mesh.get(), &fec, 1);
+  auto winding = mfem::GridFunction(&fes);
+  auto inout = mfem::GridFunction(&fes);
+  auto nodes_fes = query_mesh->GetNodalFESpace();
+
+  // Query the winding numbers at each degree of freedom (DoF) of the query mesh.
+  // The loop below independently checks (and adaptively refines) every curve for each query points.
+  // A more efficient algorithm can de defined that caches the refined curves to avoid
+  // extra refinements. We will add this in a follow-up PR.
+  for(int nidx = 0; nidx < nodes_fes->GetNDofs(); ++nidx)
+  {
+    Point2D q;
+    query_mesh->GetNode(nidx, q.data());
+
+    double wn {};
+    for(const auto& c : curves)
+    {
+      wn += axom::primal::winding_number(q, c);
+    }
+
+    winding[nidx] = wn;
+    inout[nidx] = std::round(wn);
+  }
+
+  // Save the query mesh and fields to disk using a format that can be viewed in VisIt
+  mfem::VisItDataCollection windingDC(outputPrefix, query_mesh.get());
+  windingDC.RegisterField("winding", &winding);
+  windingDC.RegisterField("inout", &inout);
+  windingDC.Save();
+
+  SLIC_INFO(axom::fmt::format("Outputting generated mesh '{}' to '{}'",
+                              windingDC.GetCollectionName(),
+                              axom::utilities::filesystem::getCWD()));
+
+  return 0;
+}

src/axom/slic/core/LogStream.cpp

@@ -69,6 +69,7 @@ std::string LogStream::getFormatedMessage(const std::string& msgLevel,
                                           const std::string& message,
                                           const std::string& tagName,
                                           const std::string& rank,
+                                          const std::string& rank_count,
                                           const std::string& fileName,
                                           int line)
 {
@@ -79,6 +80,7 @@ std::string LogStream::getFormatedMessage(const std::string& msgLevel,
   this->replaceKey(msg, "<TAG>", tagName);
   this->replaceKey(msg, "<FILE>", fileName);
   this->replaceKey(msg, "<RANK>", rank);
+  this->replaceKey(msg, "<RANK_COUNT>", rank_count);
 
   if(line != MSG_IGNORE_LINE)
   {