From 8d65bf5a8d1716bf78b8f9b004131ec9f0375a3d Mon Sep 17 00:00:00 2001 From: Chris White Date: Thu, 25 Jul 2024 14:40:46 -0700 Subject: [PATCH 01/27] split RANK_COUNT out of RANK for slic message formatting --- RELEASE-NOTES.md | 2 ++ src/axom/slic/core/LogStream.cpp | 2 ++ src/axom/slic/core/LogStream.hpp | 7 +++++- .../docs/sphinx/sections/architecture.rst | 15 ++++++++--- .../examples/basic/lumberjack_logging.cpp | 17 ++++++++++--- src/axom/slic/streams/GenericOutputStream.cpp | 1 + src/axom/slic/streams/LumberjackStream.cpp | 25 +++++++++---------- src/axom/slic/streams/SynchronizedStream.cpp | 1 + 8 files changed, 48 insertions(+), 22 deletions(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index e934f14f6e..e09e25104c 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -63,6 +63,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ external `fmt` and axom's vendored copy. - Turn off CMake finding dependencies on system paths. - `axom::Array`: trivially-copyable types with a non-trivial constructor are now initialized on the GPU. +- SLIC no longer outputs the rank count in the `RANK` format string in parallel loggers. You can access + the rank count via new format option `RANK_COUNT`. ### Removed - Removes config option `AXOM_ENABLE_ANNOTATIONS`. Annotations are now provided by `caliper` diff --git a/src/axom/slic/core/LogStream.cpp b/src/axom/slic/core/LogStream.cpp index 17e2efb5c0..42efc78bd9 100644 --- a/src/axom/slic/core/LogStream.cpp +++ b/src/axom/slic/core/LogStream.cpp @@ -69,6 +69,7 @@ std::string LogStream::getFormatedMessage(const std::string& msgLevel, const std::string& message, const std::string& tagName, const std::string& rank, + const std::string& rank_count, const std::string& fileName, int line) { @@ -79,6 +80,7 @@ std::string LogStream::getFormatedMessage(const std::string& msgLevel, this->replaceKey(msg, "", tagName); this->replaceKey(msg, "", fileName); this->replaceKey(msg, "", rank); + this->replaceKey(msg, "", rank_count); if(line != MSG_IGNORE_LINE) { diff --git a/src/axom/slic/core/LogStream.hpp b/src/axom/slic/core/LogStream.hpp index f69eb4218e..9f3b306381 100644 --- a/src/axom/slic/core/LogStream.hpp +++ b/src/axom/slic/core/LogStream.hpp @@ -58,7 +58,8 @@ class LogStream *
  • user-supplied tag
    • *
  • with the filename
    • *
  • with the line number
    • - *
  • with the MPI rank
    • + *
  • with the MPI rank(s)
    • + *
  • with the number of MPI ranks
    • *
  • date/time the message is logged
    • * * @@ -73,6 +74,7 @@ class LogStream * std::string( "* FILE=\n" ) + * std::string( "* LINE=\n" ) + * std::string( "* RANK=\n" ) + + * std::string( "* RANK_COUNT=\n" ) + * std::string( "***********************************\n" ); * \endcode */ @@ -151,6 +153,8 @@ class LogStream * \param [in] tagName user-supplied tag, may be MSG_IGNORE_TAG * \param [in] fileName filename where this message is logged, may be * MSG_IGNORE_FILE to ignore this field. + * \param [in] rank The MPI rank(s) that emitted this message + * \param [in] rank_count the number of MPI ranks that emitted this message * \param [in] line the line number within the file where the message is * logged. Likewise, may be set to MSG_IGNORE_LINE to ignore this field. * @@ -161,6 +165,7 @@ class LogStream const std::string& message, const std::string& tagName, const std::string& rank, + const std::string& rank_count, const std::string& fileName, int line); diff --git a/src/axom/slic/docs/sphinx/sections/architecture.rst b/src/axom/slic/docs/sphinx/sections/architecture.rst index 7b044f2be6..f2df63bc66 100644 --- a/src/axom/slic/docs/sphinx/sections/architecture.rst +++ b/src/axom/slic/docs/sphinx/sections/architecture.rst @@ -151,8 +151,14 @@ The list of keywords is summarized in the table below. | **** | A string tag associated with a given message, e.g., for| | | filtering during post-processing, etc. | +---------------------+--------------------------------------------------------+ -| **** | The MPI rank that emitted the message. Only applicable | -| | when the `Axom Toolkit`_ is compiled with MPI enabled | +| **** | The MPI rank(s) that emitted the message. | +| | Only applicable when Axom is compiled with MPI enabled | +| | and with MPI-aware :ref:`LogStream` instances, such as,| +| | the :ref:`SynchronizedStream` and | +| | :ref:`LumberjackStream`. | ++---------------------+--------------------------------------------------------+ +| **** | The number of MPI ranks that emitted the message. | +| | Only applicable when Axom is compiled with MPI enabled | | | and with MPI-aware :ref:`LogStream` instances, such as,| | | the :ref:`SynchronizedStream` and | | | :ref:`LumberjackStream`. | @@ -221,7 +227,7 @@ Generic Output Stream The :ref:`GenericOutputStream`, is a concrete implementation of the :ref:`LogStream` base class, that can be constructed by specifying: -#. A C++ ``std::ostream`` object instance, e.g., ``std::cout`, ``std::cerr`` for +#. A C++ ``std::ostream`` object instance, e.g., ``std::cout``, ``std::cerr`` for console output, or to a file by passing a C++ ``std::ofstream`` object, and, #. Optionally, a string that specifies the :ref:`logMessageFormat`. @@ -423,12 +429,13 @@ The ``MyStream`` class implements the ``LogStream::append()`` method of the int line, bool AXOM_UNUSED_PARAM(filtered_duplicates) ) { - assert( m_stream != nillptr ); + assert( m_stream != nullptr ); (*m_stream) << this->getFormatedMessage( message::getLevelAsString(msgLevel), message, tagName, "", + "", fileName, line ); } diff --git a/src/axom/slic/examples/basic/lumberjack_logging.cpp b/src/axom/slic/examples/basic/lumberjack_logging.cpp index 1d0b8f5a02..9bcb591a59 100644 --- a/src/axom/slic/examples/basic/lumberjack_logging.cpp +++ b/src/axom/slic/examples/basic/lumberjack_logging.cpp @@ -39,10 +39,19 @@ int main(int argc, char** argv) MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Initialize SLIC - std::string format = std::string("\n") + - std::string("\t\n") + std::string("\tLEVEL=\n") + - std::string("\tRANKS=\n") + std::string("\tFILE=\n") + - std::string("\tLINE=\n"); +constexpr const char* format = R"( + +\t +\tLEVEL= +\tRANKS= +\tRANK_COUNT= +\tFILE= +\tLINE= +)"; + // std::string format = std::string("\n") + + // std::string("\t\n") + std::string("\tLEVEL=\n") + + // std::string("\tRANKS=\n") + std::string("\tFILE=\n") + + // std::string("\tLINE=\n"); slic::initialize(); // Set SLIC logging level and Lumberjack Logging stream diff --git a/src/axom/slic/streams/GenericOutputStream.cpp b/src/axom/slic/streams/GenericOutputStream.cpp index 39de7df769..0922918657 100644 --- a/src/axom/slic/streams/GenericOutputStream.cpp +++ b/src/axom/slic/streams/GenericOutputStream.cpp @@ -74,6 +74,7 @@ void GenericOutputStream::append(message::Level msgLevel, message, tagName, "", + "", fileName, line); } diff --git a/src/axom/slic/streams/LumberjackStream.cpp b/src/axom/slic/streams/LumberjackStream.cpp index c9353cc057..42d24e3a9c 100644 --- a/src/axom/slic/streams/LumberjackStream.cpp +++ b/src/axom/slic/streams/LumberjackStream.cpp @@ -139,22 +139,21 @@ void LumberjackStream::write(bool local) if(m_lj->isOutputNode() || local) { - std::vector messages = m_lj->getMessages(); - - const int nmessages = static_cast(messages.size()); - std::string rankString; - for(int i = 0; i < nmessages; ++i) + for(const auto* curr_message : m_lj->getMessages()) { - rankString = std::to_string(messages[i]->count()) + ": " + - messages[i]->stringOfRanks(); + if (curr_message == nullptr) { + continue; + } + (*m_stream) << this->getFormatedMessage( message::getLevelAsString( - static_cast(messages[i]->level())), - messages[i]->text(), - messages[i]->tag(), - rankString, - messages[i]->fileName(), - messages[i]->lineNumber()); + static_cast(curr_message->level())), + curr_message->text(), + curr_message->tag(), + curr_message->stringOfRanks(), + std::to_string(curr_message->count()), + curr_message->fileName(), + curr_message->lineNumber()); } m_stream->flush(); diff --git a/src/axom/slic/streams/SynchronizedStream.cpp b/src/axom/slic/streams/SynchronizedStream.cpp index 12d00ed33c..0d42274d1d 100644 --- a/src/axom/slic/streams/SynchronizedStream.cpp +++ b/src/axom/slic/streams/SynchronizedStream.cpp @@ -92,6 +92,7 @@ void SynchronizedStream::append(message::Level msgLevel, message, tagName, std::to_string(rank), + "1", fileName, line)); } From c845143e133d10aefc30f8f7db3fc38b430ca8e4 Mon Sep 17 00:00:00 2001 From: Chris White Date: Thu, 25 Jul 2024 14:45:05 -0700 Subject: [PATCH 02/27] remove old code --- src/axom/slic/examples/basic/lumberjack_logging.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/axom/slic/examples/basic/lumberjack_logging.cpp b/src/axom/slic/examples/basic/lumberjack_logging.cpp index 9bcb591a59..fdf409af49 100644 --- a/src/axom/slic/examples/basic/lumberjack_logging.cpp +++ b/src/axom/slic/examples/basic/lumberjack_logging.cpp @@ -48,10 +48,6 @@ constexpr const char* format = R"( \tFILE= \tLINE= )"; - // std::string format = std::string("\n") + - // std::string("\t\n") + std::string("\tLEVEL=\n") + - // std::string("\tRANKS=\n") + std::string("\tFILE=\n") + - // std::string("\tLINE=\n"); slic::initialize(); // Set SLIC logging level and Lumberjack Logging stream From 6d66fa9c40ea4f9a6e6e50683bfe65fb322b107b Mon Sep 17 00:00:00 2001 From: format-robot Date: Thu, 25 Jul 2024 14:46:43 -0700 Subject: [PATCH 03/27] Apply style updates --- src/axom/slic/examples/basic/lumberjack_logging.cpp | 2 +- src/axom/slic/streams/LumberjackStream.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/axom/slic/examples/basic/lumberjack_logging.cpp b/src/axom/slic/examples/basic/lumberjack_logging.cpp index fdf409af49..c8d8e080bf 100644 --- a/src/axom/slic/examples/basic/lumberjack_logging.cpp +++ b/src/axom/slic/examples/basic/lumberjack_logging.cpp @@ -39,7 +39,7 @@ int main(int argc, char** argv) MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Initialize SLIC -constexpr const char* format = R"( + constexpr const char* format = R"( \t \tLEVEL= diff --git a/src/axom/slic/streams/LumberjackStream.cpp b/src/axom/slic/streams/LumberjackStream.cpp index 42d24e3a9c..0f26c430fa 100644 --- a/src/axom/slic/streams/LumberjackStream.cpp +++ b/src/axom/slic/streams/LumberjackStream.cpp @@ -141,7 +141,8 @@ void LumberjackStream::write(bool local) { for(const auto* curr_message : m_lj->getMessages()) { - if (curr_message == nullptr) { + if(curr_message == nullptr) + { continue; } From 42c7cba145c8ebe287fb85e20fefba1f3fb422cd Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 12 Apr 2024 16:18:58 -0700 Subject: [PATCH 04/27] Adds python script to convert SVGs to MFEM meshes --- src/tools/svg2contours/README.md | 12 ++ src/tools/svg2contours/requirements.txt | 3 + src/tools/svg2contours/svg2contours.ipynb | 185 ++++++++++++++++++++++ 3 files changed, 200 insertions(+) create mode 100644 src/tools/svg2contours/README.md create mode 100644 src/tools/svg2contours/requirements.txt create mode 100644 src/tools/svg2contours/svg2contours.ipynb diff --git a/src/tools/svg2contours/README.md b/src/tools/svg2contours/README.md new file mode 100644 index 0000000000..a5297b6e93 --- /dev/null +++ b/src/tools/svg2contours/README.md @@ -0,0 +1,12 @@ +### Create a virtual environment + +```shell +python3 -m venv venv + +# linux +source venv/bin/activate +# windows bash +source ./venv/Scripts/activate + +pip3 install -r requirements.txt +``` \ No newline at end of file diff --git a/src/tools/svg2contours/requirements.txt b/src/tools/svg2contours/requirements.txt new file mode 100644 index 0000000000..35f1f1c477 --- /dev/null +++ b/src/tools/svg2contours/requirements.txt @@ -0,0 +1,3 @@ +svgpathtools +jupyter +# mfem not yet available on windows \ No newline at end of file diff --git a/src/tools/svg2contours/svg2contours.ipynb b/src/tools/svg2contours/svg2contours.ipynb new file mode 100644 index 0000000000..f76b07d410 --- /dev/null +++ b/src/tools/svg2contours/svg2contours.ipynb @@ -0,0 +1,185 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "input_file = \"drawing.svg\"\n", + "# input_file = \"transform.svg\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from svgpathtools import Document, Path, Line, QuadraticBezier, CubicBezier, Arc, is_bezier_path, svg2paths, wsvg\n", + "\n", + "doc = Document(input_file)\n", + "paths = doc.paths()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def lerp(a,b,t):\n", + " return (1-t)*a+t*b\n", + "\n", + "def line_to_cubic(line : Line):\n", + " p_0,p_3 = line.bpoints()\n", + " return CubicBezier(p_0, lerp(p_0, p_3, 1/3), lerp(p_0, p_3, 2/3), p_3)\n", + "\n", + "def quadratic_to_cubic(quad : QuadraticBezier):\n", + " q_0,q_1,q_2 = quad.bpoints()\n", + " return CubicBezier(q_0, lerp(q_0,q_1, 2/3), lerp(q_1,q_2, 1/3), q_2)\n", + "\n", + "def segment_as_cubic(seg):\n", + " if isinstance(seg,Line):\n", + " return line_to_cubic(seg)\n", + " elif isinstance(seg,QuadraticBezier):\n", + " return quadratic_to_cubic(seg)\n", + " elif isinstance(seg, CubicBezier):\n", + " return seg\n", + " elif isinstance(seg,Arc):\n", + " raise Exception(\"'Arc' type not yet supported\")\n", + " else:\n", + " raise Exception(f\"'{type(seg)}' type not supported yet\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for p_idx, p in enumerate(paths):\n", + " for seg_idx, seg in enumerate(p):\n", + " try:\n", + " cubic = segment_as_cubic(seg)\n", + " print(f\"[Path {p_idx}; Seg {seg_idx}]:\\n\\t{seg}\\n\\tas cubic: {cubic}\")\n", + " except Exception as err:\n", + " print(f\"[Path {p_idx}; Seg {seg_idx}]:\\n\\t{seg}\")\n", + " print(f\"parsed unsupported type {type(seg)}. Msg={err}\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "# Create an mfem mesh\n", + "\n", + "header = \"\"\"\n", + "MFEM mesh v1.0\n", + "\n", + "# MFEM Geometry Types (see mesh/geom.hpp):\n", + "#\n", + "# POINT = 0\n", + "# SEGMENT = 1\n", + "# TRIANGLE = 2\n", + "# SQUARE = 3\n", + "# TETRAHEDRON = 4\n", + "# CUBE = 5\n", + "#\n", + "\n", + "dimension\n", + "1\n", + "\"\"\"\n", + "\n", + "vert_cnt = 0\n", + "elem_cnt = 0\n", + "verts = []\n", + "dofs = []\n", + "elts = []\n", + "\n", + "print(paths)\n", + "\n", + "for p_idx, p in enumerate(paths):\n", + "\n", + " if not is_bezier_path(p):\n", + " continue\n", + "\n", + " for seg_idx, seg in enumerate(p):\n", + " cubic = segment_as_cubic(seg)\n", + " elts.append(\" \".join(map(str,[p_idx + 1, 1, vert_cnt, vert_cnt + 1])))\n", + " verts.append(\" \".join(map(str,[cubic.start.real, cubic.start.imag])))\n", + " verts.append(\" \".join(map(str,[cubic.end.real, cubic.end.imag])))\n", + " dofs.append(\" \".join(map(str,[cubic.control1.real, cubic.control1.imag])))\n", + " dofs.append(\" \".join(map(str,[cubic.control2.real, cubic.control2.imag])))\n", + " vert_cnt += 2\n", + " elem_cnt += 1 \n", + "\n", + "mfem_file = []\n", + "mfem_file.append(header)\n", + "mfem_file.append(f\"\"\"\n", + "elements\n", + "{elem_cnt}\n", + "{\"\\n\".join(elts)}\n", + "\"\"\")\n", + "\n", + "mfem_file.append(f\"\"\"\n", + "boundary\n", + "0\n", + "\"\"\")\n", + "\n", + "mfem_file.append(f\"\"\"\n", + "vertices\n", + "{vert_cnt}\n", + "\"\"\")\n", + "\n", + "mfem_file.append(f\"\"\"\n", + "nodes\n", + "FiniteElementSpace\n", + "FiniteElementCollection: H1Pos_1D_P3\n", + "VDim: 2\n", + "Ordering: 1\n", + "\n", + "{\"\\n\".join(verts)}\n", + "{\"\\n\".join(dofs)}\n", + "\"\"\")\n", + "\n", + "fname = \"drawing.mesh\"\n", + "with open(fname, mode='w') as f:\n", + " f.write(\"\\n\".join(mfem_file))\n", + " print(f\"wrote '{fname}' with {vert_cnt} vertices and {elem_cnt} elements\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From bac696bceafe31b40dd73863d27a565a5ae76595 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 12 Apr 2024 16:23:42 -0700 Subject: [PATCH 05/27] Adds an fmt formatter to primal's BezierCurve class --- src/axom/primal/geometry/BezierCurve.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 8aa0e47c3e..c21c7e2d30 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -26,6 +26,7 @@ #include #include +#include "axom/fmt.hpp" namespace axom { @@ -919,4 +920,10 @@ std::ostream& operator<<(std::ostream& os, const BezierCurve& bCurve) } // namespace primal } // namespace axom +/// Overload to format a primal::BezierCurve using fmt +template +struct axom::fmt::formatter> + : ostream_formatter +{ }; + #endif // AXOM_PRIMAL_BEZIERCURVE_HPP_ From 570826fe906b16924c9b2fa16a3c32c97db81d46 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 12 Apr 2024 16:22:41 -0700 Subject: [PATCH 06/27] Adds a winding number driver example to quest --- src/axom/quest/examples/CMakeLists.txt | 12 + .../quest/examples/quest_winding_number.cpp | 221 ++++++++++++++++++ 2 files changed, 233 insertions(+) create mode 100644 src/axom/quest/examples/quest_winding_number.cpp diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index d462ffff2a..97242b8035 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -532,3 +532,15 @@ if (ENABLE_FORTRAN) endif() endif() endif() + +# Quest winding number example ------------------------------------------------ + +if( MFEM_FOUND) + axom_add_executable( + NAME quest_winding_number_ex + SOURCES quest_winding_number.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON axom mfem cli11 fmt + FOLDER axom/quest/examples + ) +endif() diff --git a/src/axom/quest/examples/quest_winding_number.cpp b/src/axom/quest/examples/quest_winding_number.cpp new file mode 100644 index 0000000000..9361fe9097 --- /dev/null +++ b/src/axom/quest/examples/quest_winding_number.cpp @@ -0,0 +1,221 @@ +// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file quest_winding_number.cpp + * \brief Example that computes the winding number of a grid of points + * against a collection of 2D parametric curves + */ + +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/primal.hpp" +#include "axom/quest.hpp" + +#include "axom/CLI11.hpp" +#include "axom/fmt.hpp" + +#include "mfem.hpp" + +namespace primal = axom::primal; +using Point2D = primal::Point; +using BezierCurve2D = primal::BezierCurve; +using BoundingBox2D = primal::BoundingBox; + +BezierCurve2D segment_to_curve(const mfem::Mesh* mesh, int elem_id) +{ + const auto* fes = mesh->GetNodes()->FESpace(); + const auto* fec = fes->FEColl(); + //auto* nodes = mesh->GetNodes(); + + mfem::Array vdofs; + mfem::Vector v; + + fes->GetElementVDofs(elem_id, vdofs); + mesh->GetNodes()->GetSubVector(vdofs, v); + + //for(int j = 0; j < vdofs.Size(); ++j) + //{ + // SLIC_INFO(axom::fmt::format("Element {} -- vdof {}", elem_id, vdofs[j])); + //} + + // temporarily hard-code for 3rd order + axom::Array points(4, 4); + points[0] = Point2D {v[0], v[0 + 4]}; + points[1] = Point2D {v[2], v[2 + 4]}; + points[2] = Point2D {v[3], v[3 + 4]}; + points[3] = Point2D {v[1], v[1 + 4]}; + + // // TODO: Handle interleaved points + // const int stride = vdofs.Size() / 2; + // for(int n = 0; n < stride; ++n) + // { + // points.push_back(Point2D {v[n], v[n + stride]}); + // } + + return BezierCurve2D(points, fec->GetOrder()); +} + +bool check_mesh_valid(const mfem::Mesh* mesh) +{ + const auto* fes = mesh->GetNodes()->FESpace(); + const auto* fec = fes->FEColl(); + + // TODO: Convert from Lagrange to Bernstein, if/when necessary + if(!dynamic_cast(fec)) + { + SLIC_WARNING( + "Example only currently supports meshes with nodes in the Bernstein " + "basis"); + return false; + } + + if(fes->GetVDim() != 2) + { + SLIC_WARNING("Example only currently supports 2D meshes"); + return false; + } + + return true; +} + +int main(int argc, char** argv) +{ + axom::slic::SimpleLogger raii_logger; + + axom::CLI::App app { + "Load mesh containing collection of curves" + " and optionally generate a query mesh of winding numbers."}; + + std::string filename = + std::string(AXOM_SRC_DIR) + "/tools/svg2contours/drawing.mesh"; + + bool verbose {false}; + + // Query mesh parameters + std::vector boxMins; + std::vector boxMaxs; + std::vector boxResolution; + + app.add_option("-f,--file", filename) + ->description("Mfem mesh containing contours") + ->check(axom::CLI::ExistingFile); + + app.add_flag("-v,--verbose", verbose, "verbose output")->capture_default_str(); + + auto* inline_mesh_subcommand = + app.add_subcommand("inline_mesh") + ->description("Options for setting up a query mesh") + ->fallthrough(); + + inline_mesh_subcommand->add_option("--min", boxMins) + ->description("Min bounds for box mesh (x,y)") + ->expected(2) + ->required(); + inline_mesh_subcommand->add_option("--max", boxMaxs) + ->description("Max bounds for box mesh (x,y)") + ->expected(2) + ->required(); + + inline_mesh_subcommand->add_option("--res", boxResolution) + ->description("Resolution of the box mesh (i,j)") + ->expected(2) + ->required(); + + CLI11_PARSE(app, argc, argv); + + mfem::Mesh mesh(filename); + SLIC_INFO( + axom::fmt::format("Curve mesh has a topological dimension of {}d, has {} " + "vertices and {} elements\n", + mesh.Dimension(), + mesh.GetNV(), + mesh.GetNE())); + + if(!check_mesh_valid(&mesh)) + { + return 1; + } + + axom::Array segments; + axom::Array curves; + + // only retain the 1D segments + for(int i = 0; i < mesh.GetNE(); ++i) + { + auto* el = mesh.GetElement(i); + if(el->GetGeometryType() == mfem::Geometry::SEGMENT) + { + segments.push_back(i); + } + } + + BoundingBox2D bbox; + for(int i = 0; i < segments.size(); ++i) + { + auto curve = segment_to_curve(&mesh, i); + SLIC_INFO_IF(verbose, axom::fmt::format("Element {}: {}", i, curve)); + + bbox.addBox(curve.boundingBox()); + + curves.emplace_back(std::move(curve)); + } + + SLIC_INFO(axom::fmt::format("Curve mesh bounding box: {}", bbox)); + + // Early return if user didn't set up a query mesh + if(boxMins.empty()) + { + return 0; + } + + const auto query_res = primal::NumericArray(boxResolution.data()); + const auto query_box = + BoundingBox2D(Point2D(boxMins.data()), Point2D(boxMaxs.data())); + constexpr int query_order = 1; + + auto query_mesh = std::unique_ptr( + axom::quest::util::make_cartesian_mfem_mesh_2D(query_box, + query_res, + query_order)); + auto fec = mfem::H1_FECollection(query_order, 2); + auto fes = mfem::FiniteElementSpace(query_mesh.get(), &fec, 1); + auto winding = mfem::GridFunction(&fes); + auto inout = mfem::GridFunction(&fes); + + for(int vidx = 0; vidx < query_mesh->GetNV(); ++vidx) + { + Point2D q(query_mesh->GetVertex(vidx), 2); + // SLIC_INFO(axom::fmt::format("Vertex {} -- {}", vidx, q)); + + double wn {}; + for(const auto& c : curves) + { + wn += axom::primal::winding_number(q, c); + } + + winding[vidx] = wn; + inout[vidx] = std::round(wn); + + // SLIC_INFO(axom::fmt::format( + // "Winding number for query point {} is {} -- rounded to {}", + // q, + // wn, + // std::round(wn))); + } + + std::string output_name = "winding"; + mfem::VisItDataCollection windingDC(output_name, query_mesh.get()); + windingDC.RegisterField("winding", &winding); + windingDC.RegisterField("inout", &inout); + windingDC.Save(); + + SLIC_INFO(axom::fmt::format("Outputting generated mesh '{}' to '{}'", + output_name, + axom::utilities::filesystem::getCWD())); + + return 0; +} \ No newline at end of file From 1c926677c23af70b1bfb50c51b58b71d9c2f6d28 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 12 Apr 2024 18:43:27 -0700 Subject: [PATCH 07/27] Updates the jupyter notebook on linux --- src/tools/svg2contours/README.md | 43 ++++++++++-- src/tools/svg2contours/requirements.txt | 2 +- src/tools/svg2contours/svg2contours.ipynb | 86 ++++++++++++++++------- 3 files changed, 100 insertions(+), 31 deletions(-) diff --git a/src/tools/svg2contours/README.md b/src/tools/svg2contours/README.md index a5297b6e93..70ffe725fa 100644 --- a/src/tools/svg2contours/README.md +++ b/src/tools/svg2contours/README.md @@ -1,12 +1,45 @@ ### Create a virtual environment ```shell -python3 -m venv venv +> python3 -m venv venv # linux -source venv/bin/activate +> source venv/bin/activate # windows bash -source ./venv/Scripts/activate +> source ./venv/Scripts/activate -pip3 install -r requirements.txt -``` \ No newline at end of file +> pip3 install -r requirements.txt +``` + +### Open jupyter notebook in vs code + +Need to run the notebook using a kernel from this venv. + +#### This might work + +* Open the ipynb file in vs code with jupyter extension +* Point jupyter to the ipython3 interpretter from your environment +``` +CTRL+SHIFT+P -> "python select interpretter" + +[svgpathtools@1.6.1](https://github.com/mathandy/svgpathtools/releases/tag/v1.6.1) has a bug in applying the correct rotation angle +when transforming ellipses and elliptical arcs. +This can be resolved by applying the following patch: +```shell +> patch -p1 venv/lib/python3.9/site-packages/svgpathtools/path.py -i svgpathtools-1.6.1-ellipse-rotation.patch --verbose ``` -* Note: Might need to restart vs code to get this working - -< based on https://stackoverflow.com/a/74410917 > - -#### If not, try to create a kernel +See: https://github.com/mathandy/svgpathtools/pull/221 +#### Developer's note: +The patch was generated from a git commit from the above pull request using the following command: ```shell -> python3 -m ipykernel install --user --name=svg2contour +> git format-patch -1 260a44ed2c0d114f77d57016d6d143a50729aca9 --stdout > svgpathtools-1.6.1-ellipse-rotation.patch ``` -Then point jupyter to this kernel via "select kernel" button +### Run the script on an input SVG mesh -### Run the notebook on your SVG file -This generates an mfem `*.mesh` file +To convert an SVG to an mfem NURBS mesh, run the following command: +```shell +> cd / +> ../src/tools/svg2contours/svg2contours.py -i ../data/contours/svg/shapes.svg + +SVG dimensions: width='210mm' height='297mm' viewBox='0 0 210 297' +Wrote 'drawing.mesh' with 54 vertices and NURBS 27 elements +``` +Note: This assumes your axom clone has the `data` submodule located +at `/data` ### Run the quest winding number example +Now that we have an MFEM NURBS mesh, we can run our winding number application ```shell ->./examples/quest_winding_number_ex -f ../data/contours/svg/drawing.mesh -v inline_mesh --min 0 0 --max 250 250 --res 500 500 +> cd // +> ./examples/quest_winding_number_ex \ + -f ./drawing.mesh \ + query_mesh --min 0 0 --max 250 250 --res 500 500 ``` diff --git a/src/tools/svg2contours/svgpathtools-1.6.1-ellipse-rotation.patch b/src/tools/svg2contours/svgpathtools-1.6.1-ellipse-rotation.patch new file mode 100644 index 0000000000..cc9b26b920 --- /dev/null +++ b/src/tools/svg2contours/svgpathtools-1.6.1-ellipse-rotation.patch @@ -0,0 +1,24 @@ +From 260a44ed2c0d114f77d57016d6d143a50729aca9 Mon Sep 17 00:00:00 2001 +Subject: [PATCH] Bugfix in rotation angle calculation when transforming arcs + +Use arctan to properly account for quadrant of angle. +--- + svgpathtools/path.py | 2 +- + 1 file changed, 1 insertion(+), 1 deletion(-) + +diff --git a/svgpathtools/path.py b/svgpathtools/path.py +index 8ccab5c..87fba5e 100644 +--- a/svgpathtools/path.py ++++ b/svgpathtools/path.py +@@ -334,7 +334,7 @@ def transform(curve, tf): + new_radius = complex(rx, ry) + + xeigvec = eigvecs[:, 0] +- rot = np.degrees(np.arccos(xeigvec[0])) ++ rot = np.degrees(np.arctan2(xeigvec[1], xeigvec[0])) + + if new_radius.real == 0 or new_radius.imag == 0 : + return Line(new_start, new_end) +-- +2.29.1 + From 918cbb05effb15f65231f67b562ea90929992556 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 14 Jun 2024 20:07:40 -0700 Subject: [PATCH 21/27] Improves comments and error checks in quest_winding_number example --- .../quest/examples/quest_winding_number.cpp | 110 +++++++++++++----- src/tools/svg2contours/README.md | 2 +- 2 files changed, 83 insertions(+), 29 deletions(-) diff --git a/src/axom/quest/examples/quest_winding_number.cpp b/src/axom/quest/examples/quest_winding_number.cpp index 3de293fd88..c55c38a8aa 100644 --- a/src/axom/quest/examples/quest_winding_number.cpp +++ b/src/axom/quest/examples/quest_winding_number.cpp @@ -6,7 +6,9 @@ /*! * \file quest_winding_number.cpp * \brief Example that computes the winding number of a grid of points - * against a collection of 2D parametric curves + * against a collection of 2D parametric rational curves. + * Supports MFEM meshes in the cubic positive Bernstein basis or the (rational) + * NURBS basis. */ #include "axom/config.hpp" @@ -25,18 +27,37 @@ using Point2D = primal::Point; using BezierCurve2D = primal::BezierCurve; using BoundingBox2D = primal::BoundingBox; +/*! + * Given an mfem mesh, convert element with id \a elem_id to a (rational) BezierCurve + * \pre Assumes the elements of the mfem mesh are in the positive (Bernstein) + * basis, or in the NURBS basis + */ BezierCurve2D segment_to_curve(const mfem::Mesh* mesh, int elem_id) { const auto* fes = mesh->GetNodes()->FESpace(); const auto* fec = fes->FEColl(); - //auto* nodes = mesh->GetNodes(); const bool isBernstein = dynamic_cast(fec) != nullptr; const bool isNURBS = dynamic_cast(fec) != nullptr; - SLIC_ASSERT(isBernstein || isNURBS); + SLIC_ERROR_IF( + !(isBernstein || isNURBS), + "MFEM mesh elements must be in either the Bernstein or NURBS basis"); + + const int NE = isBernstein ? mesh->GetNE() : fes->GetNURBSext()->GetNP(); + SLIC_ERROR_IF(NE < elem_id, + axom::fmt::format("Mesh does not have {} elements", elem_id)); + + const int order = + isBernstein ? fes->GetOrder(elem_id) : mesh->NURBSext->GetOrders()[elem_id]; + SLIC_ERROR_IF(order != 3, + axom::fmt::format( + "This example currently requires the input mfem mesh to " + "contain cubic elements, but the order of element {} is {}", + elem_id, + order)); mfem::Array dofs; mfem::Array vdofs; @@ -48,12 +69,7 @@ BezierCurve2D segment_to_curve(const mfem::Mesh* mesh, int elem_id) fes->GetElementVDofs(elem_id, vdofs); mesh->GetNodes()->GetSubVector(vdofs, v); - for(int j = 0; j < vdofs.Size(); ++j) - { - SLIC_INFO(axom::fmt::format("Element {} -- vdof {}", elem_id, vdofs[j])); - } - - // temporarily hard-code for 3rd order + // Currently hard-coded for 3rd order. This can easily be extended to arbitrary order axom::Array points(4, 4); if(isBernstein) { @@ -84,7 +100,18 @@ BezierCurve2D segment_to_curve(const mfem::Mesh* mesh, int elem_id) bool check_mesh_valid(const mfem::Mesh* mesh) { const auto* fes = mesh->GetNodes()->FESpace(); + if(fes == nullptr) + { + SLIC_WARNING("MFEM mesh finite element space was null"); + return false; + } + const auto* fec = fes->FEColl(); + if(fec == nullptr) + { + SLIC_WARNING("MFEM mesh finite element collection was null"); + return false; + } const bool isBernstein = dynamic_cast(fec) != nullptr; @@ -107,6 +134,27 @@ bool check_mesh_valid(const mfem::Mesh* mesh) return false; } + const int NE = isBernstein ? mesh->GetNE() : fes->GetNURBSext()->GetNP(); + int order = -1; + if(isBernstein) + { + order = NE > 0 ? fes->GetOrder(0) : 3; + } + else // isNURBS + { + //SLIC_INFO("nurbsext order :" << mesh->NURBSext->GetOrder()); + order = NE > 0 ? mesh->NURBSext->GetOrders()[0] : 3; + } + + if(order != 3) + { + SLIC_WARNING(axom::fmt::format( + "This example currently requires the input mfem mesh to contain cubic " + "elements, but the provided mesh has order {}", + order)); + return false; + } + return true; } @@ -118,8 +166,8 @@ int main(int argc, char** argv) "Load mesh containing collection of curves" " and optionally generate a query mesh of winding numbers."}; - std::string filename = - std::string(AXOM_SRC_DIR) + "/tools/svg2contours/drawing.mesh"; + std::string inputFile; + std::string outputPrefix = {"winding"}; bool verbose {false}; @@ -129,10 +177,17 @@ int main(int argc, char** argv) std::vector boxResolution; int queryOrder {1}; - app.add_option("-f,--file", filename) - ->description("Mfem mesh containing contours") + app.add_option("-i,--input", inputFile) + ->description("MFEM mesh containing contours (1D segments)") + ->required() ->check(axom::CLI::ExistingFile); + app.add_option("-o,--output-prefix", outputPrefix) + ->description( + "Prefix for output 2D query mesh (in MFEM format) mesh containing " + "winding number calculations") + ->capture_default_str(); + app.add_flag("-v,--verbose", verbose, "verbose output")->capture_default_str(); auto* query_mesh_subcommand = @@ -151,16 +206,16 @@ int main(int argc, char** argv) ->description("Resolution of the box mesh (i,j)") ->expected(2) ->required(); - query_mesh_subcommand->add_option("-o,--order", queryOrder) + query_mesh_subcommand->add_option("--order", queryOrder) ->description("polynomial order of the query mesh") ->check(axom::CLI::PositiveNumber); CLI11_PARSE(app, argc, argv); - mfem::Mesh mesh(filename); + mfem::Mesh mesh(inputFile); SLIC_INFO( - axom::fmt::format("Curve mesh has a topological dimension of {}d, has {} " - "vertices and {} elements\n", + axom::fmt::format("Curve mesh has a topological dimension of {}d, " + "has {} vertices and {} elements", mesh.Dimension(), mesh.GetNV(), mesh.GetNE())); @@ -173,7 +228,7 @@ int main(int argc, char** argv) axom::Array segments; axom::Array curves; - // only retain the 1D segments + // Loop through mesh elements, retaining the (curved) 1D segments for(int i = 0; i < mesh.GetNE(); ++i) { auto* el = mesh.GetElement(i); @@ -183,6 +238,7 @@ int main(int argc, char** argv) } } + // Extract the curves and compute their bounding boxes along the way BoundingBox2D bbox; for(int i = 0; i < segments.size(); ++i) { @@ -202,6 +258,7 @@ int main(int argc, char** argv) return 0; } + // Generate a Cartesian (high order) mesh for the query points const auto query_res = primal::NumericArray(boxResolution.data()); const auto query_box = BoundingBox2D(Point2D(boxMins.data()), Point2D(boxMaxs.data())); @@ -216,11 +273,14 @@ int main(int argc, char** argv) auto inout = mfem::GridFunction(&fes); auto nodes_fes = query_mesh->GetNodalFESpace(); + // Query the winding numbers at each degree of freedom (DoF) of the query mesh. + // The loop below independently checks (and adaptively refines) every curve for each query points. + // A more efficient algorithm can de defined that caches the refined curves to avoid + // extra refinements. We will add this in a follow-up PR. for(int nidx = 0; nidx < nodes_fes->GetNDofs(); ++nidx) { Point2D q; query_mesh->GetNode(nidx, q.data()); - // SLIC_INFO(axom::fmt::format("Node {} -- {}", nidx, q)); double wn {}; for(const auto& c : curves) @@ -230,22 +290,16 @@ int main(int argc, char** argv) winding[nidx] = wn; inout[nidx] = std::round(wn); - - // SLIC_INFO(axom::fmt::format( - // "Winding number for query point {} is {} -- rounded to {}", - // q, - // wn, - // std::round(wn))); } - std::string output_name = "winding"; - mfem::VisItDataCollection windingDC(output_name, query_mesh.get()); + // Save the query mesh and fields to disk using a format that can be viewed in VisIt + mfem::VisItDataCollection windingDC(outputPrefix, query_mesh.get()); windingDC.RegisterField("winding", &winding); windingDC.RegisterField("inout", &inout); windingDC.Save(); SLIC_INFO(axom::fmt::format("Outputting generated mesh '{}' to '{}'", - output_name, + windingDC.GetCollectionName(), axom::utilities::filesystem::getCWD())); return 0; diff --git a/src/tools/svg2contours/README.md b/src/tools/svg2contours/README.md index 7aed2a4f0a..5c4c3906c6 100644 --- a/src/tools/svg2contours/README.md +++ b/src/tools/svg2contours/README.md @@ -55,6 +55,6 @@ Now that we have an MFEM NURBS mesh, we can run our winding number application ```shell > cd // > ./examples/quest_winding_number_ex \ - -f ./drawing.mesh \ + -i ./drawing.mesh \ query_mesh --min 0 0 --max 250 250 --res 500 500 ``` From 115ada66dcd2d62190648ba57a4bf6fb1dd54eee Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 21 Jun 2024 17:16:02 -0700 Subject: [PATCH 22/27] Adds another patch file for svgpathtools w/ updated instructions in the README --- src/tools/svg2contours/README.md | 18 ++++++++----- .../svgpathtools-1.6.1-rounded-rect.patch | 26 +++++++++++++++++++ 2 files changed, 38 insertions(+), 6 deletions(-) create mode 100644 src/tools/svg2contours/svgpathtools-1.6.1-rounded-rect.patch diff --git a/src/tools/svg2contours/README.md b/src/tools/svg2contours/README.md index 5c4c3906c6..fb5319eff0 100644 --- a/src/tools/svg2contours/README.md +++ b/src/tools/svg2contours/README.md @@ -1,6 +1,6 @@ ## Running the svg2contours Python script -The `svg2countours` script converts [SVG](https://developer.mozilla.org/en-US/docs/Web/SVG) images to [MFEM NURBS meshes](https://mfem.org/mesh-format-v1.0/#nurbs-meshes) using the [svgpathtools](https://github.com/mathandy/svgpathtools) library. +The `svg2contours` script converts [SVG](https://developer.mozilla.org/en-US/docs/Web/SVG) images to [MFEM NURBS meshes](https://mfem.org/mesh-format-v1.0/#nurbs-meshes) using the [svgpathtools](https://github.com/mathandy/svgpathtools) library. The latter can be used with axom's `quest_winding_number` example application to sample the winding number field over the generated curves. @@ -20,20 +20,26 @@ Full SVG support requires a (slightly) patched copy of [svgpathtools@1.6.1](http > pip3 install -r requirements.txt ``` -### Apply patch to svgpathtools for proper treatment of rotated ellipse +### Apply patches to svgpathtools for proper treatment of rotated ellipses and rounded rectangles -[svgpathtools@1.6.1](https://github.com/mathandy/svgpathtools/releases/tag/v1.6.1) has a bug in applying the correct rotation angle -when transforming ellipses and elliptical arcs. +[svgpathtools@1.6.1](https://github.com/mathandy/svgpathtools/releases/tag/v1.6.1) has a bug in applying the correct rotation angle when transforming ellipses and elliptical arcs. This can be resolved by applying the following patch: ```shell > patch -p1 venv/lib/python3.9/site-packages/svgpathtools/path.py -i svgpathtools-1.6.1-ellipse-rotation.patch --verbose ``` See: https://github.com/mathandy/svgpathtools/pull/221 +It has another bug related to rounded rectangles, which can be resolved by applying the following pathc: +```shell +> patch -p1 venv/lib/python3.9/site-packages/svgpathtools/svg_to_paths.py -i svgpathtools-1.6.1-rounded-rect.patch --verbose +``` +See: https://github.com/mathandy/svgpathtools/pull/222 + #### Developer's note: -The patch was generated from a git commit from the above pull request using the following command: +These patches were generated from a git commit in each of the above pull requests using the following commands: ```shell -> git format-patch -1 260a44ed2c0d114f77d57016d6d143a50729aca9 --stdout > svgpathtools-1.6.1-ellipse-rotation.patch + > git format-patch -1 260a44ed2c0d114f77d57016d6d143a50729aca9 --stdout > svgpathtools-1.6.1-ellipse-rotation.patch + > git format-patch -1 ec1e1101037fcd66967caa40bc2b038c928bae4f --stdout > svgpathtools-1.6.1-rounded-rect.patch ``` ### Run the script on an input SVG mesh diff --git a/src/tools/svg2contours/svgpathtools-1.6.1-rounded-rect.patch b/src/tools/svg2contours/svgpathtools-1.6.1-rounded-rect.patch new file mode 100644 index 0000000000..bbcdfd2a98 --- /dev/null +++ b/src/tools/svg2contours/svgpathtools-1.6.1-rounded-rect.patch @@ -0,0 +1,26 @@ +From ec1e1101037fcd66967caa40bc2b038c928bae4f Mon Sep 17 00:00:00 2001 +Subject: [PATCH] Bugfix for converting rounded rect to a d-string path + +The previous algorithm supported the case where the rect was a python dictionary, +but not when it was a rect Element instance. +--- + svgpathtools/svg_to_paths.py | 3 ++- + 1 file changed, 2 insertions(+), 1 deletion(-) + +diff --git a/svgpathtools/svg_to_paths.py b/svgpathtools/svg_to_paths.py +index 65591af..c6ac6df 100644 +--- a/svgpathtools/svg_to_paths.py ++++ b/svgpathtools/svg_to_paths.py +@@ -93,7 +93,8 @@ def rect2pathd(rect): + rectangle object and proceed counter-clockwise.""" + x, y = float(rect.get('x', 0)), float(rect.get('y', 0)) + w, h = float(rect.get('width', 0)), float(rect.get('height', 0)) +- if 'rx' in rect or 'ry' in rect: ++ ++ if 'rx' in rect.keys() or 'ry' in rect.keys(): + + # if only one, rx or ry, is present, use that value for both + # https://developer.mozilla.org/en-US/docs/Web/SVG/Element/rect +-- +2.29.1 + From 5843e83d09b9b886ffbc1fed05c89f99d74bbd89 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 21 Jun 2024 17:49:42 -0700 Subject: [PATCH 23/27] Bugfix: Need to check element.keys() for attributes like `d` --- src/tools/svg2contours/svg2contours.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index cbc39a81e3..1f08d915c4 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -332,9 +332,9 @@ def main(): mfem_data = MFEMData() for p_idx, p in enumerate(paths): - # print(f"""reading {p_idx=} {p=} \n w/ {p.d()=}""") + #print(f"""reading {p_idx=} {p=} \n w/ {p.d()=}""") - is_d_path = 'd' in p.element + is_d_path = 'd' in p.element.keys() attrib = p_idx + 1 reverse_paths = True if np.linalg.det(coordinate_transform) < 0 else False @@ -343,6 +343,7 @@ def main(): continue 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: # split large elliptical arcs for easier processing @@ -371,4 +372,4 @@ def main(): if __name__ == "__main__": exitcode = main() sys.exit(exitcode) - \ No newline at end of file + From 147a06d1c7391b9975627ac30a0435c424309e9e Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Fri, 21 Jun 2024 17:56:27 -0700 Subject: [PATCH 24/27] Formats svg2contours script (using `black --line-length 120`) --- src/tools/svg2contours/svg2contours.py | 303 ++++++++++++------------- 1 file changed, 151 insertions(+), 152 deletions(-) diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index 1f08d915c4..880119191d 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -16,36 +16,50 @@ import sys import os -from svgpathtools import Document, Path, Line, QuadraticBezier, CubicBezier, Arc, is_bezier_path, is_bezier_segment, is_path_segment, svg2paths, bpoints2bezier +from svgpathtools import ( + Document, + Path, + Line, + QuadraticBezier, + CubicBezier, + Arc, + is_bezier_path, + is_bezier_segment, + is_path_segment, + svg2paths, + bpoints2bezier, +) import numpy as np import re import argparse + def get_root_transform(doc: Document): """Create transform to convert 2D coordinate system from y pointing down to y pointing up""" - + attr = doc.tree.getroot().attrib - width = attr.get('width', None) - height = attr.get('height', None) + width = attr.get("width", None) + height = attr.get("height", None) viewBox = attr.get("viewBox", None) print(f"""SVG dimensions: {width=} {height=} {viewBox=}""") tf = np.identity(3) - if height and not height.endswith('%'): + if height and not height.endswith("%"): match = re.search(r"(\d+)", height.strip()) if match: height_number = match.group(0) - tf[1,1] = -1 - tf[1,2] = height_number + tf[1, 1] = -1 + tf[1, 2] = height_number return tf -def transform_cubic(cubic : CubicBezier, tf): + +def transform_cubic(cubic: CubicBezier, tf): """Apply transformation `tf` to control points of cubic Bezier. - Adapted from internal `transform()` function in svgpathtools library + Adapted from internal `transform()` function in svgpathtools library """ def to_point(p): @@ -54,22 +68,25 @@ def to_point(p): def to_complex(v): return v.item(0) + 1j * v.item(1) - return bpoints2bezier([to_complex(tf.dot(to_point(p))) - for p in cubic.bpoints()]) + return bpoints2bezier([to_complex(tf.dot(to_point(p))) for p in cubic.bpoints()]) -def lerp(a,b,t): + +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 + return (1 - t) * a + t * b + -def line_to_cubic(line : Line): +def line_to_cubic(line: Line): """Converts from an svgpathtools Line to a (rational) CubicBezier (with unit weights)""" - q_0,q_1 = line.bpoints() - return (CubicBezier(q_0, lerp(q_0, q_1, 1/3), lerp(q_0, q_1, 2/3), q_1), [1,1,1,1]) + q_0, q_1 = line.bpoints() + return (CubicBezier(q_0, lerp(q_0, q_1, 1.0 / 3), lerp(q_0, q_1, 2.0 / 3), q_1), [1, 1, 1, 1]) + -def quadratic_to_cubic(quad : QuadraticBezier): +def quadratic_to_cubic(quad: QuadraticBezier): """Converts from an svgpathtools QuadraticBezier to a (rational) CubicBezier (with unit weights)""" - q_0,q_1,q_2 = quad.bpoints() - return (CubicBezier(q_0, lerp(q_0,q_1, 2/3), lerp(q_1,q_2, 1/3), q_2), [1,1,1,1]) + q_0, q_1, q_2 = quad.bpoints() + return (CubicBezier(q_0, lerp(q_0, q_1, 2.0 / 3), lerp(q_1, q_2, 1.0 / 3), q_2), [1, 1, 1, 1]) + def arc_to_cubic(arc: Arc): """Convertes from an svgpathtools Arc to a rational CubicBezier""" @@ -78,13 +95,13 @@ def arc_to_cubic(arc: Arc): def area(p1, p2, p3): """Computes the area of a triangle defined by vertices p1, p2 and p3""" - v21,v31 = p2-p1, p3-p1 + v21, v31 = p2 - p1, p3 - p1 return 0.5 * np.abs(v21.real * v31.imag - v21.imag * v31.real) # Notes: # (1) We have control point positions 0 and 3 and their tangents on the ellipse # as well as the midpoint; we need to find control points 1 and 2 - # (2) We are computing these as the intersections of the tangent lines, with lines connecting + # (2) We are computing these as the intersections of the tangent lines, with lines connecting # the opposite endpoint to the midpoints # (3) The weights are then derived through an isometry with the semicirle, # whose weights are proportional to [3,1,1,3] @@ -93,21 +110,21 @@ def area(p1, p2, p3): # these lines contain the internal control points c_1 and c_2 scale_fac = 10 - shoulder = arc.point(.5) + shoulder = arc.point(0.5) d_0 = arc.derivative(0) d_1 = arc.derivative(1) # print(f"""For arc {arc} w/ {d_0=} and {d_1=};\n {arc.theta=}, {arc.phi=}, {arc.rotation=}, {arc.delta=}""") # extend the line segment from q3 to shoulder point # and find the intersection w/ tangent line @ control point 0 - l_3_1 = Line(q_3, q_3 + scale_fac*(shoulder-q_3)) + l_3_1 = Line(q_3, q_3 + scale_fac * (shoulder - q_3)) l_0_1 = Line(q_0, q_0 + d_0) ints_31_01 = l_3_1.intersect(l_0_1) # print(f"""Finding intersection @ control point 1\n {shoulder=}\n {l_3_1=}\n {l_0_1=}\n {ints_31_01=}""") # extend the line segment from q0 to shoulder point # and find the intersection w/ tangent line @ control point 3 - l_0_2 = Line(q_0, q_0 + scale_fac*(shoulder-q_0)) + l_0_2 = Line(q_0, q_0 + scale_fac * (shoulder - q_0)) l_3_2 = Line(q_3, q_3 - d_1) ints_02_32 = l_0_2.intersect(l_3_2) # print(f"""Finding intersection @ control point 2\n {shoulder=}\n {l_0_2=}\n {l_3_2=}\n {ints_02_32=}""") @@ -115,24 +132,23 @@ def area(p1, p2, p3): c_1 = l_0_1.point(ints_31_01[0][1]) c_2 = l_3_2.point(ints_02_32[0][1]) # print(f"""Control points from intersections: CP 1: {c_1}; CP 2 {c_2}""") - + # print(f"""\t""") # print(f"""\t""") - # TODO: Figure out how to set the reversed flag based on arc.sweep, and possibly a "flip" parameter reversed = not arc.sweep curve = CubicBezier(q_3, c_2, c_1, q_0) if reversed else CubicBezier(q_0, c_1, c_2, q_3) - + # compute the rational weights based on areas of triangles that skip the current index # formula is from "Shape factors and shoulder points for shape control of rational Bezier curves" # https://doi.org/10.1016/j.cad.2023.103477 - b0,b1,b2,b3 = curve.bpoints() - - areas = [area(b1,b2,b3), area(b0,b2,b3) / 3., area(b0,b1,b3) / 3., area(b1,b2,b3)] - shape_fac = [areas[1]**2/(areas[0]*areas[2]), areas[2]**2/(areas[1]*areas[3])] + b0, b1, b2, b3 = curve.bpoints() + + areas = [area(b1, b2, b3), area(b0, b2, b3) / 3.0, area(b0, b1, b3) / 3.0, area(b1, b2, b3)] + shape_fac = [areas[1] ** 2 / (areas[0] * areas[2]), areas[2] ** 2 / (areas[1] * areas[3])] weights = [1, 0, 0, 1] - weights[1] = np.cbrt(shape_fac[0]**2 * shape_fac[1]) - weights[2] = weights[1]**2 / shape_fac[0] + weights[1] = np.cbrt(shape_fac[0] ** 2 * shape_fac[1]) + weights[2] = weights[1] ** 2 / shape_fac[0] # print(f""" -- {weights=} and {shape_fac=} for {curve=}""") return (curve, weights) @@ -146,31 +162,31 @@ def area(p1, p2, p3): for c in arc.as_cubic_curves(): q_0, q_1 = c.start, c.end c_1, c_2 = c.control1, c.control2 - return (CubicBezier(q_0, c_1, c_2, q_1), [3,1,1,3]) + return (CubicBezier(q_0, c_1, c_2, q_1), [3, 1, 1, 3]) -def segment_as_cubic(seg, reverse_paths : bool): - if isinstance(seg,Line): - cubic,weights = line_to_cubic(seg) - elif isinstance(seg,QuadraticBezier): - cubic,weights = quadratic_to_cubic(seg) +def segment_as_cubic(seg, reverse_paths: bool): + if isinstance(seg, Line): + cubic, weights = line_to_cubic(seg) + elif isinstance(seg, QuadraticBezier): + cubic, weights = quadratic_to_cubic(seg) elif isinstance(seg, CubicBezier): - cubic,weights = seg, [1,1,1,1] - elif isinstance(seg,Arc): - cubic,weights = arc_to_cubic(seg) + cubic, weights = seg, [1, 1, 1, 1] + elif isinstance(seg, Arc): + cubic, weights = arc_to_cubic(seg) else: raise Exception(f"'{type(seg)}' type not supported yet") - + if reverse_paths: - cubic = cubic.reversed() - weights.reverse() + cubic = cubic.reversed() + weights.reverse() + + return (cubic, weights) - return (cubic,weights) - def dist_to_ellipse(center, radius, angle, pt): - cx,cy = center.real, center.imag - rx,ry = radius.real, radius.imag + cx, cy = center.real, center.imag + rx, ry = radius.real, radius.imag rot = np.exp(-1j * np.radians(angle)) transformed_pt = rot * complex(pt.real - cx, pt.imag - cy) @@ -194,110 +210,96 @@ def __init__(self): def add_cubic_bezier(self, cubic, weights, attrib): - self.elems.append(" ".join(map(str,[attrib, 1, self.vert_cnt, self.vert_cnt + 1]))) + self.elems.append(" ".join(map(str, [attrib, 1, self.vert_cnt, self.vert_cnt + 1]))) self.vert_cnt += 2 self.edges.append(f"{self.elem_cnt} 0 1") - self.elem_cnt += 1 + self.elem_cnt += 1 # Assume for now that the order is always 3 self.knots.append("3 4 0 0 0 0 1 1 1 1") self.wgts_ends.append(f"{weights[0]} {weights[3]}") self.wgts_ints.append(f"{weights[2]} {weights[1]}") - self.dof_ends.append(" ".join(map(str,[cubic.start.real, cubic.start.imag, cubic.end.real, cubic.end.imag]))) - self.dof_ints.append(" ".join(map(str,[cubic.control2.real, cubic.control2.imag, cubic.control1.real, cubic.control1.imag]))) + self.dof_ends.append(" ".join(map(str, [cubic.start.real, cubic.start.imag, cubic.end.real, cubic.end.imag]))) + self.dof_ints.append( + " ".join(map(str, [cubic.control2.real, cubic.control2.imag, cubic.control1.real, cubic.control1.imag])) + ) def write_file(self, filename): mfem_file = [] - - mfem_file.extend([ - "MFEM NURBS mesh v1.0", - "", - "# MFEM Geometry Types (see fem/geom.hpp):", - "#", - "# SEGMENT = 1 | SQUARE = 3 | CUBE = 5", - "#", - "# element: 1 ", - "# edge: 0 1 <-- idx increases by one each time", - "# knotvector: [knots]; sizeof(knots) is 1+order+num_ctrl_pts", - "# weights: array of weights corresponding to the NURBS element", - "# FES: list of control points; vertex control points at top, then interior control points", - ""]) - - 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([ - "knotvectors", - f"{self.elem_cnt}", - "\n".join(self.knots), - ""]) - - mfem_file.extend([ - "weights", - "\n".join(self.wgts_ends), - "\n".join(self.wgts_ints), - ""]) - - mfem_file.extend([ - "FiniteElementSpace", - "FiniteElementCollection: NURBS", - "VDim: 2", - "Ordering: 1", - "", - "\n".join(self.dof_ends), - "\n".join(self.dof_ints), - ""]) - - with open(filename, mode='w') as f: - f.write("\n".join(mfem_file)) + mfem_file.extend( + [ + "MFEM NURBS mesh v1.0", + "", + "# MFEM Geometry Types (see fem/geom.hpp):", + "#", + "# SEGMENT = 1 | SQUARE = 3 | CUBE = 5", + "#", + "# element: 1 ", + "# edge: 0 1 <-- idx increases by one each time", + "# knotvector: [knots]; sizeof(knots) is 1+order+num_ctrl_pts", + "# weights: array of weights corresponding to the NURBS element", + "# FES: list of control points; vertex control points at top, then interior control points", + "", + ] + ) + + 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(["knotvectors", f"{self.elem_cnt}", "\n".join(self.knots), ""]) + + mfem_file.extend(["weights", "\n".join(self.wgts_ends), "\n".join(self.wgts_ints), ""]) + + mfem_file.extend( + [ + "FiniteElementSpace", + "FiniteElementCollection: NURBS", + "VDim: 2", + "Ordering: 1", + "", + "\n".join(self.dof_ends), + "\n".join(self.dof_ints), + "", + ] + ) + + with open(filename, mode="w") as f: + f.write("\n".join(mfem_file)) def parse_args(): - + parser = argparse.ArgumentParser(description="svg2contours: Convert the curves in an SVG to MFEM NURBS mesh") - - parser.add_argument("-i", "--input", - dest="inputfile", - required=True, - type=argparse.FileType('r', encoding='UTF-8'), - help="Input SVG image (*.svg)") - - parser.add_argument("-o", "--output", - dest="outputfile", - default="drawing.mesh", - help="Output file in mfem NURBS mesh format (*.mesh)") - - parser.add_argument("-v", "--verbose", - dest="verbose", - default=False, - action="store_true", - help="verbose output flag") + + parser.add_argument( + "-i", + "--input", + dest="inputfile", + required=True, + type=argparse.FileType("r", encoding="UTF-8"), + help="Input SVG image (*.svg)", + ) + + parser.add_argument( + "-o", + "--output", + dest="outputfile", + default="drawing.mesh", + help="Output file in mfem NURBS mesh format (*.mesh)", + ) + + parser.add_argument( + "-v", "--verbose", dest="verbose", default=False, action="store_true", help="verbose output flag" + ) opts = parser.parse_args() return vars(opts) @@ -307,10 +309,9 @@ def main(): opts = parse_args() verbose = opts["verbose"] - if verbose: print(f"Running from '{os.getcwd()}' with arguments") - for k,v in opts.items(): + for k, v in opts.items(): print(f"\t{k}: {v}") ## Load the SVG document @@ -320,11 +321,11 @@ def main(): if verbose: print(f"Attributes at root of '{input_file}':") - for k,v in doc.tree.getroot().attrib.items(): + for k, v in doc.tree.getroot().attrib.items(): print(f"\t{k}: {v}") coordinate_transform = get_root_transform(doc) - + ## Process SVG paths if verbose: print("SVG paths: \n", paths) @@ -332,9 +333,9 @@ def main(): mfem_data = MFEMData() for p_idx, p in enumerate(paths): - #print(f"""reading {p_idx=} {p=} \n w/ {p.d()=}""") + # print(f"""reading {p_idx=} {p=} \n w/ {p.d()=}""") - is_d_path = 'd' in p.element.keys() + is_d_path = "d" in p.element.keys() attrib = p_idx + 1 reverse_paths = True if np.linalg.det(coordinate_transform) < 0 else False @@ -343,23 +344,23 @@ def main(): continue for seg_idx, seg in enumerate(p): - #print(f"""processing {seg_idx=} {seg=}""") + # print(f"""processing {seg_idx=} {seg=}""") if 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(.5) - - cubic,weights = segment_as_cubic(arc1, reverse_paths) + 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) + + 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) + cubic, weights = segment_as_cubic(seg, reverse_paths) xformed_cubic = transform_cubic(cubic, coordinate_transform) mfem_data.add_cubic_bezier(xformed_cubic, weights, attrib) @@ -368,8 +369,6 @@ def main(): print(f"Wrote '{output_file}' with {mfem_data.vert_cnt} vertices and NURBS {mfem_data.elem_cnt} elements") - if __name__ == "__main__": exitcode = main() sys.exit(exitcode) - From 84dc099b0bab81a2bcd82c7f96ee105ea2eb5cdf Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sat, 27 Jul 2024 21:44:38 -0700 Subject: [PATCH 25/27] Adds option to svg2contours script to reverse orientations --- src/tools/svg2contours/svg2contours.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/tools/svg2contours/svg2contours.py b/src/tools/svg2contours/svg2contours.py index 880119191d..3d03f071f1 100755 --- a/src/tools/svg2contours/svg2contours.py +++ b/src/tools/svg2contours/svg2contours.py @@ -297,6 +297,11 @@ def parse_args(): help="Output file in mfem NURBS mesh format (*.mesh)", ) + parser.add_argument( + "-r", "--reverse", dest="reverse_paths", default=False, action="store_true", + help="reverses paths (can be helpful during coordinate system transformation from y-axis pointing down to up)" + ) + parser.add_argument( "-v", "--verbose", dest="verbose", default=False, action="store_true", help="verbose output flag" ) @@ -339,6 +344,8 @@ def main(): attrib = p_idx + 1 reverse_paths = True if np.linalg.det(coordinate_transform) < 0 else False + if opts["reverse_paths"]: + reverse_paths = not reverse_paths if not all(map(is_path_segment, p)): continue From bb938622423a4c06dc342f084a5e5007c5a40815 Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sat, 27 Jul 2024 23:17:19 -0700 Subject: [PATCH 26/27] Updates data submodule to include SVG images --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 400c4e5066..c0c33018bd 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 400c4e5066620c653b770bd84b77c326da588567 +Subproject commit c0c33018bd6796cfe2ff64e46bb2a16402f00f9c From 15fed7a971ad2da9a7abe98cc7dab51742b4856d Mon Sep 17 00:00:00 2001 From: Kenneth Weiss Date: Sat, 27 Jul 2024 23:22:03 -0700 Subject: [PATCH 27/27] Update RELEASE-NOTES --- RELEASE-NOTES.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index e934f14f6e..9dc4f023e5 100644 --- a/RELEASE-NOTES.md +++ b/RELEASE-NOTES.md @@ -42,6 +42,8 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/ - Improves support for `axom::Array` allocated in unified and pinned memory on GPU platforms. Use of GPU-based operations for Arrays allocated in a unified memory space is controlled with a new method, `Array::setDevicePreference()`. +- Adds `svg2contours` script to convert paths in an SVG file to an MFEM NURBS mesh +- Quest: Adds an example to query winding numbers on an MFEM NURBS mesh ### Changed - `MarchingCubes` masking now uses the mask field's integer values instead of