From cb5c5c1ce92c1246a25dba5cbc32909701881058 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 12 Mar 2026 12:19:56 -0700 Subject: [PATCH 01/17] Adds option to step reader example to write out patch trimming curves as mfem file Uses new nurbs patch format, which was added to mfem after the 4.9 release. --- src/axom/quest/examples/quest_step_file.cpp | 178 ++++++++++++++++++++ 1 file changed, 178 insertions(+) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 2cd720d6c6..47c8265a93 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -16,6 +16,7 @@ #include "axom/quest.hpp" #include +#include #ifdef AXOM_USE_MPI #include @@ -343,6 +344,150 @@ class PatchParametricSpaceProcessor int m_numFillZeros {0}; }; +/** + * Class that writes a patch-wise MFEM NURBS mesh containing a patch's trimming curves + * + * Each trimming curve is output as a separate 1D NURBS "patch" embedded in 2D space (u,v). + * + * Note: Support for reading these meshes was added to mfem after mfem@4.9.0 + * and is not in the current release of VisIt (visit@3.4.2) + */ +class PatchMFEMTrimmingCurveWriter +{ +public: + PatchMFEMTrimmingCurveWriter() { } + + void setOutputDirectory(const std::string& dir) { m_outputDirectory = dir; } + void setVerbosity(bool verbosityFlag) { m_verbose = verbosityFlag; } + void setNumFillZeros(int num) + { + if(num >= 0) + { + m_numFillZeros = num; + } + } + + void writeMFEMForPatch(int patchIndex, const NURBSPatch& patch) const + { + const auto& curves = patch.getTrimmingCurves(); + const int numCurves = curves.size(); + if(numCurves == 0) + { + return; + } + + using axom::utilities::filesystem::joinPath; + + const std::string outDir = joinPath(m_outputDirectory, "mfem_trim_curves"); + if(!axom::utilities::filesystem::pathExists(outDir)) + { + axom::utilities::filesystem::makeDirsForPath(outDir); + } + + const std::string meshFilename = + joinPath(outDir, + axom::fmt::format("trim_curves_patch_{:0{}}.mesh", patchIndex, m_numFillZeros)); + + std::ofstream meshFile(meshFilename); + if(!meshFile.is_open()) + { + SLIC_WARNING(axom::fmt::format("Unable to open '{}' for writing.", meshFilename)); + return; + } + + axom::fmt::memory_buffer content; + + axom::fmt::format_to(std::back_inserter(content), "MFEM NURBS mesh v1.0\n\n"); + axom::fmt::format_to(std::back_inserter(content), + "# Trim curves for STEP NURBSPatch {}\n", + patchIndex); + axom::fmt::format_to(std::back_inserter(content), + "# Parametric bbox: [{:.17g}, {:.17g}] x [{:.17g}, {:.17g}]\n", + patch.getMinKnot_u(), + patch.getMaxKnot_u(), + patch.getMinKnot_v(), + patch.getMaxKnot_v()); + axom::fmt::format_to(std::back_inserter(content), "# Number of trimming curves: {}\n\n", numCurves); + + axom::fmt::format_to(std::back_inserter(content), "dimension\n1\n\n"); + + // One element per trimming curve; each element uses its own vertex pair. + axom::fmt::format_to(std::back_inserter(content), "elements\n{}\n", numCurves); + for(int i = 0; i < numCurves; ++i) + { + const int v0 = 2 * i; + const int v1 = 2 * i + 1; + axom::fmt::format_to(std::back_inserter(content), "1 1 {} {}\n", v0, v1); + } + + axom::fmt::format_to(std::back_inserter(content), "\nboundary\n0\n\n"); + + // Edge list provides the unique knotvector index and (optionally) encodes orientation. + axom::fmt::format_to(std::back_inserter(content), "edges\n{}\n", numCurves); + for(int i = 0; i < numCurves; ++i) + { + const int v0 = 2 * i; + const int v1 = 2 * i + 1; + axom::fmt::format_to(std::back_inserter(content), "{} {} {}\n", i, v0, v1); + } + + axom::fmt::format_to(std::back_inserter(content), "\nvertices\n{}\n\n", 2 * numCurves); + + axom::fmt::format_to(std::back_inserter(content), "patches\n\n"); + for(int i = 0; i < numCurves; ++i) + { + const auto& curve = curves[i]; + SLIC_ASSERT(curve.isValidNURBS()); + + const int degree = curve.getDegree(); + const int ncp = curve.getNumControlPoints(); + const auto& knots = curve.getKnots().getArray(); + SLIC_ASSERT(knots.size() == ncp + degree + 1); + + axom::fmt::format_to(std::back_inserter(content), + "# Curve {}: {} (degree {}, {} control points)\n", + i, + curve.isRational() ? "rational" : "polynomial", + degree, + ncp); + + axom::fmt::format_to(std::back_inserter(content), "knotvectors\n1\n"); + axom::fmt::format_to(std::back_inserter(content), "{} {}", degree, ncp); + for(const auto& kv : knots) + { + axom::fmt::format_to(std::back_inserter(content), " {:.17g}", kv); + } + axom::fmt::format_to(std::back_inserter(content), "\n\n"); + + axom::fmt::format_to(std::back_inserter(content), "dimension\n2\n\n"); + axom::fmt::format_to(std::back_inserter(content), "controlpoints_cartesian\n"); + + const auto& cps = curve.getControlPoints(); + const auto& wts = curve.getWeights(); + for(int j = 0; j < ncp; ++j) + { + const double w = curve.isRational() ? wts[j] : 1.0; + axom::fmt::format_to(std::back_inserter(content), + "{:.17g} {:.17g} {:.17g}\n", + cps[j][0], + cps[j][1], + w); + } + axom::fmt::format_to(std::back_inserter(content), "\n"); + } + + meshFile << axom::fmt::to_string(content); + meshFile.close(); + + SLIC_INFO_IF(m_verbose, axom::fmt::format("MFEM trim-curve mesh generated: '{}'", meshFilename)); + } + +private: + std::string m_outputDirectory; + bool m_verbose {false}; + int m_numFillZeros {0}; +}; + #ifdef AXOM_USE_MPI // utility function to help with MPI_Allreduce calls @@ -646,6 +791,11 @@ int main(int argc, char** argv) ->description("Generate SVG files for each NURBS patch?") ->capture_default_str(); + bool output_mfem_trim_curves {false}; + app.add_flag("--output-mfem-trim-curves", output_mfem_trim_curves) + ->description("Generate one MFEM NURBS mesh per trimmed patch containing its trimming curves") + ->capture_default_str(); + app.get_formatter()->column_width(50); try @@ -844,5 +994,33 @@ int main(int argc, char** argv) } } + //--------------------------------------------------------------------------- + // Optionally output an MFEM patch-wise NURBS mesh for each patch's trimming curves, only on root rank + //--------------------------------------------------------------------------- + if(output_mfem_trim_curves && is_root) + { + const std::string outDir = joinPath(output_dir, "mfem_trim_curves"); + if(!axom::utilities::filesystem::pathExists(outDir)) + { + axom::utilities::filesystem::makeDirsForPath(outDir); + } + + SLIC_INFO( + axom::fmt::format("Generating MFEM meshes for patch trimming curves in '{}' directory", outDir)); + + const int numPatches = patches.size(); + const int numFillZeros = static_cast(std::log10(numPatches)) + 1; + + PatchMFEMTrimmingCurveWriter writer; + writer.setVerbosity(verbosity); + writer.setOutputDirectory(output_dir); + writer.setNumFillZeros(numFillZeros); + + for(int index = 0; index < numPatches; ++index) + { + writer.writeMFEMForPatch(index, patches[index]); + } + } + return rc; } From 244f6dbfdcd9b64db5e327e040d8f2a5b540a6ad Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 12 Mar 2026 13:03:44 -0700 Subject: [PATCH 02/17] Update MFEMReader to read in NURBS patches (when supported) --- src/axom/quest/io/MFEMReader.cpp | 106 ++++++++++++++++++++++++++----- 1 file changed, 91 insertions(+), 15 deletions(-) diff --git a/src/axom/quest/io/MFEMReader.cpp b/src/axom/quest/io/MFEMReader.cpp index 2f02e590b9..69bbcb9fe4 100644 --- a/src/axom/quest/io/MFEMReader.cpp +++ b/src/axom/quest/io/MFEMReader.cpp @@ -19,6 +19,15 @@ #include #include +// MFEM does not support reading patch-based 1D NURBS meshes until after the v4.9 +// release. Prefer patch-based extraction only when the MFEM version is new +// enough, otherwise fall back to element-based extraction. +#ifndef MFEM_VERSION + #define MFEM_VERSION 0 +#endif + +#define AXOM_MFEM_MIN_VERSION_PATCH_BASED_1D_NURBS 40901 + namespace axom { namespace quest @@ -68,6 +77,27 @@ int read_mfem(const std::string &fileName, return MFEMReader::READ_FAILED; } + // lambda to extract the control points for element idx from the mfem mesh as an array of primal::Point + using ControlPoint = primal::Point; + auto get_controlpoints = [nodes, fes](int idx) -> axom::Array { + mfem::Array vdofs; + mfem::Vector v; + fes->GetElementVDofs(idx, vdofs); + nodes->GetSubVector(vdofs, v); + const auto ord = fes->GetOrdering(); + + const int ncp = v.Size() / 2; + axom::Array cp(0, ncp); + for(int i = 0; i < ncp; ++i) + { + ord == mfem::Ordering::byVDIM ? cp.push_back({v[i], v[i + ncp]}) + : cp.push_back({v[2 * i], v[2 * i + 1]}); + } + + return cp; + }; + +#if MFEM_VERSION >= AXOM_MFEM_MIN_VERSION_PATCH_BASED_1D_NURBS // lambda to extract the knot vector associated with curve idx. Converts from mfem::KnotVector to primal::KnotVector auto get_knots = [&mesh](int idx) -> primal::KnotVector { mfem::Array kvs; @@ -78,16 +108,21 @@ int read_mfem(const std::string &fileName, return primal::KnotVector(knots_view, kv.GetOrder()); }; - // lambda to extract the control points for curve idx from the mfem mesh as an array of primal::Point - using ControlPoint = primal::Point; - auto get_controlpoints = [nodes, fes](int idx) -> axom::Array { - mfem::Array vdofs; + // lambda to extract the weights for curve idx from the mfem mesh + // Patch-based NURBS meshes can have multiple elements (knot spans) per patch. + // For robust extraction, build control points/weights from patch DOFs (not element DOFs). + auto get_patch_controlpoints = [nodes, fes, &mesh](int patchId) -> axom::Array { + mfem::Array dofs; + mesh->NURBSext->GetPatchDofs(patchId, dofs); + + mfem::Array vdofs(dofs); + fes->DofsToVDofs(vdofs); + mfem::Vector v; - fes->GetElementVDofs(idx, vdofs); nodes->GetSubVector(vdofs, v); const auto ord = fes->GetOrdering(); - const int ncp = v.Size() / 2; + const int ncp = dofs.Size(); axom::Array cp(0, ncp); for(int i = 0; i < ncp; ++i) { @@ -98,20 +133,39 @@ int read_mfem(const std::string &fileName, return cp; }; - // lambda to extract the weights for curve idx from the mfem mesh - auto get_weights = [&mesh, fes](int idx) -> axom::Array { + auto get_patch_weights = [&mesh](int patchId) -> axom::Array { mfem::Array dofs; - fes->GetElementDofs(idx, dofs); + mesh->NURBSext->GetPatchDofs(patchId, dofs); - const int NW = dofs.Size(); - axom::Array w(NW, NW); + const int nw = dofs.Size(); + axom::Array w(nw, nw); - // wrap our array's buffer w/ an mfem::Vector for GetSubVector - mfem::Vector mfem_vec_weights(w.data(), NW); + mfem::Vector mfem_vec_weights(w.data(), nw); mesh->NURBSext->GetWeights().GetSubVector(dofs, mfem_vec_weights); return w; }; +#endif + +#if MFEM_VERSION < AXOM_MFEM_MIN_VERSION_PATCH_BASED_1D_NURBS + auto get_element_weights = [fes, &mesh](int elemId) -> axom::Array { + mfem::Array dofs; + fes->GetElementDofs(elemId, dofs); + + const int nw = dofs.Size(); + axom::Array w(nw, nw); + + mfem::Vector mfem_vec_weights(w.data(), nw); + mesh->NURBSext->GetWeights().GetSubVector(dofs, mfem_vec_weights); + + return w; + }; + + auto get_element_degree = [fes, &mesh](int elemId) -> int { + const mfem::Array &orders = mesh->NURBSext->GetOrders(); + return (elemId < orders.Size()) ? orders[elemId] : fes->GetOrder(elemId); + }; +#endif // lambda to check if the weights correspond to a rational curve. If they are all equal it is not rational auto is_rational = [](const axom::Array &weights) -> bool { @@ -138,18 +192,40 @@ int read_mfem(const std::string &fileName, const bool isNURBS = dynamic_cast(fec) != nullptr; if(isNURBS) { +#if MFEM_VERSION >= AXOM_MFEM_MIN_VERSION_PATCH_BASED_1D_NURBS + // When MFEM can read patch-based 1D NURBS meshes, prefer reading curves from + // patches (which can contain multiple knot spans). This patch-based + // extraction is also compatible with older MFEM NURBS mesh v1.0 files. const int num_patches = fes->GetNURBSext()->GetNP(); for(int patchId = 0; patchId < num_patches; ++patchId) { const int attribute = mesh->GetPatchAttribute(patchId); const auto kv = get_knots(patchId); - const auto cp = get_controlpoints(patchId); - const auto w = get_weights(patchId); + const auto cp = get_patch_controlpoints(patchId); + const auto w = get_patch_weights(patchId); is_rational(w) ? curvemap[attribute].push_back({cp, w, kv}) : curvemap[attribute].push_back({cp, kv}); } +#else + { + // MFEM versions prior to AXOM_MFEM_MIN_VERSION_PATCH_BASED_1D_NURBS do not + // support reading patch-based 1D NURBS meshes. In that case, treat each + // MFEM element as a single (rational) Bezier span. + for(int zoneId = 0; zoneId < mesh->GetNE(); ++zoneId) + { + const int attribute = mesh->GetAttribute(zoneId); + const int degree = get_element_degree(zoneId); + + const auto cp = get_controlpoints(zoneId); + const auto w = get_element_weights(zoneId); + + is_rational(w) ? curvemap[attribute].push_back({cp, w, degree}) + : curvemap[attribute].push_back({cp, degree}); + } + } +#endif } else { From 0627a45f254920028ff0d6f926901b8dd509c138 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Thu, 12 Mar 2026 15:10:09 -0700 Subject: [PATCH 03/17] Extract the wire index of each trimming curve per patch and use as mfem attribute --- src/axom/quest/examples/quest_step_file.cpp | 32 +++++++++--- src/axom/quest/io/PSTEPReader.cpp | 16 ++++++ src/axom/quest/io/PSTEPReader.hpp | 11 +++-- src/axom/quest/io/STEPReader.cpp | 54 ++++++++++++++++----- src/axom/quest/io/STEPReader.hpp | 19 ++++++++ 5 files changed, 107 insertions(+), 25 deletions(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 47c8265a93..427dc4f5f5 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -367,7 +367,9 @@ class PatchMFEMTrimmingCurveWriter } } - void writeMFEMForPatch(int patchIndex, const NURBSPatch& patch) const + void writeMFEMForPatch(int patchId, + const NURBSPatch& patch, + axom::ArrayView trimmingCurveWireIds) const { const auto& curves = patch.getTrimmingCurves(); const int numCurves = curves.size(); @@ -376,6 +378,14 @@ class PatchMFEMTrimmingCurveWriter return; } + SLIC_WARNING_IF(trimmingCurveWireIds.size() != numCurves, + axom::fmt::format( + "Trimming curve wire id list size mismatch for patch {}: ids={}, curves={}. " + "Falling back to a single wire id.", + patchId, + trimmingCurveWireIds.size(), + numCurves)); + using axom::utilities::filesystem::joinPath; const std::string outDir = joinPath(m_outputDirectory, "mfem_trim_curves"); @@ -385,8 +395,7 @@ class PatchMFEMTrimmingCurveWriter } const std::string meshFilename = - joinPath(outDir, - axom::fmt::format("trim_curves_patch_{:0{}}.mesh", patchIndex, m_numFillZeros)); + joinPath(outDir, axom::fmt::format("trim_curves_patch_{:0{}}.mesh", patchId, m_numFillZeros)); std::ofstream meshFile(meshFilename); if(!meshFile.is_open()) @@ -400,7 +409,7 @@ class PatchMFEMTrimmingCurveWriter axom::fmt::format_to(std::back_inserter(content), "MFEM NURBS mesh v1.0\n\n"); axom::fmt::format_to(std::back_inserter(content), "# Trim curves for STEP NURBSPatch {}\n", - patchIndex); + patchId); axom::fmt::format_to(std::back_inserter(content), "# Parametric bbox: [{:.17g}, {:.17g}] x [{:.17g}, {:.17g}]\n", patch.getMinKnot_u(), @@ -415,9 +424,12 @@ class PatchMFEMTrimmingCurveWriter axom::fmt::format_to(std::back_inserter(content), "elements\n{}\n", numCurves); for(int i = 0; i < numCurves; ++i) { + // MFEM attributes are 1-based; we store the STEP 0-based wire index + 1. + const int wireId = (trimmingCurveWireIds.size() == numCurves) ? trimmingCurveWireIds[i] : 0; + const int mfem_attribute = wireId + 1; const int v0 = 2 * i; const int v1 = 2 * i + 1; - axom::fmt::format_to(std::back_inserter(content), "1 1 {} {}\n", v0, v1); + axom::fmt::format_to(std::back_inserter(content), "{} 1 {} {}\n", mfem_attribute, v0, v1); } axom::fmt::format_to(std::back_inserter(content), "\nboundary\n0\n\n"); @@ -438,6 +450,8 @@ class PatchMFEMTrimmingCurveWriter { const auto& curve = curves[i]; SLIC_ASSERT(curve.isValidNURBS()); + const int wireId = (trimmingCurveWireIds.size() == numCurves) ? trimmingCurveWireIds[i] : 0; + const int mfem_attribute = wireId + 1; const int degree = curve.getDegree(); const int ncp = curve.getNumControlPoints(); @@ -445,7 +459,10 @@ class PatchMFEMTrimmingCurveWriter SLIC_ASSERT(knots.size() == ncp + degree + 1); axom::fmt::format_to(std::back_inserter(content), - "# Curve {}: {} (degree {}, {} control points)\n", + "# Curve wireId={} (mfem_attribute={}, output index {}): {} (degree {}, " + "{} control points)\n", + wireId, + mfem_attribute, i, curve.isRational() ? "rational" : "polynomial", degree, @@ -1018,7 +1035,8 @@ int main(int argc, char** argv) for(int index = 0; index < numPatches; ++index) { - writer.writeMFEMForPatch(index, patches[index]); + const int patch_id = stepReader.getPatchIds()[index]; + writer.writeMFEMForPatch(patch_id, patches[index], stepReader.getTrimmingCurveWireIds(index)); } } diff --git a/src/axom/quest/io/PSTEPReader.cpp b/src/axom/quest/io/PSTEPReader.cpp index 00513f7566..8a5051eebc 100644 --- a/src/axom/quest/io/PSTEPReader.cpp +++ b/src/axom/quest/io/PSTEPReader.cpp @@ -78,6 +78,13 @@ int PSTEPReader::read(bool validate_model) } } } + + // Broadcast stable ids that match the input STEP enumeration. + bcast_array(m_patchIds); + for(auto& wire_ids : m_trimmingCurveWireIds) + { + bcast_array(wire_ids); + } } break; //handle other ranks @@ -145,6 +152,15 @@ int PSTEPReader::read(bool validate_model) } } } + + // Receive stable ids that match the input STEP enumeration. + bcast_array(m_patchIds); + m_trimmingCurveWireIds.clear(); + m_trimmingCurveWireIds.resize(numPatches); + for(int i = 0; i < numPatches; ++i) + { + bcast_array(m_trimmingCurveWireIds[i]); + } } break; } diff --git a/src/axom/quest/io/PSTEPReader.hpp b/src/axom/quest/io/PSTEPReader.hpp index 7a4b62e907..cefb82d377 100644 --- a/src/axom/quest/io/PSTEPReader.hpp +++ b/src/axom/quest/io/PSTEPReader.hpp @@ -113,9 +113,10 @@ class PSTEPReader : public STEPReader constexpr int ARR_DIM = axom::detail::ArrayTraits::dimension; static_assert(ARR_DIM == 1 || ARR_DIM == 2); - // Check that the value_type of the array is either double or Point + // Check that the value_type of the array is either double, int, or Point using value_type = typename ArrayType::value_type; - static_assert(std::is_same_v || primal::detail::is_point_v); + static_assert(std::is_same_v || std::is_same_v || + primal::detail::is_point_v); const bool is_root = (m_my_rank == 0); @@ -149,9 +150,9 @@ class PSTEPReader : public STEPReader } // then, send/receive the data - if constexpr(std::is_same_v) + if constexpr(std::is_same_v || std::is_same_v) { - // handles Array and Array + // handles Array, Array, and Array bcast_data(arr.view()); } else if constexpr(primal::detail::is_point_v) @@ -190,4 +191,4 @@ class PSTEPReader : public STEPReader } // namespace quest } // namespace axom -#endif // QUEST_PSTEPREADER_HPP_ \ No newline at end of file +#endif // QUEST_PSTEPREADER_HPP_ diff --git a/src/axom/quest/io/STEPReader.cpp b/src/axom/quest/io/STEPReader.cpp index 3487687156..4c9249f55e 100644 --- a/src/axom/quest/io/STEPReader.cpp +++ b/src/axom/quest/io/STEPReader.cpp @@ -69,6 +69,7 @@ struct PatchData axom::primal::BoundingBox parametricBBox; axom::primal::BoundingBox physicalBBox; axom::Array trimmingCurves_originallyPeriodic; + axom::Array trimmingCurves_wireIds; }; using PatchDataMap = std::map; @@ -870,7 +871,6 @@ class StepFileProcessor for(TopExp_Explorer edgeExp(wire, TopAbs_EDGE); edgeExp.More(); edgeExp.Next(), ++edgeIndex) { const TopoDS_Edge& edge = TopoDS::Edge(edgeExp.Current()); - const int curveIndex = patch.getNumTrimmingCurves(); TopAbs_Orientation orientation = edge.Orientation(); const bool isReversed = (orientation == TopAbs_REVERSED); @@ -878,11 +878,10 @@ class StepFileProcessor if(m_verbose) { BRepAdaptor_Curve curveAdaptor(edge); - SLIC_INFO(axom::fmt::format("[Patch {} Wire {} Edge {} Curve {}] Curve type: '{}'", + SLIC_INFO(axom::fmt::format("[Patch {} Wire {} Edge {}] Curve type: '{}'", patchIndex, wireIndex, edgeIndex, - curveIndex, curveTypeMap[curveAdaptor.GetType()])); } @@ -900,6 +899,7 @@ class StepFileProcessor auto curve = curveProcessor.nurbsCurve(); patchData.trimmingCurves_originallyPeriodic.push_back( curveProcessor.curveWasOriginallyPeriodic()); + patchData.trimmingCurves_wireIds.push_back(wireIndex); if(isReversed) // Ensure consistency of curve w.r.t. patch { @@ -911,11 +911,10 @@ class StepFileProcessor patch.addTrimmingCurve(curve); SLIC_INFO_IF(m_verbose, - axom::fmt::format("[Patch {} Wire {} Edge {} Curve {}] Added curve: {}", + axom::fmt::format("[Patch {} Wire {} Edge {}] Added curve: {}", patchIndex, wireIndex, edgeIndex, - curveIndex, curve)); // Check to ensure that curve did not change geometrically after making non-periodic @@ -924,14 +923,12 @@ class StepFileProcessor opencascade::handle origCurve = Geom2dConvert::CurveToBSplineCurve(parametricCurve); const bool withinThreshold = curveProcessor.compareToCurve(origCurve, 25); - SLIC_WARNING_IF( - !withinThreshold, - axom::fmt::format("[Patch {} Wire {} Edge {} Curve {}] Trimming curve was not " - "within threshold after clamping.", - patchIndex, - wireIndex, - edgeIndex, - curveIndex)); + SLIC_WARNING_IF(!withinThreshold, + axom::fmt::format("[Patch {} Wire {} Edge {}] Trimming curve was not " + "within threshold after clamping.", + patchIndex, + wireIndex, + edgeIndex)); } // TODO: Check that curve control points are within UV patch after adjusting periodicity @@ -1280,6 +1277,16 @@ class PatchTriangulator std::string STEPReader::getFileUnits() const { return m_stepProcessor->getFileUnits(); } +axom::ArrayView STEPReader::getTrimmingCurveWireIds(int patchArrayIndex) const +{ + if(patchArrayIndex < 0 || patchArrayIndex >= static_cast(m_trimmingCurveWireIds.size())) + { + return {}; + } + + return m_trimmingCurveWireIds[patchArrayIndex].view(); +} + std::string STEPReader::getBRepStats() const { // early return if the step file has not been loaded @@ -1569,6 +1576,27 @@ int STEPReader::read(bool validate_model) m_stepProcessor->extractPatches(m_patches); m_stepProcessor->extractTrimmingCurves(m_patches); + // Record stable ids that match the input STEP enumeration, even if consumers + // later filter/skip patches or trimming curves. + m_patchIds.clear(); + m_patchIds.resize(m_patches.size()); + m_trimmingCurveWireIds.clear(); + m_trimmingCurveWireIds.resize(m_patches.size()); + + const auto& patchDataMap = m_stepProcessor->getPatchDataMap(); + for(int patchArrayIndex = 0; patchArrayIndex < m_patches.size(); ++patchArrayIndex) + { + // Current implementation preserves patch array index == input face index. + // Keep this explicit mapping in case future logic filters patches. + m_patchIds[patchArrayIndex] = patchArrayIndex; + + auto it = patchDataMap.find(patchArrayIndex); + if(it != patchDataMap.end()) + { + m_trimmingCurveWireIds[patchArrayIndex] = it->second.trimmingCurves_wireIds; + } + } + return 0; } diff --git a/src/axom/quest/io/STEPReader.hpp b/src/axom/quest/io/STEPReader.hpp index 4570a696b2..089d67cbec 100644 --- a/src/axom/quest/io/STEPReader.hpp +++ b/src/axom/quest/io/STEPReader.hpp @@ -43,6 +43,7 @@ class STEPReader using PatchArray = axom::Array; using NURBSCurve = axom::primal::NURBSCurve; + using IndexArray = axom::Array; STEPReader() = default; virtual ~STEPReader(); @@ -69,6 +70,22 @@ class STEPReader /// Get the number of patches in the read file int numPatches() { return m_patches.size(); } + /*! + * \brief Returns a 0-based patch id per entry in \a getPatchArray() + * + * This id corresponds to the face index in the input STEP file enumeration. + * If consumers skip patches, these ids will be non-contiguous. + */ + axom::ArrayView getPatchIds() const { return m_patchIds.view(); } + + /*! + * \brief Returns the 0-based wire index for each extracted trimming curve + * + * The i-th entry corresponds to `getPatchArray()[patchArrayIndex].getTrimmingCurves()[i]` + * and stores the 0-based wire index from the input face enumeration. + */ + axom::ArrayView getTrimmingCurveWireIds(int patchArrayIndex) const; + /// Returns some information about the loaded BRep std::string getBRepStats() const; @@ -105,6 +122,8 @@ class STEPReader bool m_verbosity {false}; internal::StepFileProcessor* m_stepProcessor {nullptr}; PatchArray m_patches; + IndexArray m_patchIds; + axom::Array m_trimmingCurveWireIds; }; } // namespace quest From f33b94e7955c00c88cce4a07628107f78d09a8f7 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 11:30:31 -0700 Subject: [PATCH 04/17] Adds primal::NURBSCurve::isLinear This is a NURBS analogue to a function in the BezierCurve class. --- src/axom/primal/geometry/NURBSCurve.hpp | 48 +++++++++++++++++++ src/axom/primal/tests/primal_nurbs_curve.cpp | 50 ++++++++++++++++++++ 2 files changed, 98 insertions(+) diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index 42f95b1eac..dcd5e7eb0b 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -680,6 +680,54 @@ class NURBSCurve return OrientedBoundingBoxType(m_controlPoints.data(), static_cast(m_controlPoints.size())); } + /*! + * \brief Predicate to check if the NURBS curve is approximately linear + * + * This function checks if the interior control points of the NURBSCurve + * are approximately on the line segment defined by its two endpoints. + * + * \param [in] tol Threshold for squared distance + * \param [in] useStrictLinear If true, checks that the control points are + * evenly spaced along the line and not too far from the line + * \return True if curve is near-linear + */ + bool isLinear(double tol = 1e-8, bool useStrictLinear = false) const + { + const int npts = getNumControlPoints(); + if(npts <= 2) + { + return true; + } + + const int end_idx = npts - 1; + + if(useStrictLinear) + { + for(int p = 1; p < end_idx; ++p) + { + const double t = p / static_cast(end_idx); + PointType the_pt = PointType::lerp(m_controlPoints[0], m_controlPoints[end_idx], t); + if(squared_distance(m_controlPoints[p], the_pt) > tol) + { + return false; + } + } + } + else + { + SegmentType seg(m_controlPoints[0], m_controlPoints[end_idx]); + for(int p = 1; p < end_idx; ++p) + { + if(squared_distance(m_controlPoints[p], seg) > tol) + { + return false; + } + } + } + + return true; + } + ///@} ///@{ diff --git a/src/axom/primal/tests/primal_nurbs_curve.cpp b/src/axom/primal/tests/primal_nurbs_curve.cpp index b7a12dea0e..52f1f7ade1 100644 --- a/src/axom/primal/tests/primal_nurbs_curve.cpp +++ b/src/axom/primal/tests/primal_nurbs_curve.cpp @@ -18,6 +18,56 @@ namespace primal = axom::primal; +//------------------------------------------------------------------------------ +TEST(primal_nurbscurve, is_linear_predicate) +{ + constexpr int DIM = 2; + using CoordType = double; + using PointType = primal::Point; + using NURBSCurveType = primal::NURBSCurve; + + constexpr double tol = 1e-12; + + // Degree-1 segment + { + NURBSCurveType c(2, 1); + c[0] = PointType {0.0, 0.0}; + c[1] = PointType {1.0, 0.0}; + EXPECT_TRUE(c.isLinear(tol)); + } + + // Degree-2, collinear control polygon + { + NURBSCurveType c(3, 2); + c[0] = PointType {0.0, 0.0}; + c[1] = PointType {0.5, 0.0}; + c[2] = PointType {1.0, 0.0}; + EXPECT_TRUE(c.isLinear(tol)); + } + + // Degree-2, non-collinear interior point + { + NURBSCurveType c(3, 2); + c[0] = PointType {0.0, 0.0}; + c[1] = PointType {0.5, 1e-3}; + c[2] = PointType {1.0, 0.0}; + EXPECT_FALSE(c.isLinear(tol)); + } + + // Rational, collinear should still be linear + { + NURBSCurveType c(3, 2); + c[0] = PointType {0.0, 0.0}; + c[1] = PointType {0.5, 0.0}; + c[2] = PointType {1.0, 0.0}; + c.makeRational(); + c.setWeight(0, 1.0); + c.setWeight(1, 2.0); + c.setWeight(2, 0.5); + EXPECT_TRUE(c.isLinear(tol)); + } +} + //------------------------------------------------------------------------------ TEST(primal_nurbscurve, default_constructor) { From e88a42e9eed829af8f4640cc3032277f007ec95f Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 12:05:05 -0700 Subject: [PATCH 05/17] Adds NURBSPatch::isTriviallyTrimmed predicate This returns true if the patch has not trimming curves, or has four trimming curves that form a loop around the patch boundary in parametric space. --- src/axom/primal/geometry/NURBSPatch.hpp | 236 ++++++++++++++++++- src/axom/primal/tests/primal_nurbs_patch.cpp | 212 +++++++++++++++++ 2 files changed, 447 insertions(+), 1 deletion(-) diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index eed10b824d..3c8848327c 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -2791,6 +2791,240 @@ class NURBSPatch /// \brief Get number of trimming curves int getNumTrimmingCurves() const { return m_trimmingCurves.size(); } + /*! + * \brief Predicate to check if the patch is "trivially trimmed" in parameter space. + * + * A patch is considered trivially trimmed if either: + * - it has no trimming curves, or + * - it has exactly four trimming curves and they form an axis-aligned rectangle + * on the patch's parametric boundary: + * - Two curves are horizontal (v=min_v and v=max_v) and have opposite directions + * - Two curves are vertical (u=min_u and u=max_u) and have opposite directions + * - Each curve is approximately linear in (u,v) space + * - Curve endpoints match the patch's (min_u,max_u,min_v,max_v) corner coordinates + * + * \param [in] tol Threshold for squared distance used by NURBSCurve::isLinear and + * for the axis-alignment check on the curve endpoints. + */ + bool isTriviallyTrimmed(double tol = 1e-8) const + { + const int ncurves = getNumTrimmingCurves(); + if(ncurves == 0) + { + return true; + } + + if(ncurves != 4) + { + return false; + } + + const double min_u = static_cast(getMinKnot_u()); + const double max_u = static_cast(getMaxKnot_u()); + const double min_v = static_cast(getMinKnot_v()); + const double max_v = static_cast(getMaxKnot_v()); + + enum Flag : unsigned int + { + // Edge placement on the patch boundary + HasUminVertical = 1u << 0, + HasUmaxVertical = 1u << 1, + HasVminHorizontal = 1u << 2, + HasVmaxHorizontal = 1u << 3, + + // Edge direction coverage (opposing directions required) + HasHorizPosU = 1u << 4, + HasHorizNegU = 1u << 5, + HasVertPosV = 1u << 6, + HasVertNegV = 1u << 7, + + // Corner coverage and parity (each corner must appear exactly twice across endpoints) + SeenCornerMinUMinV = 1u << 8, + SeenCornerMaxUMinV = 1u << 9, + SeenCornerMaxUMaxV = 1u << 10, + SeenCornerMinUMaxV = 1u << 11, + + ParityCornerMinUMinV = 1u << 12, + ParityCornerMaxUMinV = 1u << 13, + ParityCornerMaxUMaxV = 1u << 14, + ParityCornerMinUMaxV = 1u << 15, + + // Use these flags to define larger checks + HasAllEdges = HasUminVertical | HasUmaxVertical | HasVminHorizontal | HasVmaxHorizontal, + HasOppositeOrientations = HasHorizPosU | HasHorizNegU | HasVertPosV | HasVertNegV, + SeenAllCorners = + SeenCornerMinUMinV | SeenCornerMaxUMinV | SeenCornerMaxUMaxV | SeenCornerMinUMaxV, + ParityFlags = + ParityCornerMinUMinV | ParityCornerMaxUMinV | ParityCornerMaxUMaxV | ParityCornerMinUMaxV, + TriviallyTrimmedFlags = HasAllEdges | HasOppositeOrientations | SeenAllCorners + }; + + auto sq = [](double x) -> double { return x * x; }; + auto near = [&](double a, double b) -> bool { return sq(a - b) <= tol; }; + + enum BoundMask : unsigned int + { + UMin = 1u << 0, + UMax = 1u << 1, + VMin = 1u << 2, + VMax = 1u << 3, + + UMask = UMin | UMax, + VMask = VMin | VMax + }; + + auto bound_mask = [&](double u, double v) -> unsigned int { + unsigned int m = 0u; + if(near(u, min_u)) + { + m |= UMin; + } + if(near(u, max_u)) + { + m |= UMax; + } + if(near(v, min_v)) + { + m |= VMin; + } + if(near(v, max_v)) + { + m |= VMax; + } + return m; + }; + + auto toggle_corner = [](unsigned int& flags, unsigned int seen_bit, unsigned int parity_bit) { + flags |= seen_bit; + flags ^= parity_bit; + }; + + auto toggle_corner_for_mask = [&](unsigned int& flags, unsigned int mask) -> bool { + switch(mask) + { + case(UMin | VMin): + toggle_corner(flags, SeenCornerMinUMinV, ParityCornerMinUMinV); + return true; + case(UMax | VMin): + toggle_corner(flags, SeenCornerMaxUMinV, ParityCornerMaxUMinV); + return true; + case(UMax | VMax): + toggle_corner(flags, SeenCornerMaxUMaxV, ParityCornerMaxUMaxV); + return true; + case(UMin | VMax): + toggle_corner(flags, SeenCornerMinUMaxV, ParityCornerMinUMaxV); + return true; + default: + return false; + } + }; + + unsigned int flags = 0u; + + for(int i = 0; i < ncurves; ++i) + { + const auto& c = m_trimmingCurves[i]; + if(!c.isValidNURBS()) + { + return false; + } + + const auto& p0 = c.getInitPoint(); + const auto& p1 = c.getEndPoint(); + const double u0 = static_cast(p0[0]); + const double v0 = static_cast(p0[1]); + const double u1 = static_cast(p1[0]); + const double v1 = static_cast(p1[1]); + + const unsigned int m0 = bound_mask(u0, v0); + const unsigned int m1 = bound_mask(u1, v1); + + // Endpoints must match a corner of the patch boundary (with tolerance) + if(!toggle_corner_for_mask(flags, m0) || !toggle_corner_for_mask(flags, m1)) + { + return false; + } + + const unsigned int u0m = (m0 & UMask); + const unsigned int v0m = (m0 & VMask); + const unsigned int u1m = (m1 & UMask); + const unsigned int v1m = (m1 & VMask); + + const double du = u1 - u0; + const double dv = v1 - v0; + const double du2 = du * du; + const double dv2 = dv * dv; + + const bool is_horizontal = (dv2 <= tol) && (du2 > tol); + const bool is_vertical = (du2 <= tol) && (dv2 > tol); + if(!(is_horizontal || is_vertical)) + { + return false; + } + + if(is_horizontal) + { + // Must lie on v=min_v or v=max_v and span u=min_u..max_u + if(v0m != v1m) + { + return false; + } + if(!(v0m == VMin || v0m == VMax)) + { + return false; + } + + const unsigned int edge_bit = (v0m == VMin) ? HasVminHorizontal : HasVmaxHorizontal; + if((flags & edge_bit) != 0u) + { + return false; + } + flags |= edge_bit; + + if(u0m == u1m) + { + return false; + } + + flags |= (du > 0.0) ? HasHorizPosU : HasHorizNegU; + } + else // is_vertical + { + // Must lie on u=min_u or u=max_u and span v=min_v..max_v + if(u0m != u1m) + { + return false; + } + if(!(u0m == UMin || u0m == UMax)) + { + return false; + } + + const unsigned int edge_bit = (u0m == UMin) ? HasUminVertical : HasUmaxVertical; + if((flags & edge_bit) != 0u) + { + return false; + } + flags |= edge_bit; + + if(v0m == v1m) + { + return false; + } + + flags |= (dv > 0.0) ? HasVertPosV : HasVertNegV; + } + + // More expensive geometric check last. + if(!c.isLinear(tol)) + { + return false; + } + } + + return flags == TriviallyTrimmedFlags; + } + /// \brief use boolean flag for trimmed-ness bool isTrimmed() const { return m_isTrimmed; } @@ -4192,4 +4426,4 @@ template struct axom::fmt::formatter> : ostream_formatter { }; -#endif // AXOM_PRIMAL_NURBSPATCH_HPP_ \ No newline at end of file +#endif // AXOM_PRIMAL_NURBSPATCH_HPP_ diff --git a/src/axom/primal/tests/primal_nurbs_patch.cpp b/src/axom/primal/tests/primal_nurbs_patch.cpp index c4bd00b3e8..8c452faa9c 100644 --- a/src/axom/primal/tests/primal_nurbs_patch.cpp +++ b/src/axom/primal/tests/primal_nurbs_patch.cpp @@ -1196,6 +1196,218 @@ TEST(primal_nurbspatch, nurbs_parameter_space_scaling) } } +//------------------------------------------------------------------------------ +TEST(primal_nurbspatch, is_trivially_trimmed_predicate) +{ + constexpr int DIM = 3; + using CoordType = double; + using PointType = primal::Point; + using NURBSPatchType = primal::NURBSPatch; + using TrimmingCurveType = primal::NURBSCurve; + + constexpr double tol = 1e-12; + + // Simple bilinear patch geometry. + PointType controlPoints[2 * 2] = {PointType {0.0, 0.0, 0.0}, + PointType {0.0, 1.0, 0.0}, + PointType {1.0, 0.0, 0.0}, + PointType {1.0, 1.0, 0.0}}; + NURBSPatchType patch(controlPoints, 2, 2, 1, 1); + + // No trimming curves -> trivially trimmed + EXPECT_TRUE(patch.isTriviallyTrimmed(tol)); + + // Wrong number of curves -> not trivially trimmed + { + NURBSPatchType p = patch; + TrimmingCurveType c(2, 1); + c[0] = primal::Point {0.0, 0.0}; + c[1] = primal::Point {1.0, 0.0}; + p.addTrimmingCurve(c); + EXPECT_FALSE(p.isTriviallyTrimmed(tol)); + } + + // Add four axis-aligned linear trimming curves + { + NURBSPatchType p = patch; + + TrimmingCurveType c0(2, 1); + c0[0] = primal::Point {0.0, 0.0}; + c0[1] = primal::Point {1.0, 0.0}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {1.0, 0.0}; + c1[1] = primal::Point {1.0, 1.0}; + + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {1.0, 1.0}; + c2[1] = primal::Point {0.0, 1.0}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.0, 1.0}; + c3[1] = primal::Point {0.0, 0.0}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_TRUE(p.isTriviallyTrimmed(tol)); + } + + // Trivially-trimmed helper generates a trivially-trimmed patch + { + NURBSPatchType p = patch; + p.makeTriviallyTrimmed(); + EXPECT_TRUE(p.isTriviallyTrimmed(tol)); + } + + // Fuzzy boundary matching should succeed within tolerance + { + constexpr double loose_tol = 1e-10; // sqrt(loose_tol) ~ 1e-5 + constexpr double eps = 1e-6; + + NURBSPatchType p = patch; + + TrimmingCurveType c0(2, 1); + c0[0] = primal::Point {0.0 + eps, 0.0 + eps}; + c0[1] = primal::Point {1.0 - eps, 0.0 + eps}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {1.0 - eps, 0.0 + eps}; + c1[1] = primal::Point {1.0 - eps, 1.0 - eps}; + + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {1.0 - eps, 1.0 - eps}; + c2[1] = primal::Point {0.0 + eps, 1.0 - eps}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.0 + eps, 1.0 - eps}; + c3[1] = primal::Point {0.0 + eps, 0.0 + eps}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_TRUE(p.isTriviallyTrimmed(loose_tol)); + } + + // Four boundary curves, but horizontal directions match -> not trivially trimmed + { + NURBSPatchType p = patch; + + TrimmingCurveType c0(2, 1); + c0[0] = primal::Point {0.0, 0.0}; + c0[1] = primal::Point {1.0, 0.0}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {1.0, 0.0}; + c1[1] = primal::Point {1.0, 1.0}; + + // Top edge also left-to-right (same direction as bottom) + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {0.0, 1.0}; + c2[1] = primal::Point {1.0, 1.0}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.0, 1.0}; + c3[1] = primal::Point {0.0, 0.0}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_FALSE(p.isTriviallyTrimmed(tol)); + } + + // Four axis-aligned linear curves, but not aligned to patch boundaries -> not trivially trimmed + { + NURBSPatchType p = patch; + + TrimmingCurveType c0(2, 1); + c0[0] = primal::Point {0.1, 0.0}; + c0[1] = primal::Point {0.9, 0.0}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {0.9, 0.0}; + c1[1] = primal::Point {0.9, 1.0}; + + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {0.9, 1.0}; + c2[1] = primal::Point {0.1, 1.0}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.1, 1.0}; + c3[1] = primal::Point {0.1, 0.0}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_FALSE(p.isTriviallyTrimmed(tol)); + } + + // Four curves, but one is diagonal -> not trivially trimmed + { + NURBSPatchType p = patch; + + TrimmingCurveType c0(2, 1); + c0[0] = primal::Point {0.0, 0.0}; + c0[1] = primal::Point {1.0, 1.0}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {1.0, 1.0}; + c1[1] = primal::Point {0.0, 1.0}; + + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {0.0, 1.0}; + c2[1] = primal::Point {0.0, 0.0}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.0, 0.0}; + c3[1] = primal::Point {1.0, 0.0}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_FALSE(p.isTriviallyTrimmed(tol)); + } + + // Four curves, but one is non-linear -> not trivially trimmed + { + NURBSPatchType p = patch; + + TrimmingCurveType c0(3, 2); + c0[0] = primal::Point {0.0, 0.0}; + c0[1] = primal::Point {0.5, 1e-3}; + c0[2] = primal::Point {1.0, 0.0}; + + TrimmingCurveType c1(2, 1); + c1[0] = primal::Point {1.0, 0.0}; + c1[1] = primal::Point {1.0, 1.0}; + + TrimmingCurveType c2(2, 1); + c2[0] = primal::Point {1.0, 1.0}; + c2[1] = primal::Point {0.0, 1.0}; + + TrimmingCurveType c3(2, 1); + c3[0] = primal::Point {0.0, 1.0}; + c3[1] = primal::Point {0.0, 0.0}; + + p.addTrimmingCurve(c0); + p.addTrimmingCurve(c1); + p.addTrimmingCurve(c2); + p.addTrimmingCurve(c3); + + EXPECT_FALSE(p.isTriviallyTrimmed(tol)); + } +} + //------------------------------------------------------------------------------ TEST(primal_nurbspatch, bezier_extraction) { From 060b0b15b0696f914c5164b6507282289d8c9a4c Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 13:24:18 -0700 Subject: [PATCH 06/17] Adds an option to output a stats file per patch in quest_step_file example This lists statistics like the curve orders and min/max parametric space bounds per patch. --- src/axom/quest/examples/quest_step_file.cpp | 200 +++++++++++++++++++- 1 file changed, 199 insertions(+), 1 deletion(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 427dc4f5f5..327ff06a8a 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -17,6 +17,8 @@ #include #include +#include +#include #ifdef AXOM_USE_MPI #include @@ -505,6 +507,144 @@ class PatchMFEMTrimmingCurveWriter int m_numFillZeros {0}; }; +/** + * Class that writes a JSON stats file per patch summarizing trimming curves. + * + * The output includes: + * - number of wires (unique STEP wire indices) + * - number of trimming curves + * - min/max curve orders (order = degree + 1) + * - patch (u,v) knot-domain bounds + */ +class PatchTrimmingCurveStatsWriter +{ +public: + PatchTrimmingCurveStatsWriter() { } + + void setOutputDirectory(const std::string& dir) { m_outputDirectory = dir; } + void setVerbosity(bool verbosityFlag) { m_verbose = verbosityFlag; } + void setNumFillZeros(int num) + { + if(num >= 0) + { + m_numFillZeros = num; + } + } + + void writeStatsForPatch(int patchId, + const NURBSPatch& patch, + axom::ArrayView trimmingCurveWireIds) const + { + const auto& curves = patch.getTrimmingCurves(); + const int numCurves = curves.size(); + + int minOrder = 0; + int maxOrder = 0; + std::map order_histogram; + for(int i = 0; i < numCurves; ++i) + { + const int order = curves[i].getDegree() + 1; + ++order_histogram[order]; + if(i == 0) + { + minOrder = order; + maxOrder = order; + } + else + { + minOrder = std::min(minOrder, order); + maxOrder = std::max(maxOrder, order); + } + } + + int numWires = 0; + if(numCurves == 0) + { + numWires = 0; + } + else if(trimmingCurveWireIds.size() == numCurves) + { + std::set uniqueWireIds; + for(int i = 0; i < numCurves; ++i) + { + uniqueWireIds.insert(trimmingCurveWireIds[i]); + } + numWires = static_cast(uniqueWireIds.size()); + } + else + { + // If we can't trust the per-curve wire ids, conservatively report a single wire. + numWires = 1; + } + + using axom::utilities::filesystem::joinPath; + + const std::string outDir = joinPath(m_outputDirectory, "trim_curve_stats"); + if(!axom::utilities::filesystem::pathExists(outDir)) + { + axom::utilities::filesystem::makeDirsForPath(outDir); + } + + const std::string statsFilename = + joinPath(outDir, + axom::fmt::format("trim_curves_patch_{:0{}}.stats.json", patchId, m_numFillZeros)); + + std::ofstream statsFile(statsFilename); + if(!statsFile.is_open()) + { + SLIC_WARNING(axom::fmt::format("Unable to open '{}' for writing.", statsFilename)); + return; + } + + axom::fmt::memory_buffer content; + axom::fmt::format_to(std::back_inserter(content), "{{\n"); + axom::fmt::format_to(std::back_inserter(content), " \"patch_id\": {},\n", patchId); + axom::fmt::format_to(std::back_inserter(content), " \"num_wires\": {},\n", numWires); + axom::fmt::format_to(std::back_inserter(content), " \"num_trimming_curves\": {},\n", numCurves); + axom::fmt::format_to(std::back_inserter(content), " \"min_curve_order\": {},\n", minOrder); + axom::fmt::format_to(std::back_inserter(content), " \"max_curve_order\": {},\n", maxOrder); + axom::fmt::format_to(std::back_inserter(content), " \"curves_by_order\": {{"); + { + bool first = true; + for(const auto& kv : order_histogram) + { + axom::fmt::format_to(std::back_inserter(content), + "{}\n \"{}\": {}", + first ? "" : ",", + kv.first, + kv.second); + first = false; + } + if(!order_histogram.empty()) + { + axom::fmt::format_to(std::back_inserter(content), "\n "); + } + } + axom::fmt::format_to(std::back_inserter(content), "}},\n"); + axom::fmt::format_to(std::back_inserter(content), + " \"uv_bbox\": {{\n" + " \"min\": [{:.17g}, {:.17g}],\n" + " \"max\": [{:.17g}, {:.17g}]\n" + " }}\n", + patch.getMinKnot_u(), + patch.getMinKnot_v(), + patch.getMaxKnot_u(), + patch.getMaxKnot_v()); + axom::fmt::format_to(std::back_inserter(content), "}}\n"); + + statsFile << axom::fmt::to_string(content); + statsFile.close(); + + SLIC_INFO_IF(m_verbose, + axom::fmt::format("Trim-curve stats JSON generated: '{}'", statsFilename)); + } + +private: + std::string m_outputDirectory; + bool m_verbose {false}; + int m_numFillZeros {0}; +}; + #ifdef AXOM_USE_MPI // utility function to help with MPI_Allreduce calls @@ -813,6 +953,19 @@ int main(int argc, char** argv) ->description("Generate one MFEM NURBS mesh per trimmed patch containing its trimming curves") ->capture_default_str(); + bool skip_trivial_trimmed_patches {false}; + app.add_flag("--skip-trivial-trimmed-patches", skip_trivial_trimmed_patches) + ->description( + "Skip patch-wise outputs for trivially-trimmed patches (4 axis-aligned linear boundary " + "curves). " + "Applies to SVG, MFEM trim-curve meshes and trim-curve stats JSON when enabled.") + ->capture_default_str(); + + bool output_trim_curve_stats_json {false}; + app.add_flag("--output-trim-curve-stats-json", output_trim_curve_stats_json) + ->description("Generate one JSON stats file per patch summarizing its trimming curves") + ->capture_default_str(); + app.get_formatter()->column_width(50); try @@ -1007,6 +1160,10 @@ int main(int argc, char** argv) for(int index = 0; index < numPatches; ++index) { + if(skip_trivial_trimmed_patches && patches[index].isTriviallyTrimmed()) + { + continue; + } patchProcessor.generateSVGForPatch(index, patches[index]); } } @@ -1036,7 +1193,48 @@ int main(int argc, char** argv) for(int index = 0; index < numPatches; ++index) { const int patch_id = stepReader.getPatchIds()[index]; - writer.writeMFEMForPatch(patch_id, patches[index], stepReader.getTrimmingCurveWireIds(index)); + const auto wire_ids = stepReader.getTrimmingCurveWireIds(index); + if(skip_trivial_trimmed_patches && patches[index].isTriviallyTrimmed()) + { + continue; + } + + writer.writeMFEMForPatch(patch_id, patches[index], wire_ids); + } + } + + //--------------------------------------------------------------------------- + // Optionally output trim-curve stats JSON per patch, only on root rank + //--------------------------------------------------------------------------- + if(output_trim_curve_stats_json && is_root) + { + const std::string outDir = joinPath(output_dir, "trim_curve_stats"); + if(!axom::utilities::filesystem::pathExists(outDir)) + { + axom::utilities::filesystem::makeDirsForPath(outDir); + } + + SLIC_INFO( + axom::fmt::format("Generating JSON trim-curve stats for patches in '{}' directory", outDir)); + + const int numPatches = patches.size(); + const int numFillZeros = static_cast(std::log10(numPatches)) + 1; + + PatchTrimmingCurveStatsWriter stats_writer; + stats_writer.setVerbosity(verbosity); + stats_writer.setOutputDirectory(output_dir); + stats_writer.setNumFillZeros(numFillZeros); + + for(int index = 0; index < numPatches; ++index) + { + const int patch_id = stepReader.getPatchIds()[index]; + const auto wire_ids = stepReader.getTrimmingCurveWireIds(index); + if(skip_trivial_trimmed_patches && patches[index].isTriviallyTrimmed()) + { + continue; + } + + stats_writer.writeStatsForPatch(patch_id, patches[index], wire_ids); } } From 97613f66421f8fcd5108cc1deae57020744438b0 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 13:38:16 -0700 Subject: [PATCH 07/17] Adds additional patch stats (isTriviallyTrimmed & num knot spans in u- and v-) --- src/axom/quest/examples/quest_step_file.cpp | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 327ff06a8a..20ffe736a4 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -537,6 +537,10 @@ class PatchTrimmingCurveStatsWriter { const auto& curves = patch.getTrimmingCurves(); const int numCurves = curves.size(); + const bool is_trivially_trimmed = patch.isTriviallyTrimmed(); + + const int num_knot_spans_u = static_cast(patch.getKnots_u().getNumKnotSpans()); + const int num_knot_spans_v = static_cast(patch.getKnots_v().getNumKnotSpans()); int minOrder = 0; int maxOrder = 0; @@ -599,6 +603,15 @@ class PatchTrimmingCurveStatsWriter axom::fmt::memory_buffer content; axom::fmt::format_to(std::back_inserter(content), "{{\n"); axom::fmt::format_to(std::back_inserter(content), " \"patch_id\": {},\n", patchId); + axom::fmt::format_to(std::back_inserter(content), + " \"is_trivially_trimmed\": {},\n", + is_trivially_trimmed ? "true" : "false"); + axom::fmt::format_to(std::back_inserter(content), + " \"num_knot_spans_u\": {},\n", + num_knot_spans_u); + axom::fmt::format_to(std::back_inserter(content), + " \"num_knot_spans_v\": {},\n", + num_knot_spans_v); axom::fmt::format_to(std::back_inserter(content), " \"num_wires\": {},\n", numWires); axom::fmt::format_to(std::back_inserter(content), " \"num_trimming_curves\": {},\n", numCurves); axom::fmt::format_to(std::back_inserter(content), " \"min_curve_order\": {},\n", minOrder); @@ -958,7 +971,7 @@ int main(int argc, char** argv) ->description( "Skip patch-wise outputs for trivially-trimmed patches (4 axis-aligned linear boundary " "curves). " - "Applies to SVG, MFEM trim-curve meshes and trim-curve stats JSON when enabled.") + "Applies to SVG and MFEM trim-curve meshes when enabled.") ->capture_default_str(); bool output_trim_curve_stats_json {false}; @@ -1229,11 +1242,6 @@ int main(int argc, char** argv) { const int patch_id = stepReader.getPatchIds()[index]; const auto wire_ids = stepReader.getTrimmingCurveWireIds(index); - if(skip_trivial_trimmed_patches && patches[index].isTriviallyTrimmed()) - { - continue; - } - stats_writer.writeStatsForPatch(patch_id, patches[index], wire_ids); } } From 2dda7c04b3254c0d7bf18fd932a9f3bc6ea918a0 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 13:48:21 -0700 Subject: [PATCH 08/17] Adds stats to indicate original periodicity of each patch When reading in the STEP file, we remove the periodicity, but users might want to know if the patch was originally periodic. --- src/axom/quest/examples/quest_step_file.cpp | 16 ++++++++++++++-- src/axom/quest/io/PSTEPReader.cpp | 12 ++++++++++-- src/axom/quest/io/STEPReader.cpp | 20 ++++++++++++++++++++ src/axom/quest/io/STEPReader.hpp | 8 ++++++++ 4 files changed, 52 insertions(+), 4 deletions(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 20ffe736a4..247e3036b0 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -533,7 +533,9 @@ class PatchTrimmingCurveStatsWriter void writeStatsForPatch(int patchId, const NURBSPatch& patch, - axom::ArrayView trimmingCurveWireIds) const + axom::ArrayView trimmingCurveWireIds, + bool was_originally_periodic_u, + bool was_originally_periodic_v) const { const auto& curves = patch.getTrimmingCurves(); const int numCurves = curves.size(); @@ -606,6 +608,12 @@ class PatchTrimmingCurveStatsWriter axom::fmt::format_to(std::back_inserter(content), " \"is_trivially_trimmed\": {},\n", is_trivially_trimmed ? "true" : "false"); + axom::fmt::format_to(std::back_inserter(content), + " \"was_originally_periodic_u\": {},\n", + was_originally_periodic_u ? "true" : "false"); + axom::fmt::format_to(std::back_inserter(content), + " \"was_originally_periodic_v\": {},\n", + was_originally_periodic_v ? "true" : "false"); axom::fmt::format_to(std::back_inserter(content), " \"num_knot_spans_u\": {},\n", num_knot_spans_u); @@ -1242,7 +1250,11 @@ int main(int argc, char** argv) { const int patch_id = stepReader.getPatchIds()[index]; const auto wire_ids = stepReader.getTrimmingCurveWireIds(index); - stats_writer.writeStatsForPatch(patch_id, patches[index], wire_ids); + stats_writer.writeStatsForPatch(patch_id, + patches[index], + wire_ids, + stepReader.patchWasOriginallyPeriodic_u(index), + stepReader.patchWasOriginallyPeriodic_v(index)); } } diff --git a/src/axom/quest/io/PSTEPReader.cpp b/src/axom/quest/io/PSTEPReader.cpp index 8a5051eebc..4c1907a535 100644 --- a/src/axom/quest/io/PSTEPReader.cpp +++ b/src/axom/quest/io/PSTEPReader.cpp @@ -79,12 +79,16 @@ int PSTEPReader::read(bool validate_model) } } - // Broadcast stable ids that match the input STEP enumeration. + // Broadcast stable ids that match the input STEP enumeration bcast_array(m_patchIds); for(auto& wire_ids : m_trimmingCurveWireIds) { bcast_array(wire_ids); } + + // Broadcast periodicity flags for each patch + bcast_array(m_patchOriginallyPeriodic_u); + bcast_array(m_patchOriginallyPeriodic_v); } break; //handle other ranks @@ -153,7 +157,7 @@ int PSTEPReader::read(bool validate_model) } } - // Receive stable ids that match the input STEP enumeration. + // Receive stable ids that match the input STEP enumeration bcast_array(m_patchIds); m_trimmingCurveWireIds.clear(); m_trimmingCurveWireIds.resize(numPatches); @@ -161,6 +165,10 @@ int PSTEPReader::read(bool validate_model) { bcast_array(m_trimmingCurveWireIds[i]); } + + // Receive periodicity flags for each patch + bcast_array(m_patchOriginallyPeriodic_u); + bcast_array(m_patchOriginallyPeriodic_v); } break; } diff --git a/src/axom/quest/io/STEPReader.cpp b/src/axom/quest/io/STEPReader.cpp index 4c9249f55e..4019db4499 100644 --- a/src/axom/quest/io/STEPReader.cpp +++ b/src/axom/quest/io/STEPReader.cpp @@ -1560,6 +1560,18 @@ STEPReader::~STEPReader() } } +bool STEPReader::patchWasOriginallyPeriodic_u(int patchArrayIndex) const +{ + SLIC_ASSERT(patchArrayIndex >= 0 && patchArrayIndex < m_patchOriginallyPeriodic_u.size()); + return m_patchOriginallyPeriodic_u[patchArrayIndex] != 0; +} + +bool STEPReader::patchWasOriginallyPeriodic_v(int patchArrayIndex) const +{ + SLIC_ASSERT(patchArrayIndex >= 0 && patchArrayIndex < m_patchOriginallyPeriodic_v.size()); + return m_patchOriginallyPeriodic_v[patchArrayIndex] != 0; +} + int STEPReader::read(bool validate_model) { m_stepProcessor = new internal::StepFileProcessor(m_fileName, m_verbosity); @@ -1582,6 +1594,12 @@ int STEPReader::read(bool validate_model) m_patchIds.resize(m_patches.size()); m_trimmingCurveWireIds.clear(); m_trimmingCurveWireIds.resize(m_patches.size()); + m_patchOriginallyPeriodic_u.clear(); + m_patchOriginallyPeriodic_u.resize(m_patches.size()); + m_patchOriginallyPeriodic_u.fill(0); + m_patchOriginallyPeriodic_v.clear(); + m_patchOriginallyPeriodic_v.resize(m_patches.size()); + m_patchOriginallyPeriodic_v.fill(0); const auto& patchDataMap = m_stepProcessor->getPatchDataMap(); for(int patchArrayIndex = 0; patchArrayIndex < m_patches.size(); ++patchArrayIndex) @@ -1594,6 +1612,8 @@ int STEPReader::read(bool validate_model) if(it != patchDataMap.end()) { m_trimmingCurveWireIds[patchArrayIndex] = it->second.trimmingCurves_wireIds; + m_patchOriginallyPeriodic_u[patchArrayIndex] = it->second.wasOriginallyPeriodic_u ? 1 : 0; + m_patchOriginallyPeriodic_v[patchArrayIndex] = it->second.wasOriginallyPeriodic_v ? 1 : 0; } } diff --git a/src/axom/quest/io/STEPReader.hpp b/src/axom/quest/io/STEPReader.hpp index 089d67cbec..f3d91f8400 100644 --- a/src/axom/quest/io/STEPReader.hpp +++ b/src/axom/quest/io/STEPReader.hpp @@ -86,6 +86,12 @@ class STEPReader */ axom::ArrayView getTrimmingCurveWireIds(int patchArrayIndex) const; + /// \brief Returns whether the input STEP surface was originally periodic in u for this patch + bool patchWasOriginallyPeriodic_u(int patchArrayIndex) const; + + /// \brief Returns whether the input STEP surface was originally periodic in v for this patch + bool patchWasOriginallyPeriodic_v(int patchArrayIndex) const; + /// Returns some information about the loaded BRep std::string getBRepStats() const; @@ -124,6 +130,8 @@ class STEPReader PatchArray m_patches; IndexArray m_patchIds; axom::Array m_trimmingCurveWireIds; + IndexArray m_patchOriginallyPeriodic_u; + IndexArray m_patchOriginallyPeriodic_v; }; } // namespace quest From 06f644b863779e6ece50bcd4fb16f7b987d896d9 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 14:07:53 -0700 Subject: [PATCH 09/17] Improves output command line options for quest_step_reader example --- src/axom/quest/examples/quest_step_file.cpp | 71 +++++++++++---------- 1 file changed, 39 insertions(+), 32 deletions(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 247e3036b0..653d485dbf 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -925,8 +925,11 @@ int main(int argc, char** argv) axom::fmt::format("Validate the model while reading it in? (default: {})", validate_model)) ->capture_default_str(); + // Output options ----------------------------------------------------------- + auto* output_opts = app.add_option_group("Output", "Parameters associated with output"); + std::string output_dir = "step_output"; - app.add_option("-o,--output-dir", output_dir) + output_opts->add_option("-o,--out,--output-dir", output_dir) ->description("Output directory for generated meshes") ->capture_default_str() ->check([](const std::string& dir) -> std::string { @@ -937,56 +940,57 @@ int main(int argc, char** argv) return std::string(); }); - TriangleMeshOutputType output_trimmed {TriangleMeshOutputType::VTK}; - app.add_option("--output-trimmed", output_trimmed) + bool output_svg {false}; + output_opts->add_flag("--svg,--output-svg", output_svg) + ->description("Generate SVG files for each NURBS patch") + ->capture_default_str(); + + bool output_mfem_trim_curves {false}; + output_opts->add_flag("--mfem-trim-curves,--output-mfem-trim-curves", output_mfem_trim_curves) + ->description("Generate one MFEM NURBS mesh per trimmed patch containing its trimming curves") + ->capture_default_str(); + + bool output_trim_curve_stats_json {false}; + output_opts->add_flag("--stats-json,--output-trim-curve-stats-json", output_trim_curve_stats_json) + ->description("Generate one JSON stats file per patch summarizing its trimming curves") + ->capture_default_str(); + + bool skip_trivial_trimmed_patches {false}; + output_opts + ->add_flag("--skip-trivial,--skip-trivial-trimmed-patches", skip_trivial_trimmed_patches) + ->description("Skip patch-wise SVG/MFEM outputs for trivially-trimmed patches") + ->capture_default_str(); + + // Triangulation options ---------------------------------------------------- + auto* tri_opts = app.add_option_group("Triangulation", "Parameters associated with triangulation"); + + TriangleMeshOutputType output_trimmed {TriangleMeshOutputType::NONE}; + tri_opts->add_option("--tri.trimmed,--output-trimmed", output_trimmed) ->description("Output format for trimmed model triangulation: 'none', 'vtk', 'stl'") ->capture_default_str() ->transform(axom::CLI::CheckedTransformer(validTriangleMeshOutputs)); TriangleMeshOutputType output_untrimmed {TriangleMeshOutputType::NONE}; - app.add_option("--output-untrimmed", output_untrimmed) + tri_opts->add_option("--tri.untrimmed,--output-untrimmed", output_untrimmed) ->description("Output format for untrimmed model triangulation: 'none', 'vtk', 'stl'") ->capture_default_str() ->transform(axom::CLI::CheckedTransformer(validTriangleMeshOutputs)); double deflection {.1}; - app.add_option("--deflection", deflection) + tri_opts->add_option("--deflection", deflection) ->description("Max distance between actual geometry and triangulated geometry") ->capture_default_str(); bool relative_deflection {false}; - app.add_flag("--relative", relative_deflection) + tri_opts->add_flag("--relative", relative_deflection) ->description("Use relative deflection instead of absolute?") ->capture_default_str(); double angular_deflection {0.5}; - app.add_option("--angular-deflection", angular_deflection) + tri_opts->add_option("--angular-deflection", angular_deflection) ->description("Angular deflection between adjacent normals when triangulating surfaces") ->capture_default_str(); - bool output_svg {false}; - app.add_flag("--output-svg", output_svg) - ->description("Generate SVG files for each NURBS patch?") - ->capture_default_str(); - - bool output_mfem_trim_curves {false}; - app.add_flag("--output-mfem-trim-curves", output_mfem_trim_curves) - ->description("Generate one MFEM NURBS mesh per trimmed patch containing its trimming curves") - ->capture_default_str(); - - bool skip_trivial_trimmed_patches {false}; - app.add_flag("--skip-trivial-trimmed-patches", skip_trivial_trimmed_patches) - ->description( - "Skip patch-wise outputs for trivially-trimmed patches (4 axis-aligned linear boundary " - "curves). " - "Applies to SVG and MFEM trim-curve meshes when enabled.") - ->capture_default_str(); - - bool output_trim_curve_stats_json {false}; - app.add_flag("--output-trim-curve-stats-json", output_trim_curve_stats_json) - ->description("Generate one JSON stats file per patch summarizing its trimming curves") - ->capture_default_str(); - app.get_formatter()->column_width(50); try @@ -1009,8 +1013,11 @@ int main(int argc, char** argv) #endif } - // Ensure output directory exists - if(is_root && !axom::utilities::filesystem::pathExists(output_dir)) + // Ensure output directory exists iff we will write anything. + const bool will_write_output = (output_trimmed != TriangleMeshOutputType::NONE) || + (output_untrimmed != TriangleMeshOutputType::NONE) || output_svg || output_mfem_trim_curves || + output_trim_curve_stats_json; + if(is_root && will_write_output && !axom::utilities::filesystem::pathExists(output_dir)) { axom::utilities::filesystem::makeDirsForPath(output_dir); } From 8b18d9480ec385d8c135895938cfb2cb42a16dd3 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 14:27:38 -0700 Subject: [PATCH 10/17] Updates svg2contours to optionally output mfem meshes using patches format This lets us avoid degree elevation of all curves to degree 3. --- src/tools/svg2contours/svg2contours.py | 191 +++++++++++++++++++++++-- 1 file changed, 177 insertions(+), 14 deletions(-) diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index 31ab1a8c0d..0dd598f4a7 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -12,6 +12,13 @@ description: Reads in an SVG document and outputs an MFEM NURBS mesh. Depends on the svgpathtools module + + notes: + - By default, all SVG curve segments are output as cubic NURBS (degree 3) for + compatibility with older MFEM/VisIt workflows. + - This script optionally supports MFEM's newer "patches" NURBS mesh format + for 1D NURBS segments embedded in 2D. Support for reading patch-based 1D + NURBS meshes was added to MFEM after MFEM 4.9.0. """ import sys @@ -127,6 +134,18 @@ def to_complex(v): return bpoints2bezier([to_complex(tf.dot(to_point(p))) for p in cubic.bpoints()]) +def transform_segment(seg, tf): + """Apply transformation `tf` to a Line/QuadraticBezier/CubicBezier segment.""" + + def to_point(p): + return np.array([[p.real], [p.imag], [1.0]]) + + def to_complex(v): + return v.item(0) + 1j * v.item(1) + + return bpoints2bezier([to_complex(tf.dot(to_point(p))) for p in seg.bpoints()]) + + def lerp(a, b, t): """linear interpolation from a to b with parameter t, typically between 0 and 1""" return (1 - t) * a + t * b @@ -240,6 +259,39 @@ def segment_as_cubic(seg, reverse_paths: bool): return (cubic, weights) +def segment_as_native_nurbs(seg, reverse_paths: bool): + """Convert an svgpathtools segment to a (Bezier) NURBS segment (degree 1/2/3) when possible. + + For elliptical arcs, this returns a rational cubic Bezier segment. + Returns (segment, weights, degree). + """ + + if isinstance(seg, Line): + out_seg = seg + weights = [1, 1] + degree = 1 + elif isinstance(seg, QuadraticBezier): + out_seg = seg + weights = [1, 1, 1] + degree = 2 + elif isinstance(seg, CubicBezier): + out_seg = seg + weights = [1, 1, 1, 1] + degree = 3 + elif isinstance(seg, Arc): + cubic, weights = arc_to_cubic(seg) + out_seg = cubic + degree = 3 + else: + raise Exception(f"'{type(seg)}' type not supported yet") + + if reverse_paths: + out_seg = out_seg.reversed() + weights.reverse() + + return (out_seg, weights, degree) + + def dist_to_ellipse(center, radius, angle, pt): cx, cy = center.real, center.imag rx, ry = radius.real, radius.imag @@ -330,6 +382,83 @@ def write_file(self, filename): f.write("\n".join(mfem_file)) +class MFEMPatchesData: + + def __init__(self): + self.elem_cnt = 0 + self.vert_cnt = 0 + self.elems = [] + self.edges = [] + self.patches = [] + + @staticmethod + def _bezier_knotvector(degree: int): + # Bezier knot vector on [0,1] + return [0] * (degree + 1) + [1] * (degree + 1) + + def add_bezier(self, seg, degree: int, weights, attrib: int): + v0 = self.vert_cnt + v1 = self.vert_cnt + 1 + self.elems.append(" ".join(map(str, [attrib, 1, v0, v1]))) + self.edges.append(f"{self.elem_cnt} {v0} {v1}") + self.vert_cnt += 2 + self.elem_cnt += 1 + + cps = seg.bpoints() + if len(cps) != degree + 1: + raise Exception( + f"Expected {degree + 1} control points for degree {degree}, got {len(cps)}") + if len(weights) != len(cps): + raise Exception(f"Expected {len(cps)} weights, got {len(weights)}") + + knots = self._bezier_knotvector(degree) + + patch_lines = [] + patch_lines.append("") + patch_lines.append( + f"# Patch {self.elem_cnt - 1}: degree {degree} ({len(cps)} control points)") + patch_lines.append("knotvectors") + patch_lines.append("1") + patch_lines.append("{} {} {}".format( + degree, + len(cps), + " ".join(str(k) for k in knots), + )) + patch_lines.append("") + patch_lines.append("dimension") + patch_lines.append("2") + patch_lines.append("") + patch_lines.append("controlpoints") + for (cp, w) in zip(cps, weights): + patch_lines.append(f"{cp.real} {cp.imag} {w}") + patch_lines.append("") + + self.patches.append("\n".join(patch_lines)) + + def write_file(self, filename): + mfem_file = [] + + mfem_file.extend([ + "MFEM NURBS mesh v1.0", + "", + "#", + "# Patch-based 1D NURBS segments embedded in 2D.", + "# NOTE: MFEM support for reading patch-based 1D NURBS meshes was added after MFEM 4.9.0.", + "#", + "", + ]) + + mfem_file.extend(["dimension", "1", ""]) + mfem_file.extend(["elements", f"{self.elem_cnt}", "\n".join(self.elems), ""]) + mfem_file.extend(["boundary", "0", ""]) + mfem_file.extend(["edges", f"{self.elem_cnt}", "\n".join(self.edges), ""]) + mfem_file.extend(["vertices", f"{self.vert_cnt}", ""]) + mfem_file.extend(["patches", "\n".join(self.patches)]) + + with open(filename, mode="w") as f: + f.write("\n".join(mfem_file)) + + def compute_svg_path_stats(paths): stats = { "paths_total": len(paths), @@ -378,6 +507,18 @@ def parse_args(): help="Output file in mfem NURBS mesh format (*.mesh)", ) + parser.add_argument( + "--mfem-patches", + dest="mfem_patches", + default=False, + action="store_true", + help= + ("Write the newer MFEM NURBS 'patches' mesh format for 1D segments embedded in 2D " + "(requires MFEM > 4.9.0 to read). When enabled, Lines/Quadratic/Cubic segments are " + "written using degree 1/2/3 respectively; elliptical arcs are written as rational cubics." + ), + ) + parser.add_argument( "--stats", dest="statsfile", @@ -468,6 +609,8 @@ def main(): print("SVG paths: \n", paths) mfem_data = MFEMData() + mfem_patches = opts.get("mfem_patches", False) + mfem_patches_data = MFEMPatchesData() if mfem_patches else None for p_idx, p in enumerate(paths): # print(f"""reading {p_idx=} {p=} \n w/ {p.d()=}""") @@ -491,23 +634,43 @@ def main(): # in `arc_to_cubic` algorithm arc1, arc2 = seg.split(0.5) - cubic, weights = segment_as_cubic(arc1, reverse_paths) - xformed_cubic = transform_cubic(cubic, coordinate_transform) - mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) - - cubic, weights = segment_as_cubic(arc2, reverse_paths) - xformed_cubic = transform_cubic(cubic, coordinate_transform) - mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) + if mfem_patches: + seg1, weights1, degree1 = segment_as_native_nurbs(arc1, reverse_paths) + xformed_seg1 = transform_segment(seg1, coordinate_transform) + mfem_patches_data.add_bezier(xformed_seg1, degree1, weights1, attrib) + + seg2, weights2, degree2 = segment_as_native_nurbs(arc2, reverse_paths) + xformed_seg2 = transform_segment(seg2, coordinate_transform) + mfem_patches_data.add_bezier(xformed_seg2, degree2, weights2, attrib) + else: + cubic, weights = segment_as_cubic(arc1, reverse_paths) + xformed_cubic = transform_cubic(cubic, coordinate_transform) + mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) + + cubic, weights = segment_as_cubic(arc2, reverse_paths) + xformed_cubic = transform_cubic(cubic, coordinate_transform) + mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) else: - cubic, weights = segment_as_cubic(seg, reverse_paths) - xformed_cubic = transform_cubic(cubic, coordinate_transform) - mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) + if mfem_patches: + seg0, weights0, degree0 = segment_as_native_nurbs(seg, reverse_paths) + xformed_seg0 = transform_segment(seg0, coordinate_transform) + mfem_patches_data.add_bezier(xformed_seg0, degree0, weights0, attrib) + else: + cubic, weights = segment_as_cubic(seg, reverse_paths) + xformed_cubic = transform_cubic(cubic, coordinate_transform) + mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) output_file = opts["outputfile"] - mfem_data.write_file(output_file) - print( - f"Wrote '{output_file}' with {mfem_data.vert_cnt} vertices and NURBS {mfem_data.elem_cnt} elements" - ) + if mfem_patches: + mfem_patches_data.write_file(output_file) + print( + f"Wrote '{output_file}' with {mfem_patches_data.vert_cnt} vertices and NURBS {mfem_patches_data.elem_cnt} elements (patches format)" + ) + else: + mfem_data.write_file(output_file) + print( + f"Wrote '{output_file}' with {mfem_data.vert_cnt} vertices and NURBS {mfem_data.elem_cnt} elements" + ) if __name__ == "__main__": From 8cdea67d78478033facb76239b4627928c6c9967 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sun, 15 Mar 2026 15:29:42 -0700 Subject: [PATCH 11/17] Output elliptic arcs as rational quadratics instead of rational cubics --- src/tools/svg2contours/README.md | 17 ++ src/tools/svg2contours/svg2contours.py | 299 ++++++++++++++++++++----- 2 files changed, 263 insertions(+), 53 deletions(-) diff --git a/src/tools/svg2contours/README.md b/src/tools/svg2contours/README.md index a24fa7a622..2647189a82 100644 --- a/src/tools/svg2contours/README.md +++ b/src/tools/svg2contours/README.md @@ -39,6 +39,23 @@ Wrote 'drawing.mesh' with 54 vertices and NURBS 27 elements ``` > :information_source: This assumes your Axom clone has the `data` submodule located at `/data` +### Optional: write MFEM patches-format NURBS meshes + +MFEM added support for reading patch-based 1D NURBS segments embedded in 2D after MFEM 4.9.0. +To write this newer format, pass `--mfem-patches`: + +```shell +> cd / +> uv run --project ../src/tools/svg2contours ../src/tools/svg2contours/svg2contours.py \ + -i ../data/contours/svg/shapes.svg --mfem-patches -o drawing_patches.mesh +``` + +In `--mfem-patches` mode, Lines/Quadratic/Cubic segments are written with degree 1/2/3 respectively, +and elliptical arcs are written as rational quadratics (degree 2). When an SVG arc is split into +multiple quadratic Bezier spans internally, the `--mfem-patches` output merges them into a single +multi-span quadratic NURBS patch (multi-knotvector) to reduce element/patch count. +Control points in the MFEM patches format are stored in homogeneous form `(x*w, y*w, w)`. + ### Run the quest winding number example Now that we have an MFEM NURBS mesh, we can run our winding number application diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index 0dd598f4a7..c07035cf82 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -24,6 +24,7 @@ import sys import os import json +from typing import Optional from svgpathtools import ( Document, Path, @@ -240,6 +241,119 @@ def area(p1, p2, p3): return (CubicBezier(q_0, c_1, c_2, q_1), [3, 1, 1, 3]) +def arc_to_quadratic_beziers(arc: Arc, *, max_sweep_angle_rad: float = np.pi / 2): + """Convert an svgpathtools Arc to rational QuadraticBezier segments. + + Notes: + - Rational quadratics represent conic sections exactly. + - We subdivide the arc into pieces (default: <= 90 degrees) to keep the + interior weight bounded away from zero. + - This function avoids relying on svgpathtools' internal angle units for + `arc.theta` / `arc.delta` (which may be degrees depending on version). + Instead, it reconstructs the start angle and sweep angle from the arc's + start/end points in the ellipse's local parameter space, then selects + the correct branch using `arc.point(0.5)` (and `arc.large_arc` as a hint). + Returns: + List[(QuadraticBezier, weights)] where weights has length 3 and the degree is 2. + """ + + if max_sweep_angle_rad <= 0: + raise ValueError("max_sweep_angle_rad must be positive") + + center = getattr(arc, "center", None) + radius = getattr(arc, "radius", None) + + if center is None or radius is None: + raise Exception("Arc is missing required attributes (center/radius)") + + rx = float(np.abs(radius.real)) + ry = float(np.abs(radius.imag)) + + if rx == 0.0 or ry == 0.0: + # Degenerate arc; let caller fall back (svgpathtools generally models these as Lines). + return [] + + # We build rational quadratic pieces by using the conic (tangent intersection) construction. + # This avoids any ambiguity about svgpathtools' internal angle units and tends to be robust. + + def cross(a: complex, b: complex) -> float: + return float(a.real * b.imag - a.imag * b.real) + + def dot(a: complex, b: complex) -> float: + return float(a.real * b.real + a.imag * b.imag) + + def try_quad_for_arc_piece(arc_piece: Arc): + p0 = arc_piece.start + p2 = arc_piece.end + m = arc_piece.point(0.5) + + d0 = arc_piece.derivative(0.0) + d2 = arc_piece.derivative(1.0) + + # Compute intersection of tangents at endpoints. + denom = cross(d0, d2) + if np.abs(denom) < 1e-14: + return None + + # p0 + s*d0 = p2 + t*d2 + s = cross((p2 - p0), d2) / denom + p1 = p0 + s * d0 + + # Solve for the middle weight w s.t. the rational quadratic matches the midpoint: + # M = (P0 + 2 w P1 + P2) / (2(1+w)) => 2 w (M - P1) = P0 + P2 - 2 M + u = m - p1 + uu = dot(u, u) + if uu <= 0.0: + return None + + rhs = p0 + p2 - 2.0 * m + w = dot(rhs, u) / (2.0 * uu) + + if not np.isfinite(w) or w <= 0.0: + return None + + # Consistency check (least-squares residual; should be near zero for true conics). + res = (p0 + p2 + 2.0 * w * p1) - (2.0 * (1.0 + w) * m) + scale2 = max(rx, ry) ** 2 + if dot(res, res) > 1e-8 * max(1.0, scale2): + return None + + return (QuadraticBezier(p0, p1, p2), [1.0, float(w), 1.0]) + + # Split until each piece's w is at least cos(max_sweep/2) (matches the circle-arc weight bound). + w_min = float(np.cos(0.5 * max_sweep_angle_rad)) + max_pieces = 1024 + + pieces = [] + work = [arc] + while work: + if len(pieces) + len(work) > max_pieces: + return [] + + a = work.pop(0) + quad = try_quad_for_arc_piece(a) + if quad is None: + a1, a2 = a.split(0.5) + work.insert(0, a2) + work.insert(0, a1) + continue + + _, wts = quad + if wts[1] < w_min - 1e-12: + a1, a2 = a.split(0.5) + work.insert(0, a2) + work.insert(0, a1) + continue + + pieces.append(quad) + + # Match legacy `arc_to_cubic` orientation handling: reverse when sweep flag is not set. + if not getattr(arc, "sweep", True): + pieces = [(q.reversed(), list(reversed(w))) for (q, w) in reversed(pieces)] + + return pieces + + def segment_as_cubic(seg, reverse_paths: bool): if isinstance(seg, Line): cubic, weights = line_to_cubic(seg) @@ -259,37 +373,38 @@ def segment_as_cubic(seg, reverse_paths: bool): return (cubic, weights) -def segment_as_native_nurbs(seg, reverse_paths: bool): - """Convert an svgpathtools segment to a (Bezier) NURBS segment (degree 1/2/3) when possible. +def segment_as_native_nurbs_segments(seg, reverse_paths: bool): + """Convert an svgpathtools segment to Bezier (NURBS) segments with native degree where possible. - For elliptical arcs, this returns a rational cubic Bezier segment. - Returns (segment, weights, degree). + Returns a list of (segment, weights, degree). + - Line: 1 segment, degree 1 + - QuadraticBezier: 1 segment, degree 2 + - CubicBezier: 1 segment, degree 3 + - Arc: 1+ rational QuadraticBezier segments (degree 2), subdivided for robustness """ + out = [] if isinstance(seg, Line): - out_seg = seg - weights = [1, 1] - degree = 1 + out = [(seg, [1, 1], 1)] elif isinstance(seg, QuadraticBezier): - out_seg = seg - weights = [1, 1, 1] - degree = 2 + out = [(seg, [1, 1, 1], 2)] elif isinstance(seg, CubicBezier): - out_seg = seg - weights = [1, 1, 1, 1] - degree = 3 + out = [(seg, [1, 1, 1, 1], 3)] elif isinstance(seg, Arc): - cubic, weights = arc_to_cubic(seg) - out_seg = cubic - degree = 3 + quads = arc_to_quadratic_beziers(seg) + if not quads: + # Fall back to a rational cubic if the arc conversion is ill-conditioned. + cubic, weights = arc_to_cubic(seg) + out = [(cubic, weights, 3)] + else: + out = [(q, w, 2) for (q, w) in quads] else: raise Exception(f"'{type(seg)}' type not supported yet") if reverse_paths: - out_seg = out_seg.reversed() - weights.reverse() + out = [(s.reversed(), list(reversed(w)), d) for (s, w, d) in reversed(out)] - return (out_seg, weights, degree) + return out def dist_to_ellipse(center, radius, angle, pt): @@ -396,27 +511,80 @@ def _bezier_knotvector(degree: int): # Bezier knot vector on [0,1] return [0] * (degree + 1) + [1] * (degree + 1) - def add_bezier(self, seg, degree: int, weights, attrib: int): + @staticmethod + def _quadratic_multispan_knotvector(num_spans: int, *, degree: int = 2): + if degree != 2: + raise ValueError("Only quadratic multispan knotvectors are supported here") + if num_spans < 1: + raise ValueError("num_spans must be >= 1") + + knots = [0.0] * (degree + 1) + if num_spans > 1: + # Use multiplicity=degree at internal knots so each span is a Bezier segment. + for i in range(1, num_spans): + t = float(i) / float(num_spans) + knots.extend([t] * degree) + knots.extend([1.0] * (degree + 1)) + return knots + + @staticmethod + def quadratic_beziers_to_multispan(quads_with_weights): + """Merge quadratic rational Bezier spans into a single quadratic multi-span NURBS patch. + + This uses internal knot multiplicity=degree (2), so each span remains a Bezier segment. + Returns (cps, weights, knots). + """ + + num_spans = len(quads_with_weights) + if num_spans < 1: + raise ValueError("Expected at least one span") + + cps = [] + weights = [] + for span_idx, (quad, wts) in enumerate(quads_with_weights): + bpts = quad.bpoints() + if len(bpts) != 3 or len(wts) != 3: + raise Exception( + "Expected quadratic Bezier spans with 3 control points and 3 weights") + if span_idx == 0: + cps.extend(bpts) + weights.extend(wts) + else: + cps.extend(bpts[1:]) + weights.extend(wts[1:]) + + knots = MFEMPatchesData._quadratic_multispan_knotvector(num_spans) + return cps, weights, knots + + def add_nurbs_patch(self, + *, + cps, + degree: int, + weights, + knots, + attrib: int, + patch_comment: Optional[str] = None): + if len(cps) != len(weights): + raise Exception(f"Expected {len(cps)} weights, got {len(weights)}") + expected_knots = 1 + degree + len(cps) + if len(knots) != expected_knots: + raise Exception(f"Expected {expected_knots} knots, got {len(knots)}") + v0 = self.vert_cnt v1 = self.vert_cnt + 1 + self.vert_cnt += 2 + self.elems.append(" ".join(map(str, [attrib, 1, v0, v1]))) self.edges.append(f"{self.elem_cnt} {v0} {v1}") - self.vert_cnt += 2 self.elem_cnt += 1 - cps = seg.bpoints() - if len(cps) != degree + 1: - raise Exception( - f"Expected {degree + 1} control points for degree {degree}, got {len(cps)}") - if len(weights) != len(cps): - raise Exception(f"Expected {len(cps)} weights, got {len(weights)}") - - knots = self._bezier_knotvector(degree) - patch_lines = [] patch_lines.append("") - patch_lines.append( - f"# Patch {self.elem_cnt - 1}: degree {degree} ({len(cps)} control points)") + if patch_comment: + patch_lines.append(f"# Patch {self.elem_cnt - 1}: {patch_comment}") + else: + patch_lines.append( + f"# Patch {self.elem_cnt - 1}: degree {degree} ({len(cps)} control points)") patch_lines.append("knotvectors") patch_lines.append("1") patch_lines.append("{} {} {}".format( @@ -430,11 +598,26 @@ def add_bezier(self, seg, degree: int, weights, attrib: int): patch_lines.append("") patch_lines.append("controlpoints") for (cp, w) in zip(cps, weights): - patch_lines.append(f"{cp.real} {cp.imag} {w}") + # MFEM expects NURBS control points in homogeneous form: (x*w, y*w, w). + ww = float(w) + patch_lines.append(f"{cp.real * ww} {cp.imag * ww} {ww}") patch_lines.append("") self.patches.append("\n".join(patch_lines)) + def add_bezier(self, seg, degree: int, weights, attrib: int): + cps = seg.bpoints() + if len(cps) != degree + 1: + raise Exception( + f"Expected {degree + 1} control points for degree {degree}, got {len(cps)}") + self.add_nurbs_patch( + cps=cps, + degree=degree, + weights=weights, + knots=self._bezier_knotvector(degree), + attrib=attrib, + ) + def write_file(self, filename): mfem_file = [] @@ -515,7 +698,7 @@ def parse_args(): help= ("Write the newer MFEM NURBS 'patches' mesh format for 1D segments embedded in 2D " "(requires MFEM > 4.9.0 to read). When enabled, Lines/Quadratic/Cubic segments are " - "written using degree 1/2/3 respectively; elliptical arcs are written as rational cubics." + "written using degree 1/2/3 respectively; elliptical arcs are written as rational quadratics." ), ) @@ -628,33 +811,43 @@ def main(): for seg_idx, seg in enumerate(p): # print(f"""processing {seg_idx=} {seg=}""") - if isinstance(seg, Arc) and seg.large_arc and is_d_path: + if (not mfem_patches) and isinstance(seg, Arc) and seg.large_arc and is_d_path: # split large elliptical arcs for easier processing # this simplifies the derivation of the internal control points # in `arc_to_cubic` algorithm arc1, arc2 = seg.split(0.5) - if mfem_patches: - seg1, weights1, degree1 = segment_as_native_nurbs(arc1, reverse_paths) - xformed_seg1 = transform_segment(seg1, coordinate_transform) - mfem_patches_data.add_bezier(xformed_seg1, degree1, weights1, attrib) - - seg2, weights2, degree2 = segment_as_native_nurbs(arc2, reverse_paths) - xformed_seg2 = transform_segment(seg2, coordinate_transform) - mfem_patches_data.add_bezier(xformed_seg2, degree2, weights2, attrib) - else: - cubic, weights = segment_as_cubic(arc1, reverse_paths) - xformed_cubic = transform_cubic(cubic, coordinate_transform) - mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) + cubic, weights = segment_as_cubic(arc1, reverse_paths) + xformed_cubic = transform_cubic(cubic, coordinate_transform) + mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) - cubic, weights = segment_as_cubic(arc2, reverse_paths) - xformed_cubic = transform_cubic(cubic, coordinate_transform) - mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) + cubic, weights = segment_as_cubic(arc2, reverse_paths) + xformed_cubic = transform_cubic(cubic, coordinate_transform) + mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) else: if mfem_patches: - seg0, weights0, degree0 = segment_as_native_nurbs(seg, reverse_paths) - xformed_seg0 = transform_segment(seg0, coordinate_transform) - mfem_patches_data.add_bezier(xformed_seg0, degree0, weights0, attrib) + out_segments = segment_as_native_nurbs_segments(seg, reverse_paths) + if isinstance(seg, Arc) and len(out_segments) > 1 and all( + d == 2 for (_, _, d) in out_segments): + quads_with_weights = [] + for out_seg, weights0, _ in out_segments: + xformed_seg0 = transform_segment(out_seg, coordinate_transform) + quads_with_weights.append((xformed_seg0, weights0)) + + cps, weights, knots = MFEMPatchesData.quadratic_beziers_to_multispan( + quads_with_weights) + mfem_patches_data.add_nurbs_patch( + cps=cps, + degree=2, + weights=weights, + knots=knots, + attrib=attrib, + patch_comment=f"quadratic (order 2, {len(out_segments)} spans)", + ) + else: + for out_seg, weights0, degree0 in out_segments: + xformed_seg0 = transform_segment(out_seg, coordinate_transform) + mfem_patches_data.add_bezier(xformed_seg0, degree0, weights0, attrib) else: cubic, weights = segment_as_cubic(seg, reverse_paths) xformed_cubic = transform_cubic(cubic, coordinate_transform) From c7aa8f7b672751df9a9348108e10209beb045d35 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Mon, 16 Mar 2026 13:11:42 -0700 Subject: [PATCH 12/17] Adds option to 2D winding number example to output degree `elevated` mfem mesh This is in support of downstream applications (such as VisIt) that do not yet support the variable-order mfem NURBS meshes (or patch-based mfem NURBS meshes) which became available after the mfem@4.9 release. --- .../examples/quest_winding_number_2d.cpp | 234 +++++++++++++++++- 1 file changed, 233 insertions(+), 1 deletion(-) diff --git a/src/axom/quest/examples/quest_winding_number_2d.cpp b/src/axom/quest/examples/quest_winding_number_2d.cpp index f6e32b7767..886e7016c0 100644 --- a/src/axom/quest/examples/quest_winding_number_2d.cpp +++ b/src/axom/quest/examples/quest_winding_number_2d.cpp @@ -5,7 +5,7 @@ // SPDX-License-Identifier: (BSD-3-Clause) /*! - * \file quest_winding_number.cpp + * \file quest_winding_number2d.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) @@ -25,6 +25,13 @@ #include "mfem.hpp" +#include +#include +#include +#include +#include +#include +#include #include namespace primal = axom::primal; @@ -35,6 +42,216 @@ using BoundingBox2D = primal::BoundingBox; using NURBSCurve2D = primal::NURBSCurve; +namespace +{ +/** + * This helper class takes in an mfem mesh (potentially with variable order curves from mfem>4.9) + * and writes out a version that is compatible with mfem@4.9 + * In particular, this allows us to visualize it with current versions of VisIt which do not yet support this feature. + * + * We can remove this class once downstream appications (such as VisIt) are updated to a version of mfem + * that support the NURBS patches format. + */ +class MFEM49ElevatedNURBSMeshWriter +{ +public: + explicit MFEM49ElevatedNURBSMeshWriter(double tol = 1e-12) : m_tol(tol) { } + + bool writeElevatedMesh(const std::string& input_file, const std::string& output_file) const + { + mfem::Mesh mesh(input_file, /*generate_edges=*/1, /*refine=*/1); + + if(mesh.NURBSext == nullptr) + { + SLIC_WARNING( + axom::fmt::format("Input mesh '{}' has no NURBS extension; skipping degree elevation", + input_file)); + return false; + } + + // Note: NURBS curve meshes in mfem@4.9 must all have the same degree, and their knotvectors must have the same length + const mfem::Array& orders = mesh.NURBSext->GetOrders(); + int max_order = 0; + for(int i = 0; i < orders.Size(); ++i) + { + max_order = std::max(max_order, orders[i]); + } + + if(max_order <= 0) + { + SLIC_WARNING( + axom::fmt::format("Input mesh '{}' has invalid NURBS orders; skipping degree elevation", + input_file)); + return false; + } + + mesh.DegreeElevate(max_order, max_order); + + makeKnotVectorsUniform(mesh); + + // Write the modified mesh to output_file + if(!writeMeshPrintFormat(mesh, output_file)) + { + SLIC_WARNING(axom::fmt::format("Failed to write elevated mesh '{}'", output_file)); + return false; + } + + SLIC_INFO(axom::fmt::format("Wrote elevated MFEM 4.9-compatible NURBS mesh '{}' (max order {})", + output_file, + max_order)); + return true; + } + +private: + struct KnotBin + { + std::int64_t key {0}; + double value {0.0}; + int max_multiplicity {0}; + }; + + // MFEM 4.9's 1D NURBS reader assumes all knotvectors have the same GetNE/GetNCP + // and uses KnotVec(0) when computing offsets. To generate MFEM 4.9/VisIt + // compatible files, insert knots so every knotvector matches the maximum + // knot multiplicity pattern present in the mesh. + void makeKnotVectorsUniform(mfem::Mesh& mesh) const + { + if(mesh.NURBSext == nullptr) + { + return; + } + if(mesh.Dimension() != 1 || mesh.NURBSext->Dimension() != 1) + { + return; + } + + const int nkv = mesh.NURBSext->GetNKV(); + if(nkv <= 1) + { + return; + } + + const double inv_tol = (m_tol > 0.0) ? (1.0 / m_tol) : 1.0e12; + auto key_for = [inv_tol](double t) -> std::int64_t { + return static_cast(std::llround(t * inv_tol)); + }; + + std::vector> counts(nkv); + std::unordered_map bins; + + for(int i = 0; i < nkv; ++i) + { + const mfem::KnotVector* kv = mesh.NURBSext->GetKnotVector(i); + if(kv == nullptr) + { + continue; + } + + auto& local = counts[i]; + for(int j = 0; j < kv->Size(); ++j) + { + const double value = (*kv)[j]; + const std::int64_t key = key_for(value); + + const int multiplicity = ++local[key]; + auto& bin = bins[key]; + bin.key = key; + bin.value = value; + bin.max_multiplicity = std::max(bin.max_multiplicity, multiplicity); + } + } + + std::vector sorted_bins; + sorted_bins.reserve(bins.size()); + for(const auto& it : bins) + { + sorted_bins.push_back(it.second); + } + std::sort(sorted_bins.begin(), sorted_bins.end(), [](const KnotBin& a, const KnotBin& b) { + if(a.value != b.value) + { + return a.value < b.value; + } + return a.key < b.key; + }); + + int total_to_insert = 0; + mfem::Array insertions(nkv); + for(int i = 0; i < nkv; ++i) + { + int need_total = 0; + for(const auto& bin : sorted_bins) + { + const auto found = counts[i].find(bin.key); + const int have = (found == counts[i].end()) ? 0 : found->second; + need_total += std::max(0, bin.max_multiplicity - have); + } + + insertions[i] = new mfem::Vector(need_total); + int pos = 0; + for(const auto& bin : sorted_bins) + { + const auto found = counts[i].find(bin.key); + const int have = (found == counts[i].end()) ? 0 : found->second; + const int need = std::max(0, bin.max_multiplicity - have); + for(int k = 0; k < need; ++k) + { + (*insertions[i])[pos++] = bin.value; + } + } + total_to_insert += need_total; + } + + if(total_to_insert > 0) + { + mesh.KnotInsert(insertions); + } + + for(int i = 0; i < nkv; ++i) + { + delete insertions[i]; + } + } + + bool writeMeshPrintFormat(const mfem::Mesh& mesh, const std::string& output_file) const + { + std::ostringstream oss; + oss.precision(16); + mesh.Print(oss); + + std::istringstream iss(oss.str()); + std::ofstream ofs(output_file); + if(!ofs.is_open()) + { + return false; + } + ofs.precision(16); + + // Note: mfem@4.9's mesh reader does not support the "NURBS2" finite element collection, + // if it occurs, we replace it with the (essentially equivalent) "NURBS" + std::string line; + while(std::getline(iss, line)) + { + constexpr const char* prefix = "FiniteElementCollection: NURBS"; + if(line.rfind(prefix, 0) == 0) + { + ofs << prefix << '\n'; + } + else + { + ofs << line << '\n'; + } + } + + return true; + } + +private: + double m_tol {1e-12}; +}; + +} // namespace + //------------------------------------------------------------------------------ // CLI input //------------------------------------------------------------------------------ @@ -50,6 +267,7 @@ class Input bool memoized {true}; bool vis {true}; bool stats {false}; + std::string elevatedMeshFile; const std::array valid_algorithms {"direct", "fast-approximation"}; std::string algorithm {valid_algorithms[1]}; // fast-approximation @@ -90,6 +308,12 @@ class Input app.add_flag("--stats,!--no-stats", stats, "Compute summary stats for query fields?") ->capture_default_str(); + app.add_option("--output-elevated-mesh", elevatedMeshFile) + ->description( + "Optional. Output MFEM mesh after elevating all NURBS curve orders to the maximum order of " + "input") + ->capture_default_str(); + // Options for query tolerances app.add_option("--edge-tol", tol.edge_tol) ->description("Relative edge tolerance for queries") @@ -219,6 +443,14 @@ int main(int argc, char** argv) axom::utilities::raii::AnnotationsWrapper annotation_raii_wrapper(input.annotationMode); AXOM_ANNOTATE_SCOPE("winding number example"); + if(!input.elevatedMeshFile.empty()) + { + AXOM_ANNOTATE_SCOPE("write_elevated_mesh"); + + MFEM49ElevatedNURBSMeshWriter writer; + writer.writeElevatedMesh(input.inputFile, input.elevatedMeshFile); + } + // Read curves from the MFEM mesh axom::Array curves; { From 327f7885b7834ec090b3bcec84e41cfb3ab97e8f Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Mon, 16 Mar 2026 17:26:19 -0700 Subject: [PATCH 13/17] Misc comments and cleanup --- src/axom/primal/geometry/NURBSCurve.hpp | 9 ++++----- src/tools/svg2contours/svg2contours.py | 7 +++---- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/axom/primal/geometry/NURBSCurve.hpp b/src/axom/primal/geometry/NURBSCurve.hpp index dcd5e7eb0b..0200324232 100644 --- a/src/axom/primal/geometry/NURBSCurve.hpp +++ b/src/axom/primal/geometry/NURBSCurve.hpp @@ -700,14 +700,13 @@ class NURBSCurve } const int end_idx = npts - 1; - if(useStrictLinear) { for(int p = 1; p < end_idx; ++p) { - const double t = p / static_cast(end_idx); - PointType the_pt = PointType::lerp(m_controlPoints[0], m_controlPoints[end_idx], t); - if(squared_distance(m_controlPoints[p], the_pt) > tol) + if(const double t = p / static_cast(end_idx); + squared_distance(m_controlPoints[p], + PointType::lerp(m_controlPoints[0], m_controlPoints[end_idx], t)) > tol) { return false; } @@ -715,7 +714,7 @@ class NURBSCurve } else { - SegmentType seg(m_controlPoints[0], m_controlPoints[end_idx]); + const SegmentType seg(m_controlPoints[0], m_controlPoints[end_idx]); for(int p = 1; p < end_idx; ++p) { if(squared_distance(m_controlPoints[p], seg) > tol) diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index c07035cf82..5b2f782d7a 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -248,9 +248,7 @@ def arc_to_quadratic_beziers(arc: Arc, *, max_sweep_angle_rad: float = np.pi / 2 - Rational quadratics represent conic sections exactly. - We subdivide the arc into pieces (default: <= 90 degrees) to keep the interior weight bounded away from zero. - - This function avoids relying on svgpathtools' internal angle units for - `arc.theta` / `arc.delta` (which may be degrees depending on version). - Instead, it reconstructs the start angle and sweep angle from the arc's + - This function reconstructs the start angle and sweep angle from the arc's start/end points in the ellipse's local parameter space, then selects the correct branch using `arc.point(0.5)` (and `arc.large_arc` as a hint). Returns: @@ -347,7 +345,7 @@ def try_quad_for_arc_piece(arc_piece: Arc): pieces.append(quad) - # Match legacy `arc_to_cubic` orientation handling: reverse when sweep flag is not set. + # Match `arc_to_cubic` orientation handling: reverse when sweep flag is not set. if not getattr(arc, "sweep", True): pieces = [(q.reversed(), list(reversed(w))) for (q, w) in reversed(pieces)] @@ -532,6 +530,7 @@ def quadratic_beziers_to_multispan(quads_with_weights): """Merge quadratic rational Bezier spans into a single quadratic multi-span NURBS patch. This uses internal knot multiplicity=degree (2), so each span remains a Bezier segment. + Note: It would be better to fix the knots/weights to keep it C1 Returns (cps, weights, knots). """ From ba6a109cd05906bdf3a07f038be7bb0aa0609bd0 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Mon, 16 Mar 2026 18:32:32 -0700 Subject: [PATCH 14/17] Fixes file filters for code checks * Only use checked in files (when using git) * Don't try to lint files from a .venv virtual environment --- src/cmake/AxomMacros.cmake | 55 +++++++++++++++++++++++++++++++------- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/src/cmake/AxomMacros.cmake b/src/cmake/AxomMacros.cmake index 192277fed0..0dcd45095f 100644 --- a/src/cmake/AxomMacros.cmake +++ b/src/cmake/AxomMacros.cmake @@ -34,20 +34,57 @@ macro(axom_add_code_checks) set(_base_dirs "axom" "examples" "thirdparty/tests" "tools") set(_ext_expressions "*.cpp" "*.hpp" "*.inl" "*.cxx" "*.hxx" "*.cc" "*.c" "*.h" "*.hh" - "*.F" "*.f" "*.f90" "*.F90" "*.py") + "*.F" "*.f" "*.f90" "*.F90" + "*.py" + "*.cmake" "CMakeLists.txt") - set(_glob_expressions) - foreach(_exp ${_ext_expressions}) - foreach(_base_dir ${_base_dirs}) - list(APPEND _glob_expressions "${PROJECT_SOURCE_DIR}/${_base_dir}/${_exp}") + set(_sources) + set(_use_git_sources FALSE) + + if(COMMAND blt_is_git_repo AND COMMAND blt_git) + blt_is_git_repo(OUTPUT_STATE _axom_is_git_repo + SOURCE_DIR ${PROJECT_SOURCE_DIR}) + if(_axom_is_git_repo) + set(_use_git_sources TRUE) + endif() + endif() + + if(_use_git_sources) + blt_git(SOURCE_DIR ${PROJECT_SOURCE_DIR} + GIT_COMMAND ls-files -- ${_base_dirs} + OUTPUT_VARIABLE _git_ls_files + RETURN_CODE _git_ls_files_result) + + if(_git_ls_files_result EQUAL 0 AND NOT _git_ls_files STREQUAL "") + string(REPLACE "\n" ";" _sources "${_git_ls_files}") + + # Keep only files handled by BLT language filters. + list(FILTER _sources INCLUDE REGEX "(\\.(cpp|hpp|inl|cxx|hxx|cc|c|h|hh|f|f90|py|cmake)|\\.F90|\\.F|CMakeLists\\.txt)$") + + # Make source paths absolute. + set(_sources_prefixed) + foreach(_tracked_file ${_sources}) + list(APPEND _sources_prefixed "${PROJECT_SOURCE_DIR}/${_tracked_file}") + endforeach() + set(_sources ${_sources_prefixed}) + else() + set(_use_git_sources FALSE) + endif() + else() + set(_glob_expressions) + foreach(_exp ${_ext_expressions}) + foreach(_base_dir ${_base_dirs}) + list(APPEND _glob_expressions "${PROJECT_SOURCE_DIR}/${_base_dir}/${_exp}") + endforeach() endforeach() - endforeach() - # Glob for list of files to run code checks on - set(_sources) - file(GLOB_RECURSE _sources ${_glob_expressions}) + # Glob for list of files to run code checks on + file(GLOB_RECURSE _sources ${_glob_expressions}) + endif() # Filter out exclusions + # Never run checks on local python environments + list(FILTER _sources EXCLUDE REGEX ".*[\\\\/]\\.venv[\\\\/].*") set(_exclude_expressions "${PROJECT_SOURCE_DIR}/axom/sidre/examples/lulesh2/*" "${PROJECT_SOURCE_DIR}/axom/slam/examples/lulesh2.0.3/*" From da15e3ea058c872300d1accb02b7542a548db497 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Tue, 17 Mar 2026 10:39:05 -0700 Subject: [PATCH 15/17] Guard against calling log(0) --- src/axom/quest/examples/quest_step_file.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/axom/quest/examples/quest_step_file.cpp b/src/axom/quest/examples/quest_step_file.cpp index 653d485dbf..8a84f3f05d 100644 --- a/src/axom/quest/examples/quest_step_file.cpp +++ b/src/axom/quest/examples/quest_step_file.cpp @@ -1178,7 +1178,7 @@ int main(int argc, char** argv) output_dir)); const int numPatches = patches.size(); - const int numFillZeros = static_cast(std::log10(numPatches)) + 1; + const int numFillZeros = (numPatches > 0) ? static_cast(std::log10(numPatches)) + 1 : 1; PatchParametricSpaceProcessor patchProcessor; patchProcessor.setUnits(stepReader.getFileUnits()); @@ -1211,7 +1211,7 @@ int main(int argc, char** argv) axom::fmt::format("Generating MFEM meshes for patch trimming curves in '{}' directory", outDir)); const int numPatches = patches.size(); - const int numFillZeros = static_cast(std::log10(numPatches)) + 1; + const int numFillZeros = (numPatches > 0) ? static_cast(std::log10(numPatches)) + 1 : 1; PatchMFEMTrimmingCurveWriter writer; writer.setVerbosity(verbosity); @@ -1246,7 +1246,7 @@ int main(int argc, char** argv) axom::fmt::format("Generating JSON trim-curve stats for patches in '{}' directory", outDir)); const int numPatches = patches.size(); - const int numFillZeros = static_cast(std::log10(numPatches)) + 1; + const int numFillZeros = (numPatches > 0) ? static_cast(std::log10(numPatches)) + 1 : 1; PatchTrimmingCurveStatsWriter stats_writer; stats_writer.setVerbosity(verbosity); From b7710461457395de83599f8ba055c8d9bd5a9e5e Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Tue, 17 Mar 2026 12:17:55 -0700 Subject: [PATCH 16/17] Updates RELEASE-NOTES --- RELEASE-NOTES.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 6135b5970a..c3c69498f0 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -23,6 +23,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Sidre: Added Conduit Node to the Python interface. - Adds yapf as a Python formatter. - Quest: Adds fast GWN methods for STL/Triangulated STEP input and linearized NURBS Curve input which leverage error-controlled approximation and a spatial index (BVH). +- Primal: Adds `NURBSCurve::isLinear()` to check if a curve is (nearly) flat (corresponding to `BezierCurve::isLinear()`) +- Primal: Adds `NURBSPatch::isTriviallyTrimmed()` to check if the trimming curves for a patch lie on the patch boundaries +- Quest: Adds support for reading mfem files with variable order NURBS curves (requires mfem>4.9). ### Changed - Primal: Axom's polygon clipping was modified to handle some corner cases. @@ -31,6 +34,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Many algorithms in Bump and Mir have been enhanced with `setAllocatorID()` methods to permit use of custom allocators. - Uberenv's spack updated to v1.1.1 - `radiuss-spack-configs` submodule removed. +- Updates CMake code check targets to only use checked in files (via `git ls-files`, when available) ### Fixed - Primal: Fixed `NURBSPatch::evaluateDerivatives` returning a matrix with empty values. The old behavior From 71a2dbae365b212670f94bd66fa60c696f859c7d Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Tue, 17 Mar 2026 16:01:59 -0700 Subject: [PATCH 17/17] MSVC appears to not like `near` as a variable name I think it is due to a `#define near` within the msvc headers. --- src/axom/primal/geometry/NURBSPatch.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/axom/primal/geometry/NURBSPatch.hpp b/src/axom/primal/geometry/NURBSPatch.hpp index 3c8848327c..48fd776f06 100644 --- a/src/axom/primal/geometry/NURBSPatch.hpp +++ b/src/axom/primal/geometry/NURBSPatch.hpp @@ -2860,7 +2860,7 @@ class NURBSPatch }; auto sq = [](double x) -> double { return x * x; }; - auto near = [&](double a, double b) -> bool { return sq(a - b) <= tol; }; + auto is_near = [&](double a, double b) -> bool { return sq(a - b) <= tol; }; enum BoundMask : unsigned int { @@ -2875,19 +2875,19 @@ class NURBSPatch auto bound_mask = [&](double u, double v) -> unsigned int { unsigned int m = 0u; - if(near(u, min_u)) + if(is_near(u, min_u)) { m |= UMin; } - if(near(u, max_u)) + if(is_near(u, max_u)) { m |= UMax; } - if(near(v, min_v)) + if(is_near(v, min_v)) { m |= VMin; } - if(near(v, max_v)) + if(is_near(v, max_v)) { m |= VMax; }