diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index c917281e1c..d0cf59f466 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -51,6 +51,7 @@ set( primal_headers operators/detail/compute_moments_impl.hpp operators/detail/intersect_bezier_impl.hpp operators/detail/intersect_bounding_box_impl.hpp + operators/detail/intersect_curved_poly_impl.hpp operators/detail/intersect_impl.hpp operators/detail/intersect_ray_impl.hpp diff --git a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp index 66925002a8..4cc84c7930 100644 --- a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp +++ b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp @@ -88,7 +88,7 @@ bool intersect_bezier_curves(const BezierCurve &c1, * As such, the we do not consider the lines to intersect if they do so * at the endpoints \a b or \d, respectively. * - * \note This function assumes the all intersections have multiplicity + * \note This function assumes that all intersections have multiplicity * one, i.e. there are no points at which the curves and their derivatives * both intersect. Thus, the function does not find tangencies. * diff --git a/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp b/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp new file mode 100644 index 0000000000..9bdc394f87 --- /dev/null +++ b/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp @@ -0,0 +1,834 @@ +// Copyright (c) 2017-2022, 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 intersect_curved_poly_impl.hpp + * + * \brief Consists of functions to test intersections among curved polygons + */ + +#ifndef AXOM_PRIMAL_INTERSECTION_CURVED_POLYGON_IMPL_HPP_ +#define AXOM_PRIMAL_INTERSECTION_CURVED_POLYGON_IMPL_HPP_ + +#include "axom/core.hpp" + +#include "axom/primal/geometry/BoundingBox.hpp" +#include "axom/primal/geometry/OrientedBoundingBox.hpp" +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Ray.hpp" +#include "axom/primal/geometry/Segment.hpp" +#include "axom/primal/geometry/Sphere.hpp" +#include "axom/primal/geometry/Triangle.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" + +#include "axom/primal/operators/squared_distance.hpp" +#include "axom/primal/operators/orientation.hpp" +#include "axom/primal/operators/detail/intersect_bezier_impl.hpp" + +#include "axom/fmt.hpp" + +namespace axom +{ +namespace primal +{ +namespace detail +{ +template +bool orient(const BezierCurve& c1, const BezierCurve& c2, T s, T t); + +template +int isContained(const CurvedPolygon& p1, + const CurvedPolygon& p2, + double sq_tol = 1e-10); + +template +class DirectionalWalk +{ +public: + static constexpr int NDIMS = 2; + using CurvedPolygonType = CurvedPolygon; + using IndexArray = std::vector; + +public: + /*! \class EdgeIntersectionInfo + * + * \brief For storing intersection points between edges of \a CurvedPolygon instances so they can be easily sorted by parameter value using std::sort + */ + struct EdgeIntersectionInfo + { + T myTime; // parameter value of intersection on curve on first CurvePolygon + int myEdge; // index of curve on first CurvedPolygon + T otherTime; // parameter value of intersection on curve on second CurvedPolygon + int otherEdge; // index of curve on second CurvedPolygon + int numinter; // unique intersection point identifier + + /// \brief Comparison operator for sorting by parameter value + bool operator<(const EdgeIntersectionInfo& other) const + { + return myTime < other.myTime; + } + }; + + /// Enum to represent different Junction states + enum class JunctionState : int + { + UNINITIALIZED = -1, + NON_JUNCTION, /// Not a junction, e.g. a vertex of the original + APPLIED, /// Use after junction has been applied + CROSS, /// Typical case: edges of polygons intersect + GRAZE /// Atypical case: polygons intersect at a terminal vertex + }; + + /*! + * Helper class for encoding information about the incident edges at an interection + */ + struct Junction + { + using JunctionIndex = int; + using EdgeIndex = int; + static constexpr EdgeIndex INVALID_EDGE_INDEX = -1; + static constexpr int INVALID_JUNCTION_INDEX = -1; + static constexpr int NON_JUNCTION_INDEX = 0; + + bool isActive() const { return junctionState > JunctionState::APPLIED; } + + EdgeIndex currentEdgeIndex(bool active) const + { + return incomingEdgeIndex[active]; + } + EdgeIndex nextEdgeIndex(bool active) const + { + return incomingEdgeIndex[active] + 1; + } + + bool operator==(const Junction& other) const + { + return index == other.index; + } + bool operator!=(const Junction& other) const { return !(*this == other); } + + public: + // indices of polygon edges leading into junction + EdgeIndex incomingEdgeIndex[2] {INVALID_EDGE_INDEX, INVALID_EDGE_INDEX}; + // describes the junction type and status + JunctionState junctionState {JunctionState::UNINITIALIZED}; + JunctionIndex index {INVALID_JUNCTION_INDEX}; + }; + + /*! + * Helper class for indexing an edge of a polygon + * Assumption is that the polygon is not empty + */ + class PolygonEdge + { + public: + using VertexIndex = int; + using EdgeIndex = int; + + PolygonEdge(int polygonID, const IndexArray& endpointJunctionIndices) + : m_polygonID(polygonID) + , m_endpointJunctionIndices(endpointJunctionIndices) + , m_numEdges(endpointJunctionIndices.size()) + { + SLIC_ASSERT_MSG(m_numEdges > 0, "Polygon " << polygonID << " was empty"); + } + + EdgeIndex getIndex() const { return m_index; } + void setIndex(EdgeIndex idx) { m_index = (idx % m_numEdges); } + + int getStartLabel() const { return m_endpointJunctionIndices[prevIndex()]; } + int getEndLabel() const { return m_endpointJunctionIndices[m_index]; } + void advance() { m_index = nextIndex(); } + + bool isJunction() const { return getEndLabel() > 0; } + + bool operator==(const PolygonEdge& other) + { + return m_index == other.m_index && m_polygonID == other.m_polygonID; + } + + bool operator!=(const PolygonEdge& other) { return !(*this == other); } + + private: + EdgeIndex lastIndex() const { return m_numEdges - 1; } + EdgeIndex nextIndex() const + { + return m_index < lastIndex() ? m_index + 1 : 0; + } + EdgeIndex prevIndex() const + { + return m_index > 0 ? m_index - 1 : lastIndex(); + } + + private: + EdgeIndex m_index; + + const int m_polygonID; + const IndexArray& m_endpointJunctionIndices; + const int m_numEdges; + }; + +public: + explicit DirectionalWalk(bool verbose = true) : m_verbose(verbose) { } + + /*! + * Splits edges of each polygon based on the intersections with the other polygon + * When the two polygons intersect, the split polygons are returned in \a psplit + * The edges are labeled by the types of intersections in \a endpointJunctionIndices + * and \a orientation contains the relative orientation of the first intersection + */ + int splitPolygonsAlongIntersections(const CurvedPolygonType& p1, + const CurvedPolygonType& p2, + double sq_tol) + { + // We store intersection information for each edge in EdgeIntersectionInfo structures + using EdgeIntersectionInfoArray = + std::vector>; + EdgeIntersectionInfoArray p1IntersectionData(p1.numEdges()); + EdgeIntersectionInfoArray p2IntersectionData(p2.numEdges()); + EdgeIntersectionInfo firstinter; // Need to do orientation test on first intersection + + // Find all intersections and store + numinters = 0; + for(int i = 0; i < p1.numEdges(); ++i) + { + for(int j = 0; j < p2.numEdges(); ++j) + { + std::vector p1times; + std::vector p2times; + intersect_bezier_curves(p1[i], + p2[j], + p1times, + p2times, + sq_tol, + p1[i].getOrder(), + p2[j].getOrder(), + 0., + 1., + 0., + 1.); + const int edgeIntersections = p1times.size(); + if(edgeIntersections > 0) + { + if(numinters == 0) + { + firstinter = {p1times[0], i, p2times[0], j, 1}; + } + + for(int k = 0; k < edgeIntersections; ++k, ++numinters) + { + p1IntersectionData[i].push_back( + {p1times[k], i, p2times[k], j, numinters + 1}); + p2IntersectionData[j].push_back( + {p2times[k], j, p1times[k], i, numinters + 1}); + + if(m_verbose) + { + SLIC_INFO(fmt::format( + "Found intersection {} -- on edge {} of polygon1 at t={}" + " and edge {} of polygon2 at t={}; intersection point {}", + numinters + 1, + i, + p1times[k], + j, + p2times[k], + p1[i].evaluate(p1times[k]))); + } + } + } + } + } + + if(numinters > 0) + { + // Orient the first intersection point to be sure we get the intersection + orientation = detail::orient(p1[firstinter.myEdge], + p2[firstinter.otherEdge], + firstinter.myTime, + firstinter.otherTime); + + for(int i = 0; i < p1.numEdges(); ++i) + { + std::sort(p1IntersectionData[i].begin(), p1IntersectionData[i].end()); + } + for(int i = 0; i < p2.numEdges(); ++i) + { + std::sort(p2IntersectionData[i].begin(), p2IntersectionData[i].end()); + } + + psplit[0] = p1; + psplit[1] = p2; + + endpointJunctionIndices[0].reserve(p1.numEdges() + numinters); + endpointJunctionIndices[1].reserve(p2.numEdges() + numinters); + + junctions = std::vector(numinters + 1); + + if(m_verbose) + { + SLIC_INFO("Poly1 before split: " << psplit[0]); + SLIC_INFO("Poly2 before split: " << psplit[1]); + } + + this->splitPolygon(0, p1IntersectionData); + this->splitPolygon(1, p2IntersectionData); + + if(m_verbose) + { + SLIC_INFO("Poly1 after split: " << psplit[0]); + SLIC_INFO("Poly2 after split: " << psplit[1]); + } + } + + return numinters; + } + +public: + /*! + * Finds all intersection regions of the two input polygons + * + * This function must be called after splitPolygonsAlongIntersections + * which inserts the intersection points into each polygon and generates + * the \a endpointJunctionIndices structure. + * + * The results will be a vector of polygons inserted into OUT parameter \a pnew + */ + void findIntersectionRegions(std::vector& pnew) + { + PolygonEdge currentEdge[2] = {PolygonEdge(0, endpointJunctionIndices[0]), + PolygonEdge(1, endpointJunctionIndices[1])}; + + if(m_verbose) + { + const int nJunctions = junctions.size(); + SLIC_INFO("At start of 'findIntersectionRegions', there are " + << nJunctions << " junctions:"); + for(int i = 0; i < nJunctions; ++i) + { + SLIC_INFO( + axom::fmt::format("\tJunction {}: incomingEdge[0]: {}, " + "incomingEdge[1]: {}, state: {}, index: {} ", + i, + junctions[i].incomingEdgeIndex[0], + junctions[i].incomingEdgeIndex[1], + junctions[i].junctionState, + junctions[i].index)); + } + } + + // We use junctionIndex to loop through the active junctions starting with index 1 + int junctionIndex = 1; + + // Find all connected regions of the intersection (we're done when we've found all junctions) + const int numJunctions = junctions.size(); + while(junctionIndex < numJunctions) + { + // Attempt to get an initial "active" junction for this polygon + Junction* startJunction = nullptr; + do + { + startJunction = (junctionIndex < numJunctions) + ? &(junctions[junctionIndex++]) + : nullptr; + } while(startJunction != nullptr && !startJunction->isActive()); + + // If we've found all the active junctions, we're done + if(startJunction == nullptr) + { + break; + } + // switch(startJunction->junctionState) + // { + // case JunctionState::CROSS: + // break; + // case JunctionState::GRAZE: + // break; + // default: + // break; + // } + else + { + startJunction->junctionState = JunctionState::APPLIED; + } + + // This variable allows us to switch between the two polygons + bool active = orientation; + + // Set the index of the active edge to that of the start index + PolygonEdge* activeEdge = &(currentEdge[active]); + const auto startEdgeIndex = startJunction->nextEdgeIndex(active); + activeEdge->setIndex(startEdgeIndex); + + if(m_verbose) + { + const int p0Edges = psplit[0].numEdges(); + const int p1Edges = psplit[1].numEdges(); + const int p0CurrIdx = startJunction->currentEdgeIndex(0); + const int p0NextIdx = startJunction->nextEdgeIndex(0) % p0Edges; + const int p1CurrIdx = startJunction->currentEdgeIndex(1); + const int p1NextIdx = startJunction->nextEdgeIndex(1) % p1Edges; + + SLIC_INFO( + axom::fmt::format("Starting with junction {}" + "\n\t edges[0] {} -> {} {{in: {}; out: {} }}" + "\n\t edges[1] {} -> {} {{in: {}; out: {} }}" + "\n\t active: {}", + startJunction->index, + p0CurrIdx, + p0NextIdx, + psplit[0][p0CurrIdx], + psplit[0][p0NextIdx], + p1CurrIdx, + p1NextIdx, + psplit[1][p1CurrIdx], + psplit[1][p1NextIdx], + (active ? 1 : 0))); + } + + CurvedPolygonType aPart; // Tracks the edges of the current intersection polygon + + // Each polygon iterates until it returns to the startJunctions + Junction* junction = nullptr; + while(junction == nullptr || *junction != *startJunction) + { + // Add all edges until we find a junction edge + while(!activeEdge->isJunction()) + { + const auto edgeIndex = activeEdge->getIndex(); + if(m_verbose) + { + SLIC_INFO("Adding edge (non-junction): " << psplit[active][edgeIndex]); + } + aPart.addEdge(psplit[active][edgeIndex]); + activeEdge->advance(); + } + + // Add last leg of previous segment + { + const auto edgeIndex = activeEdge->getIndex(); + if(m_verbose) + { + SLIC_INFO("Adding edge (end of last): " << psplit[active][edgeIndex]); + } + aPart.addEdge(psplit[active][edgeIndex]); + } + + // Handle junction + junction = &(junctions[activeEdge->getEndLabel()]); + SLIC_ASSERT(junction != nullptr); + + if(m_verbose) + { + SLIC_INFO("" + << "Swapped to junction " << junction->index + << "\n\t edges[0] {in: " << junction->currentEdgeIndex(0) + << " -- " << psplit[0][junction->currentEdgeIndex(0)] + << "; out: " << junction->nextEdgeIndex(0) << " -- " + << psplit[0][junction->nextEdgeIndex(0)] << "}" + << "\n\t edges[1] {in: " << junction->currentEdgeIndex(1) + << " -- " << psplit[1][junction->currentEdgeIndex(1)] + << "; out: " << junction->nextEdgeIndex(1) << " -- " + << psplit[1][junction->nextEdgeIndex(1)] << "}" + << "\n\t active: " << (active ? 1 : 0)); + } + + if(junction->isActive()) + { + // swap active edge using junction data + active = !active; + const auto nextEdgeIndex = junction->nextEdgeIndex(active); + activeEdge = &(currentEdge[active]); + activeEdge->setIndex(nextEdgeIndex); + + junction->junctionState = JunctionState::APPLIED; + + if(m_verbose) + { + SLIC_INFO("Swapped to other polygon, edge index: " + << nextEdgeIndex + << "; edge: " << psplit[active][activeEdge->getIndex()]); + } + } + } + // Finalize polygon + pnew.push_back(aPart); + } + } + +private: + /// Splits the polygon with id \a polygonID at every intersection point in \a edgeIntersections. + /// Uses internal array \a edgeLabels to store ids for vertices (junction indices or 0 for original vertices) + void splitPolygon( + int polygonID, + std::vector>& allEdgeIntersections) + { + using axom::utilities::isNearlyEqual; + + CurvedPolygonType& polygon = psplit[polygonID]; + IndexArray& junctionIndices = endpointJunctionIndices[polygonID]; + + bool fixupBeginning = false; + int fixupBeginningJunctionIdx = Junction::INVALID_JUNCTION_INDEX; + + int addedIntersections = 0; + const int numOrigEdges = polygon.numEdges(); + for(int i = 0; i < numOrigEdges; ++i) // foreach edge + { + // mark this edge's endpoint as 'original' + junctionIndices.push_back(Junction::NON_JUNCTION_INDEX); + double previous_tj = 0.; + auto& curEdgeIntersections = allEdgeIntersections[i]; + const int nIntersect = curEdgeIntersections.size(); + for(int j = 0; j < nIntersect; ++j) // foreach intersection on this edge + { + // split edge at parameter t_j + const double t_j = curEdgeIntersections[j].myTime; + const int edgeIndex = i + addedIntersections; + + if(m_verbose) + { + SLIC_INFO( + fmt::format("i {}, j {}, added {}, t_j {}, previous t_j {}, " + "edgeIndex {} (sz: {}), polygon size {}", + i, + j, + addedIntersections, + t_j, + previous_tj, + edgeIndex, + junctionIndices.size(), + polygon.numEdges())); + } + + const bool nearlyZero = isNearlyEqual(t_j, 0.); + + // TODO: Handle case where t_j is nearly 0. + if(!nearlyZero) + { + polygon.splitEdge(edgeIndex, t_j); + + // update edge label + const int idx = curEdgeIntersections[j].numinter; + junctionIndices.insert(junctionIndices.begin() + edgeIndex, idx); + junctions[idx].incomingEdgeIndex[polygonID] = edgeIndex; + //if(junctions[idx].junctionState == JunctionState::UNINITIALIZED) + //{ + junctions[idx].junctionState = JunctionState::CROSS; + //} + junctions[idx].index = idx; + + ++addedIntersections; + previous_tj = t_j; + } + else + { + const int idx = curEdgeIntersections[j].numinter; + + if(edgeIndex > 0) + { + junctionIndices[edgeIndex] = idx; + junctions[idx].incomingEdgeIndex[polygonID] = edgeIndex; + } + else + { + fixupBeginning = true; + fixupBeginningJunctionIdx = idx; + } + + junctions[idx].junctionState = + JunctionState::CROSS; //JunctionState::GRAZE; + junctions[idx].index = idx; + } + + if(m_verbose) + { + SLIC_INFO(axom::fmt::format("junctionIndices: {}\n junctions:", + axom::fmt::join(junctionIndices, " "))); + for(const auto& jj : junctions) + { + SLIC_INFO(fmt::format( + "\t{{index: {}, state: {}, edgeIndex[0]: {}, edgeIndex[1]: {} }}", + jj.index, + jj.junctionState, + jj.incomingEdgeIndex[0], + jj.incomingEdgeIndex[1])); + } + } + + // update remaining intersections on this edge; special case if already at end of curve + if(!isNearlyEqual(1., t_j)) + { + for(int k = j + 1; k < nIntersect; ++k) + { + const double t_k = curEdgeIntersections[k].myTime; + curEdgeIntersections[k].myTime = (t_k - t_j) / (1.0 - t_j); + } + } + else + { + for(int k = j + 1; k < nIntersect; ++k) + { + curEdgeIntersections[k].myTime = 1.; + } + } + } + } + + // If the first vertex of the first edge was a junction, we need to fix it up using the current last edge + if(fixupBeginning) + { + const int edgeIndex = junctionIndices.size() - 1; + junctionIndices[edgeIndex] = fixupBeginningJunctionIdx; + junctions[fixupBeginningJunctionIdx].incomingEdgeIndex[polygonID] = + edgeIndex; + } + } + +public: + // Objects to store completely split polygons (split at every intersection point) and vector with unique id for each intersection and zeros for corners of original polygons. + CurvedPolygonType psplit[2]; // The two completely split polygons will be stored in this array + IndexArray endpointJunctionIndices[2]; // 0 for curves that end in original vertices, unique id for curves that end in intersection points + std::vector junctions; + int numinters {0}; + bool orientation {false}; + + bool m_verbose {false}; +}; + +/*! + * \brief Test whether the regions bounded by CurvedPolygons \a p1 and \a p2 intersect + * \return status true iff \a p1 intersects with \a p2, otherwise false. + * + * \param [in] p1, p2 CurvedPolygon objects to intersect + * \param [in] sq_tol tolerance parameter for the base case of intersect_bezier_curve + * \param [out] pnew vector of type CurvedPolygon holding CurvedPolygon objects representing boundaries of intersection regions. + */ +template +bool intersect_polygon(const CurvedPolygon& p1, + const CurvedPolygon& p2, + std::vector>& pnew, + double sq_tol) +{ + // Intersection is empty if either of the two polygons are empty + if(p1.empty() || p2.empty()) + { + return false; + } + + DirectionalWalk walk; + + int numinters = walk.splitPolygonsAlongIntersections(p1, p2, sq_tol); + + if(numinters > 0) + { + // This performs the directional walking method using the completely split polygons + walk.findIntersectionRegions(pnew); + return true; + } + else // If there are no intersection points, check for containment + { + int containment = isContained(p1, p2, sq_tol); + switch(containment) + { + case 0: + return false; + case 1: + pnew.push_back(p1); + return true; + case 2: + pnew.push_back(p2); + return true; + } + return false; // Catch + } +} + +/*! \brief Checks if two polygons are mutually exclusive or if one includes the other, assuming that they have no intersection points + * + * \param [in] p1, p2 CurvedPolygons to be tested + * \param [in] sq_tol tolerance parameter for the base case of intersect_bezier_curves + * \return 0 if mutually exclusive, 1 if p1 is in p2, 2 if p2 is in p1 + */ +template +int isContained(const CurvedPolygon& p1, + const CurvedPolygon& p2, + double sq_tol) +{ + const int NDIMS = 2; + using BCurve = BezierCurve; + using PointType = typename BCurve::PointType; + + int p1c = 0; + int p2c = 0; + T p1t = .5; + T p2t = .5; + PointType controlPoints[2] = {p1[p1c].evaluate(p1t), p2[p2c].evaluate(p2t)}; + BCurve LineGuess = BCurve(controlPoints, 1); + T line1s = 0.0; + T line2s = 1.0; + for(int j = 0; j < p1.numEdges(); ++j) + { + std::vector temps; + std::vector tempt; + intersect_bezier_curves(LineGuess, + p1[j], + temps, + tempt, + sq_tol, + 1, + p1[j].getOrder(), + 0., + 1., + 0., + 1.); + + for(int i = 0; i < static_cast(temps.size()); ++i) + { + if(temps[i] > line1s) + { + line1s = temps[i]; + p1c = j; + p1t = tempt[i]; + } + } + } + for(int j = 0; j < p2.numEdges(); ++j) + { + std::vector temps; + std::vector tempt; + intersect_bezier_curves(LineGuess, + p2[j], + temps, + tempt, + sq_tol, + 1, + p2[j].getOrder(), + 1., + 0., + 1., + 0.); + intersect(LineGuess, p2[j], temps, tempt); + for(int i = 0; i < static_cast(temps.size()); ++i) + { + if(temps[i] < line2s && temps[i] > line1s) + { + line2s = temps[i]; + p2c = j; + p2t = tempt[i]; + } + } + } + + using Vec3 = primal::Vector; + const bool E1inE2 = + Vec3::cross_product(p1[p1c].dt(p1t), LineGuess.dt(line1s))[2] < 0; + const bool E2inE1 = + Vec3::cross_product(p2[p2c].dt(p2t), LineGuess.dt(line2s))[2] < 0; + if(E1inE2 && E2inE1) + { + return 1; + } + else if(!E1inE2 && !E2inE1) + { + return 2; + } + else + { + return 0; + } +} + +/*! + * \brief Determines orientation of a bezier curve \a c1 with respect to another bezier curve \a c2, + * given that they intersect at parameter values \a s and \a t, respectively + * + * \param [in] c1 the first bezier curve + * \param [in] c2 the second bezier curve + * \param [in] s the parameter value of intersection on \a c1 + * \param [in] t the parameter value of intersection on \a c2 + * \return True if \a c1's positive direction is counterclockwise from \a c2's positive direction + */ +template +bool orient(const BezierCurve& c1, const BezierCurve& c2, T s, T t) +{ + const auto orientation = + primal::Vector::cross_product(c1.dt(s), c2.dt(t))[2]; + return (orientation > 0); +} + +enum class JunctionIntersectionType : int +{ + Non_Intersection = -1, // + Cross, + Type1, + Type2, + Type3, + Type4 +}; + +// when interesection is at a vertex of the original polygon(s), +// determine if the in/out edges of second polygon is on the same +// side of the first polygon +template +JunctionIntersectionType getJunctionIntersectionType( + const BezierCurve& p0In, + const BezierCurve& p0Out, + const BezierCurve& p1In, + const BezierCurve& p1Out) +{ + using SegmentType = primal::Segment; + using VectorType = primal::Vector; + using PointType = primal::Point; + + // TODO: Error checking to ensure the four curves have a common intersection point + // P0In(1) == P0Out(0) == P1In(1) == P1Out(0) + // If not, return JunctionIntersectionType::Non_Intersection + + const VectorType p0Tangents[2] = {-p0In.dt(1.).unitVector(), + p0Out.dt(0.).unitVector()}; + const SegmentType p0Seg(PointType {p0Tangents[0][0], p0Tangents[0][1]}, + PointType {p0Tangents[1][0], p0Tangents[1][1]}); + + const VectorType p1Tangents[2] = {-p1In.dt(1.).unitVector(), + p1Out.dt(0.).unitVector()}; + const SegmentType p1Seg(PointType {p1Tangents[0][0], p1Tangents[0][1]}, + PointType {p1Tangents[1][0], p1Tangents[1][1]}); + + const int orientIn = primal::orientation(p1Seg[0], p0Seg); + const int orientOut = primal::orientation(p1Seg[1], p0Seg); + + if(orientIn == orientOut) + { + const int orient_p1_p0In = primal::orientation(p0Seg[0], p1Seg); + if(orientIn == ON_POSITIVE_SIDE) + { + return orient_p1_p0In == ON_POSITIVE_SIDE + ? JunctionIntersectionType::Type3 + : JunctionIntersectionType::Type1; + } + else + { + { + return orient_p1_p0In == ON_POSITIVE_SIDE + ? JunctionIntersectionType::Type4 + : JunctionIntersectionType::Type2; + } + } + } + else + { + return JunctionIntersectionType::Cross; + } +} + +template +constexpr int DirectionalWalk::Junction::NON_JUNCTION_INDEX; + +} // namespace detail +} // namespace primal +} // namespace axom + +#endif // AXOM_PRIMAL_INTERSECTION_CURVED_POLYGON_IMPL_HPP_ diff --git a/src/axom/primal/operators/intersect.hpp b/src/axom/primal/operators/intersect.hpp index f20eca17d0..2b66bfab9a 100644 --- a/src/axom/primal/operators/intersect.hpp +++ b/src/axom/primal/operators/intersect.hpp @@ -25,11 +25,13 @@ #include "axom/primal/geometry/Sphere.hpp" #include "axom/primal/geometry/Triangle.hpp" #include "axom/primal/geometry/BezierCurve.hpp" +#include "axom/primal/geometry/CurvedPolygon.hpp" #include "axom/primal/operators/detail/intersect_impl.hpp" #include "axom/primal/operators/detail/intersect_ray_impl.hpp" #include "axom/primal/operators/detail/intersect_bounding_box_impl.hpp" #include "axom/primal/operators/detail/intersect_bezier_impl.hpp" +#include "axom/primal/operators/detail/intersect_curved_poly_impl.hpp" namespace axom { @@ -477,6 +479,44 @@ bool intersect(const OrientedBoundingBox& b1, /// \name Bezier Curve Intersection Routines /// @{ +/*! + * \brief Tests if two CurvedPolygon \a p1 and \a p2 intersect. + * \return status true iff \a p1 intersects \a p2, otherwise false. + * + * \param [in] p1 the first CurvedPolygon + * \param [in] p2 the second CurvedPolygon + * \param [out] pnew vector of resulting CurvedPolygons + * \return True if the curvedPolygons intersect, false otherwise. + * + * Finds the set of curved polygons that bound the region of intersection between p1 and p2. + * + * \note This function assumes two dimensional curved polygons in a plane. + * + * \note This function assumes that the curved polygons are in general + * position. Specifically, we assume that all intersections are at points + * and that the component curves don't overlap. + * + * \note This function assumes the all intersections have multiplicity + * one, i.e. there are no points at which the component curves and their + * derivatives both intersect. Thus, the function does not find tangencies. + * + * \note This function assumes that the component curves are half-open, + * i.e. they contain their first endpoint, but not their last endpoint. + * Thus, component curves do not intersect at \f$ s==1 \f$ or at + * \f$ t==1 \f$. + */ +template +bool intersect(const CurvedPolygon& p1, + const CurvedPolygon& p2, + std::vector>& pnew, + double tol = 1E-8) +{ + // for efficiency, linearity check actually uses a squared tolerance + const double sq_tol = tol * tol; + + return detail::intersect_polygon(p1, p2, pnew, sq_tol); +} + /*! * \brief Tests if two Bezier Curves \a c1 and \a c2 intersect. * \return status true iff \a c1 intersects \a c2, otherwise false. diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index 6b11631efd..992874d297 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -15,6 +15,7 @@ set( primal_tests primal_compute_bounding_box.cpp primal_compute_moments.cpp primal_curved_polygon.cpp + primal_curved_polygon_intersect.cpp primal_in_sphere.cpp primal_intersect.cpp primal_intersect_impl.cpp diff --git a/src/axom/primal/tests/primal_curved_polygon.cpp b/src/axom/primal/tests/primal_curved_polygon.cpp index ef0ef95373..00a627c4bf 100644 --- a/src/axom/primal/tests/primal_curved_polygon.cpp +++ b/src/axom/primal/tests/primal_curved_polygon.cpp @@ -12,13 +12,7 @@ #include "axom/config.hpp" #include "axom/slic.hpp" -#include "axom/primal/geometry/Point.hpp" -#include "axom/primal/geometry/Segment.hpp" -#include "axom/primal/geometry/CurvedPolygon.hpp" -#include "axom/primal/geometry/OrientationResult.hpp" -#include "axom/primal/operators/intersect.hpp" -#include "axom/primal/operators/compute_moments.hpp" -#include "axom/primal/operators/orientation.hpp" +#include "axom/primal.hpp" #include #include diff --git a/src/axom/primal/tests/primal_curved_polygon_intersect.cpp b/src/axom/primal/tests/primal_curved_polygon_intersect.cpp new file mode 100644 index 0000000000..9611f63ee4 --- /dev/null +++ b/src/axom/primal/tests/primal_curved_polygon_intersect.cpp @@ -0,0 +1,985 @@ +// Copyright (c) 2017-2022, 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 primal_curved_polygon_intersect.cpp + * \brief This file tests intersections of CurvedPolygon instances + */ + +#include "gtest/gtest.h" + +#include "axom/slic.hpp" +#include "axom/primal.hpp" + +#include "axom/fmt.hpp" + +#include +#include + +namespace primal = axom::primal; + +template +void outputAsSVG( + const std::string& filename, + const primal::CurvedPolygon& polygon1, + const primal::CurvedPolygon& polygon2, + const std::vector> intersectionPolygons) +{ + // Find the bounding box of the set of polygons + primal::BoundingBox bbox; + bbox.addBox(polygon1.boundingBox()); + bbox.addBox(polygon2.boundingBox()); + for(const auto& cp : intersectionPolygons) + { + bbox.addBox(cp.boundingBox()); + } + bbox.scale(1.1); + + std::string header = axom::fmt::format( + "", + bbox.getMin()[0], + bbox.getMin()[1], + bbox.range()[0], + bbox.range()[1]); + std::string footer = ""; + + // lambda to convert a CurvedPolygon to an SVG path string + auto cpToSVG = [](const primal::CurvedPolygon& cp) { + axom::fmt::memory_buffer out; + bool is_first = true; + + for(auto& curve : cp.getEdges()) + { + // Only write out first point for first edge + if(is_first) + { + axom::fmt::format_to(out, "M {} {} ", curve[0][0], curve[0][1]); + is_first = false; + } + + switch(curve.getOrder()) + { + case 1: + axom::fmt::format_to(out, "L {} {} ", curve[1][0], curve[1][1]); + break; + case 2: + axom::fmt::format_to(out, + "Q {} {}, {} {} ", + curve[1][0], + curve[1][1], + curve[2][0], + curve[2][1]); + break; + case 3: + axom::fmt::format_to(out, + "C {} {}, {} {}, {} {} ", + curve[1][0], + curve[1][1], + curve[2][0], + curve[2][1], + curve[3][0], + curve[3][1]); + break; + default: + SLIC_WARNING( + "Unsupported case: can only output up to cubic curves as SVG."); + } + } + return axom::fmt::format(" \n", + axom::fmt::to_string(out)); + }; + + std::string poly1Group; + std::string poly2Group; + std::string intersectionGroup; + + // render polygon1 as SVG + { + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + axom::fmt::format_to(out, cpToSVG(polygon1)); + axom::fmt::format_to(out, " \n"); + poly1Group = axom::fmt::to_string(out); + } + + // render polygon2 as SVG + { + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + axom::fmt::format_to(out, cpToSVG(polygon2)); + axom::fmt::format_to(out, " \n"); + poly2Group = axom::fmt::to_string(out); + } + + //render intersection polygons as SVG + { + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + for(auto& cp : intersectionPolygons) + { + axom::fmt::format_to(out, cpToSVG(cp)); + } + axom::fmt::format_to(out, " \n"); + intersectionGroup = axom::fmt::to_string(out); + } + + // Write the file + { + std::ofstream fs(filename); + fs << header << std::endl; + fs << poly1Group << std::endl; + fs << poly2Group << std::endl; + fs << intersectionGroup << std::endl; + fs << footer << std::endl; + } +} + +/*! + * Helper function to compute the set of intersection polygons given two input polygons + * and to check that they match expectations, stored in \a expbPolygon. + * Intersection polygon is computed to within tolerance \a eps and checks use \a test_eps. + */ +template +void checkIntersection( + const primal::CurvedPolygon& bPolygon1, + const primal::CurvedPolygon& bPolygon2, + const std::vector> expbPolygon, + const double eps = 1e-15, + const double test_eps = 1e-13) +{ + constexpr int DIM = 2; + using CurvedPolygonType = primal::CurvedPolygon; + using BezierCurveType = typename CurvedPolygonType::BezierCurveType; + using PointType = typename BezierCurveType::PointType; + + std::vector intersectionPolys; + + //Compute intersection using algorithm with tolerance of eps + intersect(bPolygon1, bPolygon2, intersectionPolys, eps); + //Check that expected number of intersection regions are found + ASSERT_EQ(expbPolygon.size(), intersectionPolys.size()); + + //Check that expected intersection curves are found to within test_eps + const int nPolygons = expbPolygon.size(); + for(int p = 0; p < nPolygons; ++p) + { + const CurvedPolygonType& polyExp = expbPolygon[p]; + const CurvedPolygonType& polyActual = intersectionPolys[p]; + EXPECT_EQ(polyExp.numEdges(), polyActual.numEdges()); + + const int nEdges = polyExp.numEdges(); + for(int e = 0; e < nEdges; ++e) + { + const BezierCurveType& curveExp = polyExp[e]; + const BezierCurveType& curveActual = polyActual[e]; + EXPECT_EQ(curveExp.getOrder(), curveActual.getOrder()); + + const int nPts = curveExp.getOrder() + 1; + for(int idx = 0; idx < nPts; ++idx) + { + const PointType& ptExp = curveExp[idx]; + const PointType& ptActual = curveActual[idx]; + + for(int d = 0; d < DIM; ++d) + { + EXPECT_NEAR(ptExp[d], ptActual[d], test_eps) + << "Difference in polygon " << p << " edge " << e + << " control point " << idx << " dimension " << d; + } + } + } + } +} + +/*! + * Helper function to create a CurvedPolygon from a list of control points and a list + * of orders of component curves. Control points should be given as a list of Points + * in order of orientation with no duplicates except that the first control point + * should also be the last control point (if the polygon is closed). Orders should + * be given as a list of ints in order of orientation, representing the orders of the component curves. + */ +template +primal::CurvedPolygon createPolygon( + const std::vector> ControlPoints, + const std::vector orders) +{ + using PointType = primal::Point; + using CurvedPolygonType = primal::CurvedPolygon; + using BezierCurveType = primal::BezierCurve; + + const int num_edges = orders.size(); + const int num_unique_control_points = ControlPoints.size(); + + //checks if the orders and control points given will yield a valid curved polygon + { + const int sum_of_orders = std::accumulate(orders.begin(), orders.end(), 0); + EXPECT_EQ(sum_of_orders + 1, num_unique_control_points); + } + + //Converts the control points to BezierCurves of specified orders and stores them in a CurvedPolygon object. + CurvedPolygonType bPolygon; + int iter = 0; + for(int j = 0; j < num_edges; ++j) + { + std::vector subCP; + subCP.assign(ControlPoints.begin() + iter, + ControlPoints.begin() + iter + orders[j] + 1); + BezierCurveType addCurve(subCP, orders[j]); + bPolygon.addEdge(addCurve); + iter += (orders[j]); + } + return bPolygon; +} + +//---------------------------------------------------------------------------------- + +TEST(primal_curvedpolygon, detail_intersection_type) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + using primal::detail::getJunctionIntersectionType; + using primal::detail::JunctionIntersectionType; + + SLIC_INFO("Tests various junction intersections"); + + // the first pair is a horizontal line from left to right + std::vector orders = {1, 1}; + std::vector CP = {PointType {-1, 0}, + PointType {0, 0}, + PointType {1, 0}}; + CurvedPolygonType bPolygon1 = createPolygon(CP, orders); + + // A type2 case + { + std::vector CP2 = {PointType {1, -1}, + PointType {0, 0}, + PointType {-1, -1}}; + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + auto xType = getJunctionIntersectionType(bPolygon1[0], + bPolygon1[1], + bPolygon2[0], + bPolygon2[1]); + EXPECT_EQ(JunctionIntersectionType::Type2, xType); + } +} + +//---------------------------------------------------------------------------------- + +TEST(primal_curvedpolygon, intersection_squares) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of two squares (single region)"); + + std::vector orders = {1, 1, 1, 1}; + + // Unit square scaled by 2 + std::vector CP = {PointType {0, 0}, + PointType {2, 0}, + PointType {2, 2}, + PointType {0, 2}, + PointType {0, 0}}; + + CurvedPolygonType bPolygon1 = createPolygon(CP, orders); + + // Unit square scaled by 2 and offset by (-1,-1) + std::vector CP2 = {PointType {-1, -1}, + PointType {1, -1}, + PointType {1, 1}, + PointType {-1, 1}, + PointType {-1, -1}}; + + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + // Intersection should be a unit square + std::vector expCP = {PointType {1, 0}, + PointType {1, 1}, + PointType {0, 1}, + PointType {0, 0}, + PointType {1, 0}}; + std::vector exporders = {1, 1, 1, 1}; + CurvedPolygonType expbPolygon = createPolygon(expCP, exporders); + + std::vector expbPolygons {expbPolygon}; + checkIntersection(bPolygon1, bPolygon2, expbPolygons); + + // Output intersections as SVG + { + std::vector intersections; + intersect(bPolygon1, bPolygon2, intersections, 1e-15); + outputAsSVG("curved_polygon_intersections_linear_squares.svg", + bPolygon1, + bPolygon2, + intersections); + } +} + +//---------------------------------------------------------------------------------- + +TEST(primal_curvedpolygon, intersection_triangle_linear) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of two linear triangles (single region)"); + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.3, 2.0}, + PointType {0.0, 1.6}, + PointType {0.6, 1.2}}; + std::vector orders = {1, 1, 1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP, orders); + + std::vector CP2 = {PointType {0.71, 1.31}, + PointType {0.41, 2.11}, + PointType {0.11, 1.71}, + PointType {0.71, 1.31}}; + + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + std::vector expCP = { + PointType {0.3091666666666666666666, 1.9755555555555555555}, + PointType {0.11, 1.71}, + PointType {0.5083333333333333333, 1.44444444444444444444}, + PointType {0.3091666666666666666666, 1.9755555555555555555}}; + std::vector exporders = {1, 1, 1}; + CurvedPolygonType expbPolygon = createPolygon(expCP, exporders); + + std::vector expbPolygons {expbPolygon}; + checkIntersection(bPolygon1, bPolygon2, expbPolygons); + + // Output intersections as SVG + { + std::vector intersections; + intersect(bPolygon1, bPolygon2, intersections, 1e-15); + outputAsSVG("curved_polygon_intersections_linear_triangles.svg", + bPolygon1, + bPolygon2, + intersections); + } +} + +//---------------------------------------------------------------------------------- + +TEST(primal_curvedpolygon, intersections_triangle_rectangle) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO( + "Test several intersection cases b/w a linear triangle and rectangle"); + + // Rectangle with bounds, -5 <= x < 5 ; and 0 <= y <= 2 + std::vector rectanglePts = {PointType {-5, 0}, + PointType {5, 0}, + PointType {5, 2}, + PointType {-5, 2}, + PointType {-5, 0}}; + std::vector rectangleOrders = {1, 1, 1, 1}; + + // Equilateral triangle with base from -3 <= x <= 3 and height 1 + CurvedPolygonType bRectangle = createPolygon(rectanglePts, rectangleOrders); + + std::vector triPts = {PointType {-3, -2}, + PointType {3, -2}, + PointType {0, -1}, + PointType {-3, -2}}; + std::vector triOrders = {1, 1, 1}; + + const bool bVerbose = true; + const double EPS = 1e-7; + const double SQ_EPS = EPS * EPS; + + // No intersection case: Triangle is below the rectangle + { + CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); + primal::detail::DirectionalWalk walk(bVerbose); + int nIntersections = + walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); + EXPECT_EQ(0, nIntersections); + + std::vector intersections; + EXPECT_FALSE(primal::intersect(bRectangle, bTriangle, intersections, EPS)); + + outputAsSVG("intersection_test_tri_rect_none.svg", + bRectangle, + bTriangle, + intersections); + } + + // No intersection case: Triangle apex grazes rectangle + { + triPts[2][1] = 0; + CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); + + primal::detail::DirectionalWalk walk(bVerbose); + int nIntersections = + walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); + EXPECT_EQ(1, nIntersections); + + std::vector walkIntersections; + walk.findIntersectionRegions(walkIntersections); + + // EXPECT_EQ(0, walkIntersections.size()); // Warning: Not properly handled yet + + outputAsSVG("intersection_test_tri_rect_lower_graze.svg", + bRectangle, + bTriangle, + walkIntersections); + + std::vector intersections; + // Warning: Not properly handled yet + // EXPECT_FALSE( + primal::intersect(bRectangle, bTriangle, intersections, EPS) + //) + ; + } + + // Simple intersection case: Triangle apex inside rectangle + { + triPts[2][1] = 1; + CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); + + primal::detail::DirectionalWalk walk(bVerbose); + int nIntersections = + walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); + EXPECT_EQ(2, nIntersections); + + std::vector walkIntersections; + walk.findIntersectionRegions(walkIntersections); + EXPECT_EQ(1, walkIntersections.size()); + + outputAsSVG("intersection_test_tri_rect_intersect_tri.svg", + bRectangle, + bTriangle, + walkIntersections); + + std::vector expCP = {PointType {1, 0}, + PointType {0, 1}, + PointType {-1, 0}, + PointType {1, 0}}; + std::vector exporders = {1, 1, 1}; + + std::vector expbPolygons = { + createPolygon(expCP, exporders)}; + + checkIntersection(bRectangle, bTriangle, expbPolygons); + } + + // Grazing intersection case: Triangle apex intersects top + { + triPts[2][1] = 2; + CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); + + primal::detail::DirectionalWalk walk(bVerbose); + int nIntersections = + walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); + EXPECT_EQ(3, nIntersections); + + std::vector walkIntersections; + walk.findIntersectionRegions(walkIntersections); + EXPECT_EQ(1, walkIntersections.size()); + + //EXPECT_EQ(3, walkIntersections[0].numEdges()); // Warning: Not properly handled yet + + outputAsSVG("intersection_test_tri_rect_upper_graze.svg", + bRectangle, + bTriangle, + walkIntersections); + + std::vector intersections; + // Warning: Not properly handled yet + EXPECT_TRUE(primal::intersect(bRectangle, bTriangle, intersections, EPS)); + } + + // intersection case: Triangle apex above rectangle + { + triPts[2][1] = 3; + CurvedPolygonType bTriangle = createPolygon(triPts, triOrders); + + primal::detail::DirectionalWalk walk(bVerbose); + int nIntersections = + walk.splitPolygonsAlongIntersections(bRectangle, bTriangle, SQ_EPS); + EXPECT_EQ(4, nIntersections); + + std::vector walkIntersections; + walk.findIntersectionRegions(walkIntersections); + EXPECT_EQ(1, walkIntersections.size()); + + outputAsSVG("intersection_test_tri_rect_intersect_rect.svg", + bRectangle, + bTriangle, + walkIntersections); + + std::vector expCP = {PointType {1.8, 0}, + PointType {0.6, 2}, + PointType {-0.6, 2}, + PointType {-1.8, 0}, + PointType {1.8, 0}}; + std::vector exporders = {1, 1, 1, 1}; + + std::vector expbPolygons = { + createPolygon(expCP, exporders)}; + + checkIntersection(bRectangle, bTriangle, expbPolygons); + } +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, intersection_triangle_quadratic) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersecting two quadratic triangles (single region)"); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + std::vector CP = {PointType {0.6, 1.2}, + PointType {0.4, 1.3}, + PointType {0.3, 2.0}, + PointType {0.27, 1.5}, + PointType {0.0, 1.6}, + PointType {0.1, 1.5}, + PointType {0.6, 1.2}}; + std::vector orders = {2, 2, 2}; + CurvedPolygonType bPolygon1 = createPolygon(CP, orders); + + std::vector CP2 = {PointType {0.71, 1.31}, + PointType {0.51, 1.41}, + PointType {0.41, 2.11}, + PointType {0.38, 1.61}, + PointType {0.11, 1.71}, + PointType {0.21, 1.61}, + PointType {0.71, 1.31}}; + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + std::vector expCP = { + PointType {0.335956890729522, 1.784126953773395}, + PointType {0.297344765794753, 1.718171485335525}, + PointType {0.2395677533016981, 1.700128235793371}, + PointType {0.221884203146682, 1.662410644580941}, + PointType {0.199328465398189, 1.636873522352205}, + PointType {0.277429214338182, 1.579562422716502}, + PointType {0.408882616650578, 1.495574996394597}, + PointType {0.368520120719339, 1.616453177259694}, + PointType {0.335956890729522, 1.784126953773394}}; + std::vector exporders = {2, 2, 2, 2}; + CurvedPolygonType expbPolygon = createPolygon(expCP, exporders); + std::vector expbPolygons = {expbPolygon}; + + checkIntersection(bPolygon1, bPolygon2, expbPolygons, 1e-15, 1e-13); +} + +//---------------------------------------------------------------------------------- +TEST(primal_curvedpolygon, intersection_triangle_quadratic_two_regions) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of two quadratic triangles (two regions)"); + + std::vector CP1 = {PointType {0.6, 1.2}, + PointType {0.4, 1.3}, + PointType {0.3, 2.0}, + PointType {0.27, 1.5}, + PointType {0.0, 1.6}, + PointType {0.1, 1.5}, + PointType {0.6, 1.2}}; + + std::vector CP2 = {PointType {1.0205, 1.6699}, + PointType {0.8339, 1.5467}, + PointType {0.1777, 1.8101}, + PointType {0.5957, 1.5341}, + PointType {0.3741, 1.3503}, + PointType {0.5107, 1.3869}, + PointType {1.0205, 1.6699}}; + std::vector orders = {2, 2, 2}; + + std::vector expCP1 = { + PointType {0.343364196589264, 1.747080669655736}, + PointType {0.305984025190458, 1.760433098612141}, + PointType {0.266743999290327, 1.775316659915674}, + PointType {0.263419346128088, 1.763343410502168}, + PointType {0.259796003065908, 1.752116885838515}, + PointType {0.320641367919239, 1.705796408318085}, + PointType {0.362111919147859, 1.662268860466508}, + PointType {0.352450139541348, 1.702947255097842}, + PointType {0.343364196589264, 1.747080669655736}}; + + std::vector expCP2 = { + PointType {0.435276730907216, 1.423589798138227}, + PointType {0.416268597450954, 1.385275578571685}, + PointType {0.374100000000000, 1.350300000000000}, + PointType {0.404839872482010, 1.358536305511285}, + PointType {0.454478985809487, 1.379250566393211}, + PointType {0.444689566319939, 1.400290430035245}, + PointType {0.435276730907216, 1.423589798138227}}; + + std::vector exporder1 = {2, 2, 2, 2}; + std::vector exporder2 = {2, 2, 2}; + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + CurvedPolygonType expbPolygon1 = createPolygon(expCP1, exporder1); + CurvedPolygonType expbPolygon2 = createPolygon(expCP2, exporder2); + std::vector expIntersections = {expbPolygon1, expbPolygon2}; + + checkIntersection(bPolygon1, bPolygon2, expIntersections); + + // Output intersections as SVG + { + std::vector intersections; + intersect(bPolygon1, bPolygon2, intersections, 1e-15); + outputAsSVG( + "curved_polygon_intersections_quadratic_triangles_two_regions.svg", + bPolygon1, + bPolygon2, + intersections); + } +} + +TEST(primal_curvedpolygon, area_intersection_triangle_inclusion) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of two quadratic triangles (inclusion)"); + + std::vector CP1 = {PointType {0.0, 0.0}, + PointType {0.5, 0.0}, + PointType {1.0, 0.0}, + PointType {0.5, 0.5}, + PointType {0.0, 1.0}, + PointType {0.0, 0.5}, + PointType {0.0, 0.0}}; + + std::vector CP2 = {PointType {0.05, 0.05}, + PointType {0.30, 0.05}, + PointType {0.55, 0.05}, + PointType {0.30, 0.30}, + PointType {0.05, 0.55}, + PointType {0.05, 0.30}, + PointType {0.05, 0.05}}; + std::vector orders = {2, 2, 2}; + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + std::vector expIntersections = {bPolygon2}; + + checkIntersection(bPolygon1, bPolygon2, expIntersections); + checkIntersection(bPolygon2, bPolygon1, expIntersections); +} + +TEST(primal_curvedpolygon, doubleIntersection) +{ + const double EPS = 1e-8; + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Tests multiple intersections along an edge"); + + { + // Unit square + std::vector CP1 = {PointType {0, 0}, + PointType {1, 0}, + PointType {1, 1}, + PointType {0, 1}, + PointType {0, 0}}; + std::vector orders1 = {1, 1, 1, 1}; + + // Bi-gon defined by a quadratic edge and a straight line + std::vector CP2 = {PointType {0.8, .25}, + PointType {2.0, .50}, + PointType {0.8, .75}, + PointType {0.8, .25}}; + std::vector orders2 = {2, 1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders1); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders2); + + // Find intersections using DirectionalWalk class + double walkIntersectionArea = 0.; + { + const bool verbose = true; + primal::detail::DirectionalWalk walk(verbose); + int numIntersections = + walk.splitPolygonsAlongIntersections(bPolygon1, bPolygon2, EPS * EPS); + + EXPECT_EQ(2, numIntersections); + EXPECT_NEAR(primal::area(bPolygon1), primal::area(walk.psplit[0]), EPS); + EXPECT_NEAR(primal::area(bPolygon2), primal::area(walk.psplit[1]), EPS); + + SLIC_INFO("Checking for intersections between " << bPolygon1 << " and " + << bPolygon2); + + std::vector regions; + walk.findIntersectionRegions(regions); + EXPECT_EQ(1, regions.size()); + + if(!regions.empty()) + { + EXPECT_EQ(4, regions[0].numEdges()); + walkIntersectionArea = primal::area(regions[0]); + } + else + { + FAIL() << "Expected the two polygons to intersect"; + } + + SLIC_INFO("Found intersections (directional walk): "); + for(auto cp : regions) + { + SLIC_INFO("\t" << cp); + } + } + + // Find interesections using primal::intersect + double directIntersectionArea = 0.; + { + std::vector regions; + bool intersects = intersect(bPolygon1, bPolygon2, regions, EPS); + EXPECT_TRUE(intersects); + EXPECT_EQ(1, regions.size()); + + if(!regions.empty()) + { + EXPECT_EQ(4, regions[0].numEdges()); + directIntersectionArea = primal::area(regions[0]); + } + else + { + FAIL() << "Expected the two polygons to intersect"; + } + + SLIC_INFO("Found intersections (direct): "); + for(auto cp : regions) + { + SLIC_INFO("\t" << cp); + } + } + + EXPECT_NEAR(walkIntersectionArea, directIntersectionArea, EPS); + } +} + +TEST(primal_curvedpolygon, regression) +{ + const double EPS = 1e-8; + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of pairs of polygons from regression data"); + + // First test: Intersecting a pair of linear quadrilaterals + // Note: One of the intersections happens betwen a vertex and edge + // and is not yet properly handled by primal::intersect + { + std::vector CP1 = {PointType {1, -2}, + PointType {0, -1}, + PointType {0, 1}, + PointType {1, 2}, + PointType {1, -2}}; + + std::vector CP2 = {PointType {-.9, -2}, + PointType {0.1, -1}, + PointType {2.1, -1}, + PointType {3.1, -2}, + PointType {-.9, -2}}; + std::vector orders = {1, 1, 1, 1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + std::vector expIntersections; + + intersect(bPolygon1, bPolygon2, expIntersections, EPS); + + SLIC_INFO("There were " << expIntersections.size() + << " intersection polygons"); + for(auto cp : expIntersections) + { + SLIC_INFO("\t" << cp); + } + } +} + +TEST(primal_curvedpolygon, adjacent_squares) +{ + const double EPS = 1e-8; + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon; + using PointType = primal::Point; + + SLIC_INFO("Test intersection of pairs of adjacent squares"); + + // Sharing an edge + { + std::vector CP1 = {PointType {-1, 0}, + PointType {0, 0}, + PointType {0, 1}, + PointType {-1, 1}, + PointType {-1, 0}}; + + std::vector CP2 = {PointType {0, 0}, + PointType {0, 1}, + PointType {1, 1}, + PointType {1, 0}, + PointType {0, 0}}; + std::vector orders = {1, 1, 1, 1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + { + std::vector expIntersections; + intersect(bPolygon1, bPolygon2, expIntersections, EPS); + + EXPECT_TRUE(expIntersections.empty()); + + if(!expIntersections.empty()) + { + SLIC_INFO("There were " << expIntersections.size() + << " intersection polygons"); + for(auto cp : expIntersections) + { + SLIC_INFO("\t" << cp); + } + outputAsSVG("adjacent_squares_edge_LR.svg", + bPolygon1, + bPolygon2, + expIntersections); + } + } + + { + std::vector expIntersections; + intersect(bPolygon2, bPolygon1, expIntersections, EPS); + + EXPECT_TRUE(expIntersections.empty()); + + if(!expIntersections.empty()) + { + SLIC_INFO("There were " << expIntersections.size() + << " intersection polygons"); + for(auto cp : expIntersections) + { + SLIC_INFO("\t" << cp); + } + outputAsSVG("adjacent_squares_edge_RL.svg", + bPolygon2, + bPolygon1, + expIntersections); + } + } + } + + // Sharing a vertex + if(false) + { + std::vector CP1 = {PointType {-1, 0}, + PointType {0, 0}, + PointType {0, 1}, + PointType {-1, 1}, + PointType {-1, 0}}; + + std::vector CP2 = {PointType {0, 1}, + PointType {0, 2}, + PointType {1, 2}, + PointType {1, 1}, + PointType {0, 1}}; + std::vector orders = {1, 1, 1, 1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP1, orders); + CurvedPolygonType bPolygon2 = createPolygon(CP2, orders); + + { + std::vector expIntersections; + intersect(bPolygon1, bPolygon2, expIntersections, EPS); + + EXPECT_TRUE(expIntersections.empty()); + + if(!expIntersections.empty()) + { + SLIC_INFO("There were " << expIntersections.size() + << " intersection polygons"); + for(auto cp : expIntersections) + { + SLIC_INFO("\t" << cp); + } + outputAsSVG("adjacent_squares_corner_LR.svg", + bPolygon1, + bPolygon2, + expIntersections); + } + } + + if(false) + { + std::vector expIntersections; + intersect(bPolygon2, bPolygon1, expIntersections, EPS); + + EXPECT_TRUE(expIntersections.empty()); + + if(!expIntersections.empty()) + { + SLIC_INFO("There were " << expIntersections.size() + << " intersection polygons"); + for(auto cp : expIntersections) + { + SLIC_INFO("\t" << cp); + } + outputAsSVG("adjacent_squares_corner_RL.svg", + bPolygon2, + bPolygon1, + expIntersections); + } + } + } +} + +//---------------------------------------------------------------------------------- +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::SimpleLogger logger; + + result = RUN_ALL_TESTS(); + + return result; +} diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 353ab529af..67c855fed3 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -209,3 +209,23 @@ if (ENABLE_FORTRAN) endif() endif() endif() + +if(MFEM_FOUND) + blt_add_executable( + NAME quest_high_order_remap_ex + SOURCES quest_high_order_remap.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON ${quest_example_depends} mfem + FOLDER axom/quest/examples + ) + + #string(STRIP "${AXOM_DISABLE_UNUSED_PARAMETER_WARNINGS}" MFEM_COMPILE_FLAGS) + + #if (ENABLE_CUDA) + # set(MFEM_COMPILE_FLAGS "-Xcompiler=${MFEM_COMPILE_FLAGS}") + #endif() + + #blt_add_target_compile_flags( TO quest_high_order_remap_ex FLAGS "${MFEM_COMPILE_FLAGS}" ) + +endif() + diff --git a/src/axom/quest/examples/quest_high_order_remap.cpp b/src/axom/quest/examples/quest_high_order_remap.cpp new file mode 100644 index 0000000000..462ed05cb1 --- /dev/null +++ b/src/axom/quest/examples/quest_high_order_remap.cpp @@ -0,0 +1,877 @@ +// Copyright (c) 2017-2022, 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_high_order_remap.cpp + * \brief Demonstrates conservative field remap on 2D high order meshes + */ + +// Axom includes +#include "axom/config.hpp" +#include "axom/core.hpp" +#include "axom/slic.hpp" +#include "axom/primal.hpp" +#include "axom/spin.hpp" + +#ifdef AXOM_USE_MFEM + #include "mfem.hpp" +#else + #error "This example requires mfem" +#endif + +#include "axom/fmt.hpp" +#include "axom/CLI11.hpp" + +#include + +namespace primal = axom::primal; +namespace spin = axom::spin; + +/** + * \brief Wrapper for a 2D mfem mesh + * + * Helps with conversion to BezierCurves and CurvedPolygons + * + * The meshwrapper assumes ownership of the wrapped mesh + */ +class MeshWrapper +{ +public: + constexpr static int DIM = 2; + using BBox = primal::BoundingBox; + using SpacePoint = primal::Point; + + using CurvedPolygonType = primal::CurvedPolygon; + using BCurve = CurvedPolygonType::BezierCurveType; + +private: + /// \brief Checks if the mesh's nodes are in the Bernstein basis + bool isBernsteinBasis() const + { + auto* fec = m_mesh->GetNodalFESpace()->FEColl(); + + if(fec == nullptr) + { + return false; + } + + if(const mfem::H1_FECollection* h1Fec = + dynamic_cast(fec)) + { + return h1Fec->GetBasisType() == mfem::BasisType::Positive; + } + + if(const mfem::L2_FECollection* l2Fec = + dynamic_cast(fec)) + { + return l2Fec->GetBasisType() == mfem::BasisType::Positive; + } + + if(dynamic_cast(fec) || + dynamic_cast(fec) || + dynamic_cast(fec)) + { + return true; + } + + return false; + } + +public: + MeshWrapper() : m_mesh(nullptr) { } + + MeshWrapper(mfem::Mesh* mesh) { setMesh(mesh); } + + ~MeshWrapper() + { + if(m_mesh != nullptr) + { + delete m_mesh; + m_mesh = nullptr; + } + } + + /// Sets the mfem mesh pointer for this MeshWrapper instance + void setMesh(mfem::Mesh* mesh) + { + SLIC_ERROR_IF(mesh == nullptr, "Mesh was null"); + m_mesh = mesh; + + // Check that mesh is high order + SLIC_ERROR_IF(m_mesh->GetNodes() == nullptr, "The mesh must be high order."); + + // Check that mesh nodes are using Bernstein basis + SLIC_ERROR_IF(!isBernsteinBasis(), + "The mesh must be in the Bernstein basis"); + + const double tol = 1E-8; + computeBoundingBoxes(1 + tol); + } + + int numVertices() const { return m_mesh->GetNV(); } + int numElements() const { return m_mesh->GetNE(); } + + bool hasMesh() const { return m_mesh != nullptr; } + mfem::Mesh* getMesh() { return m_mesh; } + const mfem::Mesh* getMesh() const { return m_mesh; } + + const BBox& elementBoundingBox(int i) const { return m_boundingBoxes[i]; } + const BBox& meshBoundingBox() const { return m_meshBBox; } + + /*! + * \brief Transform the mfem element into a primal::CurvedPolygon + * + * \param elemId The index of the element + * \return The element as a CurvedPolygon composed of BezierCurves + */ + CurvedPolygonType elemAsCurvedPolygon(int elemId) + { + SLIC_ASSERT(elemId >= 0 && elemId < numElements()); + + auto* fes = m_mesh->GetNodalFESpace(); + auto* nodes = m_mesh->GetNodes(); + + // Get the edge Ids for this element + mfem::Array edgeIds, edgeOrients; + m_mesh->GetElementEdges(elemId, edgeIds, edgeOrients); + const int nEdges = edgeIds.Size(); + + CurvedPolygonType poly(nEdges); + const int order = fes->GetOrder(0); + + // Get the grid function data associated with this edge + mfem::Array dofIndices; + BCurve curve(order); + for(int e = 0; e < nEdges; ++e) + { + // get the dof (degree of freedom) indices for this edge + fes->GetEdgeDofs(edgeIds[e], dofIndices); + + // possibly reverse the dofs, based on the orientation + // Note: The dofs are ordered by vertices, then by edge + const bool bReverse = (edgeOrients[e] > 0); + if(bReverse) + { + curve[0] = spacePointFromDof(dofIndices[1], fes, nodes); + for(int p = 1; p < order; ++p) + { + curve[p] = spacePointFromDof(dofIndices[order - (p - 1)], fes, nodes); + } + curve[order] = spacePointFromDof(dofIndices[0], fes, nodes); + } + else + { + curve[0] = spacePointFromDof(dofIndices[0], fes, nodes); + for(int p = 1; p < order; ++p) + { + curve[p] = spacePointFromDof(dofIndices[p + 1], fes, nodes); + } + curve[order] = spacePointFromDof(dofIndices[1], fes, nodes); + } + + // Note: mfem's orientation is reversed w.r.t. primal's + //SLIC_INFO("Elem " << elemId << " edge " << e << " -- curve: " << curve); + //curve.reverseOrientation(); + poly[nEdges - e - 1] = curve; + } + // std::cout << poly << std::endl; + return poly; + } + +private: + /// Get the coordinates of the point from the dof index + SpacePoint spacePointFromDof(int idx, + const mfem::FiniteElementSpace* fes, + const mfem::GridFunction* nodes) + { + return SpacePoint {(*nodes)(fes->DofToVDof(idx, 0)), + (*nodes)(fes->DofToVDof(idx, 1))}; + } + + /*! + * \brief Compute the element and mesh bounding boxes + * + * \param[in] bboxScaleFac Scale factor to increase the bounding boxes + * \pre bboxScaleFac >= 1. + */ + void computeBoundingBoxes(double bboxScaleFac) + { + SLIC_ASSERT(bboxScaleFac >= 1.); + + m_boundingBoxes.resize(numElements()); + m_meshBBox.clear(); + + auto* nodes = m_mesh->GetNodes(); + auto* fes = m_mesh->GetNodalFESpace(); + mfem::Array dofIndices; + + const int NE = numElements(); + for(int elem = 0; elem < NE; ++elem) + { + auto& bbox = m_boundingBoxes[elem]; + bbox.clear(); + + // Add each dof of the element to the bbox + // Note: positivity of Bernstein bases ensures that convex + // hull of element nodes contain entire element + fes->GetElementDofs(elem, dofIndices); + for(int i = 0; i < dofIndices.Size(); ++i) + { + int nIdx = dofIndices[i]; + bbox.addPoint(spacePointFromDof(nIdx, fes, nodes)); + } + + // Slightly scale the bbox to account for numerical noise + bbox.scale(bboxScaleFac); + + m_meshBBox.addBox(bbox); + } + } + +private: + mfem::Mesh* m_mesh; + + std::vector m_boundingBoxes; + BBox m_meshBBox; +}; + +struct Remapper +{ +public: + constexpr static int DIM = 2; + using CandidateList = std::vector; + +private: + using GridType = spin::ImplicitGrid; + +public: + Remapper() = default; + + ~Remapper() + { + fecMap.DeleteData(true); + fesMap.DeleteData(true); + } + + mfem::GridFunction* project_to_pos_basis(const mfem::GridFunction* gf, + bool& is_new) + { + mfem::GridFunction* out_pos_gf = nullptr; + is_new = false; + + SLIC_ASSERT(gf != nullptr); + + const mfem::FiniteElementSpace* nodal_fe_space = gf->FESpace(); + SLIC_ERROR_IF(nodal_fe_space == nullptr, + "project_to_pos_basis(): nodal_fe_space is NULL!"); + + const mfem::FiniteElementCollection* nodal_fe_coll = nodal_fe_space->FEColl(); + SLIC_ERROR_IF(nodal_fe_coll == nullptr, + "project_to_pos_basis(): nodal_fe_coll is NULL!"); + + mfem::Mesh* gf_mesh = nodal_fe_space->GetMesh(); + SLIC_ERROR_IF(gf_mesh == nullptr, + "project_to_pos_basis(): gf_mesh is NULL!"); + + int order = nodal_fe_space->GetOrder(0); + int dim = gf_mesh->Dimension(); + + auto* pos_fe_coll = + new mfem::H1_FECollection(order, dim, mfem::BasisType::Positive); + + // { + // mfem::Geometry::Type geom_type = gf_mesh->GetElementBaseGeometry(0); + // int map_type = (nodal_fe_coll != nullptr) + // ? nodal_fe_coll->FiniteElementForGeometry(geom_type)->GetMapType() + // : static_cast(mfem::FiniteElement::VALUE); + + // fecMap.Register("src_fec", fec, true); + // mfem::FiniteElementCollection pos_fe_coll = *nodal_fe_coll; + // detail::get_pos_fec(nodal_fe_coll, order, dim, map_type); + // } + + SLIC_ERROR_IF(pos_fe_coll == nullptr, + "Problem generating a positive finite element collection " + << "corresponding to the mesh's '" << nodal_fe_coll->Name() + << "' finite element collection."); + + //SLIC_INFO("Good so far... pos_fe_coll is not null. Making FESpace and GridFunction."); + const int dims = nodal_fe_space->GetVDim(); + + // Create a positive (Bernstein) grid function for the nodes + mfem::FiniteElementSpace* pos_fe_space = + new mfem::FiniteElementSpace(gf_mesh, pos_fe_coll, dims); + mfem::GridFunction* pos_nodes = new mfem::GridFunction(pos_fe_space); + + // m_pos_nodes takes ownership of pos_fe_coll's memory (and pos_fe_space's memory) + pos_nodes->MakeOwner(pos_fe_coll); + + // Project the nodal grid function onto this + pos_nodes->ProjectGridFunction(*gf); + + out_pos_gf = pos_nodes; + is_new = true; + + SLIC_WARNING_IF( + out_pos_gf == nullptr, + "project_to_pos_basis(): Construction failed; out_pos_gf is NULL!"); + + return out_pos_gf; + } + + /*! + * \brief Loads a mesh (source or target) and applies some simple transformations + * + * \param isSource Determines is we're loading the source mesh (true) or target mesh (false) + * \param fname File pointing to an mfem mesh. If empty, we'll generate a Cartesian mesh over the unit square + * \param offset_x Offset for translating the mesh in the x direction + * \param offset_y Offset for translating the mesh in the y direction + * \param scale Factor for uniformly scaling the mesh + * \param mref Number of uniform refinements to apply to the mesh + * \param order Polynomial order for the mesh's nodal grid function + */ + void loadMesh(bool isSource, + const std::string& fname, + double offset_x, + double offset_y, + double scale, + int mref, + int order) + { + mfem::Mesh* mesh = nullptr; + + // Load mesh from file or as Cartesian mesh + if(!fname.empty()) + { + mesh = new mfem::Mesh(fname.c_str(), 1, 1); + } + else + { + const auto quadType = mfem::Element::QUADRILATERAL; + const bool generateEdges = true; + const int res = 1; + mesh = new mfem::Mesh(res, res, quadType, generateEdges); + } + + // Apply uniform refinement + for(int i = 0; i < mref; ++i) + { + mesh->UniformRefinement(); + } + + // Ensure that mesh has high order nodes + mesh->SetCurvature(order); + + // Scale and offset mesh + xformMesh(mesh, scale, offset_x, offset_y); + + // dump original mesh + { + std::ofstream file; + file.open(axom::fmt::format("{}_mesh_orig.mfem", isSource ? "src" : "tgt")); + mesh->Print(file); + } + + // project nodes to Berstein basis, if necessary + { + bool is_mesh_gf_new {false}; + mfem::GridFunction* mesh_nodes = mesh->GetNodes(); + mfem::GridFunction* pos_mesh_nodes_ptr = + project_to_pos_basis(mesh_nodes, is_mesh_gf_new); + + mfem::GridFunction& pos_mesh_nodes = + (is_mesh_gf_new ? *pos_mesh_nodes_ptr : *mesh_nodes); + mesh->NewNodes(pos_mesh_nodes, true); + } + + // dump modified mesh + { + std::ofstream file; + file.open(axom::fmt::format("{}_mesh_set.mfem", isSource ? "src" : "tgt")); + mesh->Print(file); + } + + // set as active source/target mesh + if(isSource) + { + srcMesh.setMesh(mesh); + } + else + { + tgtMesh.setMesh(mesh); + } + + SLIC_INFO(axom::fmt::format( + "Loaded {} mesh from {} w/ {} elements." + "\n\t(Slightly inflated) mesh bounding box: {}", + isSource ? "source" : "target", + fname.empty() ? "Cartesian grid" : fname, + mesh->GetNE(), + isSource ? srcMesh.meshBoundingBox() : tgtMesh.meshBoundingBox())); + } + + /// Setup the implicit grid spatial index over the source mesh + void initSpatialIndex() + { + const int NE = srcMesh.numElements(); + grid.initialize(srcMesh.meshBoundingBox(), nullptr, NE); + + for(int i = 0; i < NE; ++i) + { + grid.insert(srcMesh.elementBoundingBox(i), i); + } + } + + /*! + * Computes the overlap areas between the elements of the target mesh + * to the elements of the source mesh + */ + double computeOverlapAreas() + { + double totalArea = 0.0; + double correctArea = 0.0; + const int nTargetElems = tgtMesh.numElements(); + + double calcE = 0.0; + for(int i = 0; i < nTargetElems; ++i) + { + // Finds the candidates from the source mesh that can intersect this target element + auto candidates = getSourceCandidates(i); + if(candidates.empty()) + { + continue; + } + + auto tgtPoly = tgtMesh.elemAsCurvedPolygon(i); + //SLIC_INFO("Target Element: \n\t" << tgtPoly); + correctArea += primal::area(tgtPoly); + //SLIC_INFO("Target elem " << i << " -- area " << primal::area(tgtPoly) + // << " -- bbox " << tgtMesh.elementBoundingBox(i) + //); + + double A = 0.0; + for(int srcElem : candidates) + { + auto srcPoly = srcMesh.elemAsCurvedPolygon(srcElem); + //SLIC_INFO("\tSource Element: " << srcPoly); + //SLIC_INFO( + // "\tSource elem \n\t\t" << srcElem << "\n\t\t -- area " << primal::area(srcPoly) + // << " -- bbox " << srcMesh.elementBoundingBox(srcElem) + //); + + std::vector> pnew; + tgtPoly.reverseOrientation(); + srcPoly.reverseOrientation(); + if(primal::intersect(tgtPoly, srcPoly, pnew, 1e-8)) + { + for(int i = 0; i < static_cast(pnew.size()); ++i) + { + A -= primal::area(pnew[i]); + // SLIC_INFO("** Intersection area :" << -primal::area(pnew[i]) + // ); + } + } + srcPoly.reverseOrientation(); + tgtPoly.reverseOrientation(); + // SLIC_INFO("* Calculated area: " << srcElem + // << " -- area " << A + // //<< " -- bbox " << srcMesh.elementBoundingBox(srcElem) + // ); + } + calcE += abs(primal::area(tgtPoly) - A); + totalArea += A; + // SLIC_INFO("Calculated Area :" << A); + } + // std::cout << inclusion << ", "; + // std::cout << calcE << ", "; + // const double tgt_scale = .999712378102150; + // const double tgt_trans = .0001345747181586; + // const double tgt_scale = .9999712378102150; + // const double tgt_trans = .00001345747181586; + // const double tgt_scale = .99999999999712378102150; + // const double tgt_trans = .000000000001345747181586; + // double trueError = (tgt_scale*tgt_scale-totalArea); + // std::cout << trueError << ", "; + std::cout << "Calculated area (supermesh): " << std::fixed << totalArea + << std::endl; + std::cout << "Calculated area (target mesh): " << correctArea << std::endl; + return (totalArea - correctArea); + } + + /*! + * Gets the IDs of the candidate elements from the source mesh + * that might intersect with \a targetElemID from the target mesh + */ + CandidateList getSourceCandidates(int targetElemID) const + { + SLIC_ASSERT(targetElemID < tgtMesh.numElements()); + using BitsetType = GridType::BitsetType; + + CandidateList filteredCandidates; + + auto& targetBBox = tgtMesh.elementBoundingBox(targetElemID); + auto candidateBits = grid.getCandidates(targetBBox); + + // Filter the candidates; return as std vec + for(auto idx = candidateBits.find_first(); idx != BitsetType::npos; + idx = candidateBits.find_next(idx)) + { + if(primal::intersect(targetBBox, srcMesh.elementBoundingBox(idx))) + filteredCandidates.push_back(idx); + } + + return filteredCandidates; + } + + void outputAsSVG() + { + std::string header = R"html()html"; + std::string footer = ""; + + // lambda to convert a CurvedPolygon to an SVG path string + auto cpToSVG = [](const MeshWrapper::CurvedPolygonType& cp) { + axom::fmt::memory_buffer out; + bool is_first = true; + + for(auto& curve : cp.getEdges()) + { + // Only write out first point for first edge + if(is_first) + { + axom::fmt::format_to(out, "M {} {} ", curve[0][0], curve[0][1]); + is_first = false; + } + + switch(curve.getOrder()) + { + case 1: + axom::fmt::format_to(out, "L {} {} ", curve[1][0], curve[1][1]); + break; + case 2: + axom::fmt::format_to(out, + "Q {} {}, {} {} ", + curve[1][0], + curve[1][1], + curve[2][0], + curve[2][1]); + break; + case 3: + axom::fmt::format_to(out, + "C {} {}, {} {}, {} {} ", + curve[1][0], + curve[1][1], + curve[2][0], + curve[2][1], + curve[3][0], + curve[3][1]); + break; + default: + SLIC_WARNING( + "Unsupported case: can only output up to cubic curves as SVG."); + } + } + return axom::fmt::format(" \n", + axom::fmt::to_string(out)); + }; + + std::string srcGroup; + std::string tgtGroup; + std::string intersectionGroup; + + // output src mesh + { + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + auto& meshWrapper = srcMesh; + for(int i = 0; i < meshWrapper.numElements(); ++i) + { + auto cp = meshWrapper.elemAsCurvedPolygon(i); + axom::fmt::format_to(out, cpToSVG(cp)); + } + axom::fmt::format_to(out, " \n"); + srcGroup = axom::fmt::to_string(out); + } + + //output tgt mesh + { + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + auto& meshWrapper = tgtMesh; + for(int i = 0; i < meshWrapper.numElements(); ++i) + { + auto cp = meshWrapper.elemAsCurvedPolygon(i); + axom::fmt::format_to(out, cpToSVG(cp)); + } + axom::fmt::format_to(out, " \n"); + tgtGroup = axom::fmt::to_string(out); + } + + //output intersection elements + { + double EPS = 1e-8; + axom::fmt::memory_buffer out; + axom::fmt::format_to( + out, + " \n"); + + // foreach target element, find src intersections and add to group + const int nTargetElems = tgtMesh.numElements(); + for(int i = 0; i < nTargetElems; ++i) + { + auto candidates = getSourceCandidates(i); + if(candidates.empty()) + { + continue; + } + + auto tgtPoly = tgtMesh.elemAsCurvedPolygon(i); + for(int srcElem : candidates) + { + auto srcPoly = srcMesh.elemAsCurvedPolygon(srcElem); + + srcPoly.reverseOrientation(); + tgtPoly.reverseOrientation(); + + std::vector pnew; + if(primal::intersect(tgtPoly, srcPoly, pnew, EPS)) + { + for(auto& cp : pnew) + { + axom::fmt::format_to(out, cpToSVG(cp)); + } + } + + srcPoly.reverseOrientation(); + tgtPoly.reverseOrientation(); + } + } + + axom::fmt::format_to(out, " \n"); + intersectionGroup = axom::fmt::to_string(out); + } + + // Write the file + { + std::string fname = "high_order_intersections.svg"; + std::ofstream fs(fname); + fs << header << std::endl; + fs << srcGroup << std::endl; + fs << tgtGroup << std::endl; + fs << intersectionGroup << std::endl; + fs << footer << std::endl; + } + } + +public: + MeshWrapper srcMesh; + MeshWrapper tgtMesh; + GridType grid; + +private: + // Named fields map manage memory associated with the + // finite element collection and spaces + mfem::NamedFieldsMap fecMap; + mfem::NamedFieldsMap fesMap; + +private: + // scale and translate the vertices of the given mesh + void xformMesh(mfem::Mesh* mesh, double sc, double off1, double off2) + { + // std::cout << "Transforming Mesh now" << std::endl; + //for (int v = 0 ; v < NumVerts ; ++v) + // { + // double pt* = mesh->GetVertex(v); + // pt[0] = sc*pt[0] + off1; + // pt[1] = sc*pt[1] + off2; + // } + mfem::GridFunction* mesh_nodes = mesh->GetNodes(); + //std::cout << *mesh_nodes[0] << std::endl; + int NumDofs = mesh_nodes->Size(); + for(int e = 0; e < (NumDofs / 2); ++e) + { + (*mesh_nodes)[2 * e] = sc * (*mesh_nodes)[2 * e] + off1; + (*mesh_nodes)[2 * e + 1] = sc * (*mesh_nodes)[2 * e + 1] + off2; + } + } +}; + +/// Simple struct to hold properties for initializing a mesh instance +struct MeshProps +{ + std::string file; + std::vector offset; + double scale {1.}; + int order {2}; + int refinement {0}; + + bool isDefault() const { return file.empty() && offset.empty(); } + bool hasOffset() const { return offset.size() == 2; } + + friend std::ostream& operator<<(std::ostream& os, const MeshProps& props) + { + os << axom::fmt::format( + "{{\n" + " file: {} \n" + " offset({}): {} \n" + " order: {} \n" + " scale: {} \n" + " refinement: {} \n" + "}}", + props.file.empty() ? "<>" : axom::fmt::format("'{}'", props.file), + props.offset.size(), + axom::fmt::join(props.offset, " "), + props.order, + props.scale, + props.refinement); + return os; + } +}; + +//------------------------------------------------------------------------------ +int main(int argc, char** argv) +{ + axom::slic::SimpleLogger logger; + + SLIC_INFO("The application conservatively maps fields from a source \n" + << "high order mesh to a target high order mesh!"); + + MeshProps srcMesh; + MeshProps tgtMesh; + + // Setup default meshes if user doesn't pass in data +#ifdef AXOM_DATA_DIR + // Use a mesh from a file if we have the data directory + namespace fs = axom::utilities::filesystem; + MeshProps defaultSrcMesh; + defaultSrcMesh.file = fs::joinPath(AXOM_DATA_DIR, "mfem/disc-nurbs.mesh"); + defaultSrcMesh.scale = 1.5; + defaultSrcMesh.offset = {.01517288412347, .02571238506182}; + defaultSrcMesh.order = 3; + + MeshProps defaultTgtMesh; + defaultTgtMesh.file = fs::joinPath(AXOM_DATA_DIR, "mfem/disc-nurbs.mesh"); + defaultTgtMesh.offset = {.001237586, -.06172376}; + defaultTgtMesh.order = 3; +#else + // parameters for a Cartesian mesh + MeshProps defaultSrcMesh; + defaultSrcMesh.refinement = 1; + defaultSrcMesh.order = 5; + + // quad mesh partially covering unit square + MeshProps defaultTgtMesh; + defaultSrcMesh.refinement = 1; + defaultTgtMesh.scale = 712378102150; + defaultTgtMesh.offset = {.1345747181586, .1345747181586}; + defaultTgtMesh.order = 5; +#endif + + // Set up and parse command line args + axom::CLI::App app {"High order mesh intersection application"}; + { + app.add_option("--srcFile", srcMesh.file) + ->description("mfem mesh file for source mesh") + ->check(axom::CLI::ExistingFile); + app.add_option("--srcOffset", srcMesh.offset) + ->description("offset for source mesh") + ->expected(2); + app.add_option("--srcScale", srcMesh.scale) + ->description("scale for source mesh") + ->capture_default_str(); + app.add_option("--srcOrder", srcMesh.order) + ->description("polynomial order for source mesh") + ->capture_default_str(); + app.add_option("--srcRef", srcMesh.refinement) + ->description("refinement levels for source mesh") + ->capture_default_str(); + + app.add_option("--tgtFile", tgtMesh.file) + ->description("mfem mesh file for source mesh") + ->check(axom::CLI::ExistingFile); + app.add_option("--tgtOffset", tgtMesh.offset) + ->description("offset for target mesh") + ->expected(2); + app.add_option("--tgtScale", tgtMesh.scale) + ->description("scale for target mesh") + ->capture_default_str(); + app.add_option("--tgtOrder", tgtMesh.order) + ->description("polynomial order for target mesh") + ->capture_default_str(); + app.add_option("--tgtRef", tgtMesh.refinement) + ->description("refinement levels for target mesh") + ->capture_default_str(); + } + CLI11_PARSE(app, argc, argv); + + SLIC_INFO("Running intersection tests..."); + axom::utilities::Timer timer(false); + + Remapper remap; + + // Load the source mesh + { + const bool isSource = true; + if(srcMesh.isDefault()) + { + srcMesh = defaultSrcMesh; + } + if(!srcMesh.hasOffset()) + { + srcMesh.offset = {0, 0}; + } + SLIC_INFO("Loading source mesh: " << srcMesh); + remap.loadMesh(isSource, + srcMesh.file, + srcMesh.offset[0], + srcMesh.offset[1], + srcMesh.scale, + srcMesh.refinement, + srcMesh.order); + } + // Load the target mesh + { + const bool isSource = false; + if(tgtMesh.isDefault()) + { + tgtMesh = defaultTgtMesh; + } + if(!tgtMesh.hasOffset()) + { + tgtMesh.offset = {0, 0}; + } + SLIC_INFO("Loading target mesh: " << tgtMesh); + remap.loadMesh(isSource, + tgtMesh.file, + tgtMesh.offset[0], + tgtMesh.offset[1], + tgtMesh.scale, + tgtMesh.refinement, + tgtMesh.order); + } + + // Set up the spatial index + remap.initSpatialIndex(); + + // Computes the overlaps between elements of the target and source meshes + timer.start(); + double area = remap.computeOverlapAreas(); + timer.stop(); + + remap.outputAsSVG(); + + SLIC_INFO(axom::fmt::format("Intersecting meshes: area: {}, time: {}", + area, + timer.elapsed())); + + return 0; +}