diff --git a/src/axom/core/numerics/Matrix.hpp b/src/axom/core/numerics/Matrix.hpp index b93ee04d87..6820657a57 100644 --- a/src/axom/core/numerics/Matrix.hpp +++ b/src/axom/core/numerics/Matrix.hpp @@ -118,6 +118,11 @@ class Matrix { public: + /*! + * \brief Default constructor + */ + Matrix() : m_rows(0), m_cols(0), m_data(nullptr), m_usingExternal(false) {} + /*! * \brief Constructor, creates a Matrix with the given rows and columns. * @@ -530,12 +535,6 @@ class Matrix private: - /*! - * \brief Default constructor. Does nothing. - * \note Made private to prevent host-code from calling this. - */ - Matrix() : m_rows(0), m_cols(0), m_data(nullptr) { }; - /// \name Private Helper Methods /// @{ diff --git a/src/axom/core/tests/utils_utilities.cpp b/src/axom/core/tests/utils_utilities.cpp index 85c092b7e8..b321462ec7 100644 --- a/src/axom/core/tests/utils_utilities.cpp +++ b/src/axom/core/tests/utils_utilities.cpp @@ -116,3 +116,90 @@ TEST(core_Utilities,minmax) EXPECT_EQ(a, axom::utilities::max(a,b)); } } + +TEST(core_Utilities, binomial_coefficient) +{ + std::cout<<"Testing binomial coefficient function."<< std::endl; + + // test n less than zero + { + const int n = -1; + const int exp = 0; + for(int k=-1 ; k < 10 ; ++k) + { + auto binom_k_n = axom::utilities::binomialCoefficient(n,k); + EXPECT_EQ( exp, binom_k_n); + } + } + + // test n := 0 + { + const int n = 0; + + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,0)); + + EXPECT_EQ( 0, axom::utilities::binomialCoefficient(n,-1)); + EXPECT_EQ( 0, axom::utilities::binomialCoefficient(n, 1)); + + } + + // test n := 1 + { + const int n = 1; + + EXPECT_EQ( 0, axom::utilities::binomialCoefficient(n,-1)); + + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n, 0)); + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n, 1)); + + EXPECT_EQ( 0, axom::utilities::binomialCoefficient(n, 2)); + + } + + // test n := 2 + { + const int n = 2; + + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,0)); + EXPECT_EQ( 2, axom::utilities::binomialCoefficient(n,1)); + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,2)); + + } + + // test n := 3 + { + const int n = 3; + + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,0)); + EXPECT_EQ( 3, axom::utilities::binomialCoefficient(n,1)); + EXPECT_EQ( 3, axom::utilities::binomialCoefficient(n,2)); + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,3)); + } + + // test n := 4 + { + const int n = 4; + + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,0)); + EXPECT_EQ( 4, axom::utilities::binomialCoefficient(n,1)); + EXPECT_EQ( 6, axom::utilities::binomialCoefficient(n,2)); + EXPECT_EQ( 4, axom::utilities::binomialCoefficient(n,3)); + EXPECT_EQ( 1, axom::utilities::binomialCoefficient(n,4)); + } + + // test recurrence relation nCk = (n-1)C(k-1) + (n-1)C(k) + { + for(int n = 1 ; n < 10 ; ++n) + { + for(int k=1 ; k <= n ; ++k) + { + auto binom_n_k = axom::utilities::binomialCoefficient(n,k); + auto binom_n1_k1 = axom::utilities::binomialCoefficient(n-1,k-1); + auto binom_n1_k = axom::utilities::binomialCoefficient(n-1,k); + + EXPECT_EQ(binom_n_k, binom_n1_k1 + binom_n1_k ); + } + } + } + +} diff --git a/src/axom/core/utilities/Utilities.cpp b/src/axom/core/utilities/Utilities.cpp index 66a05ceaa3..f7573fd25b 100644 --- a/src/axom/core/utilities/Utilities.cpp +++ b/src/axom/core/utilities/Utilities.cpp @@ -40,5 +40,29 @@ void processAbort() #endif } +int binomialCoefficient(int n, int k) +{ + if(k > n || k < 0) // check if out-of-bounds + { + return 0; + } + if(k == n || k == 0) // early return + { + return 1; + } + if (k > n-k) // exploit symmetry to reduce work + { + k= n-k; + } + + int val = 1; + for (int i=1 ; i<=k ; ++i) + { + val*=(n-k+i); + val/=i; + } + return val; +} + } // end namespace utilities } // end namespace axom diff --git a/src/axom/core/utilities/Utilities.hpp b/src/axom/core/utilities/Utilities.hpp index d891db13e2..3ef614607c 100644 --- a/src/axom/core/utilities/Utilities.hpp +++ b/src/axom/core/utilities/Utilities.hpp @@ -144,6 +144,14 @@ T clampLower(T val, T lower) return val < lower ? lower : val; } +/*! + * \brief Computes the binomial coefficient `n choose k` + * + * \return \f$ {n\choose k} = n! / (k! * (n-k)!)\f$ + * when \f$ n \ge k \ge 0 \f$, 0 otherwise. + */ +int binomialCoefficient(int n, int k); + /*! * \brief Returns a random real number within the specified interval * diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index b76851f80f..7d32e207b1 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -20,6 +20,7 @@ set( primal_headers ## geometry geometry/BezierCurve.hpp geometry/BoundingBox.hpp + geometry/CurvedPolygon.hpp geometry/OrientedBoundingBox.hpp geometry/OrientationResult.hpp geometry/NumericArray.hpp @@ -44,6 +45,7 @@ set( primal_headers operators/detail/clip_impl.hpp operators/detail/intersect_bezier_impl.hpp + operators/detail/intersect_curved_poly_impl.hpp operators/detail/intersect_impl.hpp ) diff --git a/src/axom/primal/geometry/BezierCurve.hpp b/src/axom/primal/geometry/BezierCurve.hpp index 97c1fe1354..d2eeca7459 100644 --- a/src/axom/primal/geometry/BezierCurve.hpp +++ b/src/axom/primal/geometry/BezierCurve.hpp @@ -24,8 +24,8 @@ #include "axom/primal/operators/squared_distance.hpp" -#include "fmt/fmt.hpp" #include +#include #include namespace axom @@ -42,6 +42,35 @@ template < typename T,int NDIMS > std::ostream& operator<<(std::ostream & os, const BezierCurve< T,NDIMS > & bCurve); +/*! + * \brief Computes the weights for BezierCurve's sectorArea() function + * + * \param order The polynomial order of the curve + * \return An anti-symmetric matrix with (order+1)*{order+1) entries + * containing the integration weights for entry (i,j) + * + * The derivation is provided in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ +template < typename T> +numerics::Matrix generateBezierCurveSectorWeights(int order); + +/*! + * \brief Computes the weights for BezierCurve's sectorCentroid() function + * + * \param order The polynomial order of the curve + * \return An anti-symmetric matrix with (order+1)*{order+1) entries + * containing the integration weights for entry (i,j) + * + * The derivation is provided in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ +template < typename T> +std::vector> generateBezierCurveSectorCentroidsWeights(int order); + + /*! * \class BezierCurve * @@ -66,6 +95,15 @@ class BezierCurve using BoundingBoxType = BoundingBox< T, NDIMS >; using OrientedBoundingBoxType = OrientedBoundingBox< T, NDIMS >; +private: + // Caches of precomputed coefficients for sector area calculations + // on a given polynomial order + using SectorWeights = numerics::Matrix; + using WeightsMap = std::map; + using CentroidsWeightsMap = std::map>; + static WeightsMap s_sectorWeightsMap; + static CentroidsWeightsMap s_sectorCentroidsWeightsMap; + public: /*! @@ -135,6 +173,24 @@ class BezierCurve } } + /*! + * \brief Constructor for a Bezier Curve from an vector of coordinates + * + * \param [in] pts a vector with ord+1 control points + * \param [in] ord The Curve's polynomial order + * \pre order is greater than or equal to zero + * + */ + + BezierCurve(const std::vector& pts, int ord) + { + SLIC_ASSERT(ord >= 0); + + const int sz = utilities::max(0, ord+1); + m_controlPoints.resize(sz); + m_controlPoints = pts; + } + /*! Sets the order of the Bezier Curve*/ void setOrder( int ord) { m_controlPoints.resize(ord+1); } @@ -176,6 +232,17 @@ class BezierCurve return m_controlPoints; } + /*! Reverses the order of the Bezier curve's control points */ + void reverseOrientation() + { + const int ord = getOrder(); + CoordsVec old_controlPoints= m_controlPoints; + for (int i=0 ; i<=ord ; ++i) + { + m_controlPoints[i] = old_controlPoints[ord-i]; + } + } + /*! Returns an axis-aligned bounding box containing the Bezier curve */ BoundingBoxType boundingBox() const { @@ -227,6 +294,45 @@ class BezierCurve return ptval; } + /*! + * \brief Computes the tangent of a Bezier curve at a particular parameter value \a t + * + * \param [in] t parameter value at which to compute tangent + * \return p the tangent vector of the Bezier curve at t + * + * \note We typically find the tangent of the curve at \a t between 0 and 1 + */ + + VectorType dt(T t) const + { + VectorType val; + + const int ord = getOrder(); + std::vector dCarray(ord+1); + + // Run de Casteljau algorithm on each dimension + for ( int i=0 ; i < NDIMS ; ++i) + { + for ( int p=0 ; p <= ord ; ++p) + { + dCarray[p] = m_controlPoints[p][i]; + } + + // stop one step early and take difference of last two values + for ( int p=1 ; p <= ord-1 ; ++p) + { + const int end = ord-p; + for ( int k=0 ; k <= end ; ++k) + { + dCarray[k]=(1-t)*dCarray[k] + t*dCarray[k+1]; + } + } + val[i]=ord*(dCarray[1]-dCarray[0]); + } + + return val; + } + /*! * \brief Splits a Bezier curve into two Bezier curves at particular parameter * value between 0 and 1 @@ -266,6 +372,76 @@ class BezierCurve return; } + + /* + * \brief Calculates the sector moment of a Bezier Curve + * + * The sector moment is the moment between the curve and the origin. + * The equation and derivation are generalizations of: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ + PointType sectorCentroid() const + { + T Mx = 0; + T My = 0; + const int ord = getOrder(); + + // Compute and cache the weights if they are not already available + if (s_sectorCentroidsWeightsMap.find(ord)==s_sectorCentroidsWeightsMap.end()) + { + auto wts = generateBezierCurveSectorCentroidsWeights(ord); + s_sectorCentroidsWeightsMap.emplace(std::make_pair(ord,wts)); + } + + const auto& weights = s_sectorCentroidsWeightsMap[ord]; + for (int r=0 ; r<=ord ; ++r) + { + for (int p=0 ; p<=ord ; ++p) + { + for (int q=0 ; q<=ord ; ++q) + { + Mx+= weights[r](p,q)*m_controlPoints[p][1]* m_controlPoints[q][0]* m_controlPoints[r][0]; + My+= weights[r](p,q)*m_controlPoints[p][1]* m_controlPoints[q][0]* m_controlPoints[r][1]; + } + } + } + PointType M= PointType::make_point(Mx,My); + return M; + } + + /* + * \brief Calculates the sector area of a Bezier Curve + * + * The sector area is the area between the curve and the origin. + * The equation and derivation is described in: + * Ueda, K. "Signed area of sectors between spline curves and the origin" + * IEEE International Conference on Information Visualization, 1999. + */ + T sectorArea() const + { + T A = 0; + const int ord = getOrder(); + + // Compute and cache the weights if they are not already available + if (s_sectorWeightsMap.find(ord)==s_sectorWeightsMap.end()) + { + auto wts = generateBezierCurveSectorWeights(ord); + s_sectorWeightsMap.emplace(std::make_pair(ord,wts)); + } + + const auto& weights = s_sectorWeightsMap[ord]; + + for (int p=0 ; p<=ord ; ++p) + { + for (int q=0 ; q<=ord ; ++q) + { + A+= weights(p,q)*m_controlPoints[p][1]* m_controlPoints[q][0]; + } + } + return A; + } + /*! * \brief Predicate to check if the Bezier curve is approximately linear * @@ -318,8 +494,19 @@ class BezierCurve CoordsVec m_controlPoints; }; +// Declaration of sectorArea weights map +template < typename T, int NDIMS > +typename BezierCurve::WeightsMap +BezierCurve::s_sectorWeightsMap; + +// Declaration of sectorCentroids weights map +template < typename T, int NDIMS > +typename BezierCurve::CentroidsWeightsMap +BezierCurve::s_sectorCentroidsWeightsMap; + + //------------------------------------------------------------------------------ -/// Free functions implementing BezierCurve's operators +/// Free functions related to BezierCurve //------------------------------------------------------------------------------ template < typename T, int NDIMS > std::ostream& operator<<(std::ostream & os, @@ -329,6 +516,73 @@ std::ostream& operator<<(std::ostream & os, return os; } +template +numerics::Matrix generateBezierCurveSectorWeights(int ord) +{ + numerics::Matrix weights(ord+1,ord+1); + T binom_2n_n = static_cast(utilities::binomialCoefficient(2*ord,ord)); + for (int i=0 ; i<=ord ; ++i) + { + weights(i,i) = 0.; // zero on the diagonal + for (int j=i+1 ; j<=ord ; ++j) + { + double val = 0.; + if (i != j) + { + T binom_ij_i = + static_cast(utilities::binomialCoefficient(i+j,i)); + T binom_2nij_nj = + static_cast(utilities::binomialCoefficient(2*ord-i-j,ord-j)); + + val = ((j-i) * ord) / binom_2n_n + * (binom_ij_i / static_cast(i+j)) + * (binom_2nij_nj / (2.*ord -j-i)); + } + weights(i,j) = val; // antisymmetric + weights(j,i) = -val; + } + } + return weights; +} + +template +std::vector> generateBezierCurveSectorCentroidsWeights(int ord) +{ + std::vector> weights; + weights.resize(ord+1); + numerics::Matrix weightsr(ord+1,ord+1); + for (int k=0 ; k<=ord ; ++k) + { + for (int i=0 ; i<=ord ; ++i) + { + weightsr(i,i) = 0.; // zero on the diagonal + for (int j=i+1 ; j<=ord ; ++j) + { + double val = 0.; + if (i != j) + { + T binom_n_i = + static_cast(utilities::binomialCoefficient(ord,i)); + T binom_n_j = + static_cast(utilities::binomialCoefficient(ord,j)); + T binom_n_k = + static_cast(utilities::binomialCoefficient(ord,k)); + T binom_3n2_ijk1 = + static_cast(utilities::binomialCoefficient(3*ord-2,i+j+k-1)); + + val = (1.*(j-i) ) / (3.*(3*ord-1)) + * (1.*binom_n_i*binom_n_j*binom_n_k/ (1.*binom_3n2_ijk1)); + } + weightsr(i,j) = val; // antisymmetric + weightsr(j,i) = -val; + } + } + weights[k]=weightsr; + } + return weights; +} + + } // namespace primal } // namespace axom diff --git a/src/axom/primal/geometry/CurvedPolygon.hpp b/src/axom/primal/geometry/CurvedPolygon.hpp new file mode 100644 index 0000000000..7a4385d28b --- /dev/null +++ b/src/axom/primal/geometry/CurvedPolygon.hpp @@ -0,0 +1,271 @@ +// Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file CurvedPolygon.hpp + * + * \brief A CurvedPolygon primitive whose edges are Bezier Curves + */ + +#ifndef PRIMAL_CURVEDPOLYGON_HPP_ +#define PRIMAL_CURVEDPOLYGON_HPP_ + +#include "axom/slic.hpp" + +#include "axom/primal/geometry/Point.hpp" +#include "axom/primal/geometry/Vector.hpp" +#include "axom/primal/geometry/NumericArray.hpp" +#include "axom/primal/geometry/BezierCurve.hpp" + +#include +#include // for std::ostream + +namespace axom +{ +namespace primal +{ + +// Forward declare the templated classes and operator functions +template < typename T, int NDIMS > +class CurvedPolygon; + +/*! \brief Overloaded output operator for polygons */ +template < typename T,int NDIMS > +std::ostream& operator<<(std::ostream & os, const CurvedPolygon< T,NDIMS > & poly); + +/*! + * \class CurvedPolygon + * + * \brief Represents a polygon with curved edges defined by BezierCurves + * \tparam T the coordinate type, e.g., double, float, etc. + * \tparam NDIMS the number of dimensions + * \note The component curves should be ordered in a counter clockwise + * orientation with respect to the polygon's normal vector + */ +template < typename T,int NDIMS > +class CurvedPolygon +{ +public: + using PointType = Point< T,NDIMS >; + using VectorType = Vector< T,NDIMS >; + using NumArrayType = NumericArray< T,NDIMS >; + using BezierCurveType = BezierCurve; +public: + /*! Default constructor for an empty polygon */ + CurvedPolygon() = default; + + /*! + * \brief Constructor for an empty CurvedPolygon that reserves space for + * the given number of Edges + * + * \param [in] numExpectedEdges number of edges for which to reserve space + * \pre numExpectedEdges is at least 1 + * + */ + explicit CurvedPolygon(int nEdges) + { + SLIC_ASSERT(nEdges >= 1); + m_edges.reserve(nEdges); + m_edges.resize(nEdges); + } + + CurvedPolygon(BezierCurveType* curves, int nEdges) + { + SLIC_ASSERT(curves != nullptr); + SLIC_ASSERT(nEdges >=1); + + m_edges.reserve(nEdges); + + for (int e =0; eaddEdge(curves[e]); + } + } + + /*! Return the number of edges in the polygon */ + int numEdges() const { return m_edges.size(); } + + void setNumEdges(int ngon) + { + SLIC_ASSERT(ngon >= 0); + m_edges.resize(ngon); + } + + /* Checks equality of two Bezier Curve */ + friend inline bool operator==(const CurvedPolygon< T, NDIMS>& lhs, + const CurvedPolygon< T, NDIMS>& rhs) + { + return lhs.m_edges == rhs.m_edges; + } + + friend inline bool operator!=(const CurvedPolygon< T, NDIMS>& lhs, + const CurvedPolygon< T, NDIMS>& rhs) + { + return !(lhs == rhs); + } + + /*! Appends a BezierCurve to the list of edges */ + void addEdge(const BezierCurveType& c1) + { + m_edges.push_back(c1); + } + + /*! Splits an edge "in place" */ + void splitEdge(int idx, T t) + { + SLIC_ASSERT(idx> getEdges() const + { + return m_edges; + } + + /*! Retrieves the Bezier Curve at index idx */ + BezierCurveType& operator[](int idx) { return m_edges[idx]; } + /*! Retrieves the vertex at index idx */ + const BezierCurveType& operator[](int idx) const { return m_edges[idx]; } + + /*! + * \brief Simple formatted print of a CurvedPolygon instance + * + * \param os The output stream to write to + * \return A reference to the modified ostream + */ + std::ostream& print(std::ostream& os) const + { + const int sz = numEdges(); + + os <<"{" << sz <<"-sided Bezier polygon:"; + for (int i=0 ; i< sz-1 ; ++i) + { + os << m_edges[i] << ","; + } + if (sz >= 2) + { + os<> old_edges=m_edges; + for (int i=0 ; i > m_edges; +}; + + +//------------------------------------------------------------------------------ +/// Free functions implementing Polygon's operators +//------------------------------------------------------------------------------ +template < typename T, int NDIMS > +std::ostream& operator<<(std::ostream & os, const CurvedPolygon< T,NDIMS > & poly) +{ + poly.print(os); + return os; +} + +} // namespace primal +} // namespace axom + +#endif // PRIMAL_CURVEDPOLYGON_HPP_ diff --git a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp index 11db161e95..511e8f0bba 100644 --- a/src/axom/primal/operators/detail/intersect_bezier_impl.hpp +++ b/src/axom/primal/operators/detail/intersect_bezier_impl.hpp @@ -86,6 +86,12 @@ bool intersect_bezier_curves( const BezierCurve< T, NDIMS>& 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 + * one, i.e. there are no points at which the curves and their derivatives + * both intersect. Thus, the function does not find tangencies. + * + * \note This function assumes two dimensional curves in a plane. + * * \note This function does not properly handle collinear lines */ @@ -111,6 +117,7 @@ bool intersect_bezier_curves( const BezierCurve< T, NDIMS>& c1, double t_offset, double t_scale) { using BCurve = BezierCurve< T, NDIMS>; + SLIC_ASSERT(NDIMS==2); // Check bounding boxes to short-circuit the intersection if(!intersect(c1.boundingBox(), c2.boundingBox()) ) @@ -120,7 +127,7 @@ bool intersect_bezier_curves( const BezierCurve< T, NDIMS>& c1, bool foundIntersection = false; - if ( c1.isLinear(sq_tol) && c2.isLinear(sq_tol)) + if ( c1.isLinear(sq_tol) && c2.isLinear(sq_tol)) { T s,t; if(intersect_2d_linear(c1[0],c1[order1],c2[0],c2[order2],s,t)) 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..bcde15cedf --- /dev/null +++ b/src/axom/primal/operators/detail/intersect_curved_poly_impl.hpp @@ -0,0 +1,333 @@ +// Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file intersect.hpp + * + * \brief Consists of functions to test intersection among geometric primitives. + */ + +#ifndef INTERSECTION_CURVED_POLYGON_HPP_ +#define INTERSECTION_CURVED_POLYGON_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/core/utilities/Utilities.hpp" + +#include "axom/primal/operators/squared_distance.hpp" +#include "axom/primal/operators/detail/intersect_bezier_impl.hpp" + +namespace axom +{ +namespace primal +{ +namespace detail +{ +template < typename T > +class IntersectionInfo; + +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); + +/*! + * \brief Test whether the regions within CurvedPolygons p1 and p2 intersect. + * \return status true iff p1 intersects with 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 < typename T, int NDIMS> +bool intersect_polygon( const CurvedPolygon< T, NDIMS>& p1, + const CurvedPolygon< T, NDIMS>& p2, + std::vector>& pnew, + double sq_tol + ) +{ + // Object to store intersections + std::vector>> E1IntData(p1.numEdges()); + std::vector>> E2IntData(p2.numEdges()); + IntersectionInfo firstinter; // Need to do orientation test on first intersection + + // Find all intersections and store + int 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.); + for (int k =0; k< static_cast(p1times.size()); ++k) + { + E1IntData[i].push_back({p1times[k],i,p2times[k],j,numinters+k+1}); + E2IntData[j].push_back({p2times[k],j,p1times[k],i,numinters+k+1}); + if (numinters==0) + { + firstinter={p1times[0],i,p2times[0],j,1}; + } + } + numinters+=p1times.size(); + } + } + if (numinters>0) + { + for (int i = 0; i < p1.numEdges(); ++i) + { + std::sort( E1IntData[i].begin(), E1IntData[i].end() ); + std::sort( E2IntData[i].begin(), E2IntData[i].end() ); + } + + // Orient the first intersection point to be sure we get the intersection + bool orientation = detail::orient(p1[firstinter.myEdge],p2[firstinter.otherEdge],firstinter.myTime,firstinter.otherTime); + + // 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. + std::vector edgelabels[2]; // 0 for curves that end in original vertices, unique id for curves that end in intersection points + CurvedPolygon psplit[2]; // The two completely split polygons will be stored in this array + psplit[0]=p1; + psplit[1]=p2; + + splitPolygon(psplit[0],E1IntData,edgelabels[0]); + splitPolygon(psplit[1],E2IntData,edgelabels[1]); + + // This performs the directional walking method using the completely split polygons + std::vector::iterator> usedlabels; + bool addingcurves=true; // When this is false, we are walking between intersection regions + int startvertex=1; // Start at the vertex with "unique id" 1 + int nextvertex; // The next vertex id is unknown at this time + bool currentelement=orientation; // This variable allows us to switch between the two elements + int currentit= std::find(edgelabels[currentelement].begin(),edgelabels[currentelement].end(),startvertex)-edgelabels[currentelement].begin(); // This is the iterator pointing to the end vertex of the edge of the completely split polygon we are on + int startit = currentit; // This is the iterator to the end vertex of the starting edge on the starting polygon + int nextit = (currentit+1)%edgelabels[0].size(); // This is the iterator to the end vertex of the next edge of whichever polygon we will be on next + nextvertex=edgelabels[currentelement][nextit]; // This is the next vertex id + while (numinters>0) + { + CurvedPolygon aPart; // Object to store the current intersection polygon (could be multiple) + while (!(nextit==startit && currentelement==orientation) || addingcurves==false ) // Once the end vertex of the current edge is the start vertex, we need to switch regions + { + if (nextit==currentit ) + { + nextit=(currentit+1)%edgelabels[0].size(); + } + nextvertex=edgelabels[currentelement][nextit]; + while (nextvertex==0) + { + currentit=nextit; + if (addingcurves) + { + aPart.addEdge(psplit[currentelement][nextit]); + } + nextit=(currentit+1)%edgelabels[0].size(); + nextvertex=edgelabels[currentelement][nextit]; + } + if (edgelabels[currentelement][nextit]>0) + { + if (addingcurves) + { + aPart.addEdge(psplit[currentelement][nextit]); + edgelabels[currentelement][nextit]=-edgelabels[currentelement][nextit]; + currentelement=!currentelement; + nextit = std::find(edgelabels[currentelement].begin(),edgelabels[currentelement].end(),nextvertex)-edgelabels[currentelement].begin(); + edgelabels[currentelement][nextit]=-edgelabels[currentelement][nextit]; + currentit=nextit; + numinters-=1; + } + else + { + addingcurves=true; + startit=nextit; + currentit=nextit; + nextit=(currentit+1)%edgelabels[0].size(); + orientation=currentelement; + } + } + else + { + currentelement=!currentelement; + nextit = std::find(edgelabels[currentelement].begin(),edgelabels[currentelement].end(),nextvertex)-edgelabels[currentelement].begin(); + edgelabels[currentelement][nextit]=-edgelabels[currentelement][nextit]; + currentit=nextit; + } + } + pnew.push_back(aPart); + currentelement=!currentelement; + currentit = std::find(edgelabels[currentelement].begin(),edgelabels[currentelement].end(),-nextvertex)-edgelabels[currentelement].begin(); + nextit=(currentit+1)%edgelabels[0].size(); + if (numinters>0) + { + addingcurves=false; + } + } + 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< T, NDIMS >; + using PointType = typename BCurve::PointType; + using VectorType = typename BCurve::VectorType; + + 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 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(temps.size()) ; ++i) + { + if (temps[i]>line1s) + { + line1s = temps[i]; + p1c = j; + p1t = tempt[i]; + } + } + } + for (int j=0 ; j 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(temps.size()) ; ++i) + { + if (temps[i]line1s) + { + line2s = temps[i]; + p2c = j; + p2t = tempt[i]; + } + } + } + + bool E1inE2 = VectorType::cross_product(p1[p1c].dt(p1t),LineGuess.dt(line1s))[2]<0; + bool E2inE1 = VectorType::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 Splits a CurvedPolygon p1 at every intersection point stored in IntersectionData + */ +template +void splitPolygon( CurvedPolygon< T, NDIMS>& p1, + std::vector>>& IntersectionData, + std::vector& edgelabels) +{ + const int nEd = p1.numEdges(); + //split polygon 2 at all the intersection points and store as psplit[1] + int addedints=0; + for (int i = 0; i < nEd; ++i) + { + edgelabels.push_back(0); + for (int j = 0; j < static_cast(IntersectionData[i].size()); ++j) + { + p1.splitEdge(i+addedints,IntersectionData[i][j].myTime); + edgelabels.insert(edgelabels.begin()+i+addedints,IntersectionData[i][j].numinter); + addedints+=1; + for (int k = j+1; k < static_cast(IntersectionData[i].size()); ++k) + { + IntersectionData[i][k].myTime=(IntersectionData[i][k].myTime-IntersectionData[i][j].myTime)/(1-IntersectionData[i][j].myTime); + } + } + } +} + +/* \brief Determines orientation of a bezier curve c1 with respect to another bezier curve c2, given that they intersect at parameter values s and t + * + * \param [in] c1 the first bezier curve + * \param [in] c2 the second bezier curve + * \param [in] s the parameter value of intersection on c1 + * \param [in] t the parameter value of intersection on c2 + * \return True if the c1's positive direction is counterclockwise from c2's positive direction + */ +template +bool orient(const BezierCurve& c1, const BezierCurve& c2, T s, T t) +{ + using VectorType = primal::Vector; + + auto orientation = VectorType::cross_product(c1.dt(s), c2.dt(t))[2]; + return (orientation>0); +} + +/*! \class IntersectionInfo + * + * \brief For storing intersection points between CurvedPolygons so they can be easily sorted by parameter value using std::sort + */ +template +class IntersectionInfo +{ + public: + 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 IntersectionInfo& other) const + { + return myTime& b1, } +/*! + * \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 < typename T, int NDIMS> +bool intersect( const CurvedPolygon< T, NDIMS>& p1, + const CurvedPolygon< T, NDIMS>& p2, + std::vector< CurvedPolygon< T, NDIMS>>& 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. @@ -322,10 +363,16 @@ bool intersect(const OrientedBoundingBox< T, 3 >& b1, * * Finds all intersection points between the two curves. * + * \note This function assumes two dimensional curves in a plane. + * * \note This function assumes that the curves are in general position. * Specifically, we assume that all intersections are at points and that * the curves don't overlap. * + * \note This function assumes the 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. + * * \note This function assumes that the curves are half-open, i.e. they * contain their first endpoint, but not their last endpoint. Thus, the * curves do not intersect at \f$ s==1 \f$ or at \f$ t==1 \f$. @@ -337,6 +384,7 @@ bool intersect( const BezierCurve< T, NDIMS>& c1, std::vector< T >& tp, double tol = 1E-8) { + const double offset = 0.; const double scale = 1.; @@ -347,7 +395,6 @@ bool intersect( const BezierCurve< T, NDIMS>& c1, c1.getOrder(), c2.getOrder(), offset, scale, offset, scale); } - } /* namespace primal */ } /* namespace axom */ diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index fd1c3b8ca7..be2b00af83 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -26,6 +26,8 @@ set( primal_tests primal_vector.cpp primal_bezier_curve.cpp primal_bezier_intersect.cpp + primal_curved_polygon.cpp + primal_curved_polygon_intersect.cpp ) set(primal_test_depends fmt primal slic mint gtest) diff --git a/src/axom/primal/tests/primal_bezier_curve.cpp b/src/axom/primal/tests/primal_bezier_curve.cpp index a99bc79869..923f2f45e1 100644 --- a/src/axom/primal/tests/primal_bezier_curve.cpp +++ b/src/axom/primal/tests/primal_bezier_curve.cpp @@ -56,7 +56,7 @@ TEST( primal_beziercurve, set_order ) EXPECT_EQ(-1, bCurve.getOrder()); const int order = 1; - PointType controlPoints[] = { + PointType controlPoints[2] = { PointType::make_point(0.6, 1.2, 1.0), PointType::make_point(0.0, 1.6, 1.8) }; @@ -167,6 +167,127 @@ TEST( primal_beziercurve, evaluate) } } +//------------------------------------------------------------------------------ +TEST( primal_beziercurve_, tangent) +{ + SLIC_INFO("Testing Bezier tangent calculation"); + + const int DIM = 3; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using VectorType = primal::Vector< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + const int order = 3; + PointType data[order+1] = { PointType::make_point(0.6, 1.2, 1.0), + PointType::make_point(1.3, 1.6, 1.8), + PointType::make_point(2.9, 2.4, 2.3), + PointType::make_point(3.2, 3.5, 3.0) }; + + BezierCurveType b2Curve(data, order); + + VectorType midtval = VectorType::make_vector(3.15,2.325,1.875); + VectorType starttval = VectorType::make_vector(2.1,1.2,2.4); + VectorType endtval = VectorType::make_vector(.9,3.3,2.1); + + // Evaluate the curve at several parameter values + // Curve should be tangent to control net at endpoints + VectorType eval0 = b2Curve.dt(0.0); + VectorType eval1 = b2Curve.dt(1.0); + VectorType evalMid = b2Curve.dt(0.5); + + for ( int i=0 ; i; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + { + SLIC_INFO("Testing Bezier sector area calculation for a cubic"); + const int order = 3; + PointType data[order+1] = { PointType::make_point(0.6, 1.2), + PointType::make_point(1.3, 1.6), + PointType::make_point(2.9, 2.4), + PointType::make_point(3.2, 3.5) }; + + BezierCurveType bCurve(data, order); + EXPECT_TRUE(axom::utilities::isNearlyEqual(bCurve.sectorArea(),.1455)); + } +} + +//------------------------------------------------------------------------------ +TEST( primal_beziercurve, sector_moment_cubic) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + { + SLIC_INFO("Testing Bezier sector moment calculation for a cubic"); + const int order = 3; + PointType data[order+1] = { PointType::make_point(0.6, 1.2), + PointType::make_point(1.3, 1.6), + PointType::make_point(2.9, 2.4), + PointType::make_point(3.2, 3.5) }; + + BezierCurveType bCurve(data, order); + PointType M = bCurve.sectorCentroid(); + EXPECT_NEAR(M[0],-.429321428571429,2e-15); + EXPECT_NEAR(M[1],-.354010714285715,2e-15); + } +} + + +//------------------------------------------------------------------------------ +TEST( primal_beziercurve, sector_area_point ) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + { + SLIC_INFO("Testing Bezier sector area calculation for a point"); + const int order = 0; + PointType data[order+1] = { PointType::make_point(0.6, 1.2)}; + + BezierCurveType bCurve(data, order); + EXPECT_DOUBLE_EQ(bCurve.sectorArea(),0.0); + } +} + + +//------------------------------------------------------------------------------ +TEST( primal_beziercurve, sector_moment_point ) +{ + const int DIM = 2; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + { + SLIC_INFO("Testing Bezier sector moment calculation for a point"); + const int order = 0; + PointType data[order+1] = { PointType::make_point(0.6, 1.2)}; + + BezierCurveType bCurve(data, order); + PointType M = bCurve.sectorCentroid(); + EXPECT_DOUBLE_EQ(M[0],0.0); + EXPECT_DOUBLE_EQ(M[1],0.0); + } +} + //------------------------------------------------------------------------------ TEST( primal_beziercurve, split_cubic ) { @@ -408,6 +529,123 @@ TEST( primal_beziercurve, isLinear) } +TEST( primal_beziercurve, sector_weights) +{ + SLIC_INFO("Testing weights for BezierCurve::sectorArea()"); + + // NOTE: Expected weights are provided in the reference paper [Ueda99] + // See doxygen comment for BezierCurve::sectorArea() + + using CoordType = double; + + // order 1 + { + const int ord = 1; + auto weights = primal::generateBezierCurveSectorWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(2,1); + axom::numerics::Matrix exp(ord+1,ord+1); + exp(0,0) = 0; exp(0,1) = 1; + exp(1,0) = -1; exp(1,1) = 0; + + for (int i=0 ; i<=ord ; ++i) + { + for (int j=0 ; j<=ord ; ++j) + { + EXPECT_DOUBLE_EQ( exp(i,j) * binomInv, weights(i,j) ); + } + } + } + + // order 2 + { + const int ord = 2; + auto weights = primal::generateBezierCurveSectorWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(4,2); + axom::numerics::Matrix exp(ord+1,ord+1); + exp(0,0) = 0; exp(0,1) = 2; exp(0,2) = 1; + exp(1,0) = -2; exp(1,1) = 0; exp(1,2) = 2; + exp(2,0) = -1; exp(2,1) = -2; exp(2,2) = 0; + + for (int i=0 ; i<=ord ; ++i) + { + for (int j=0 ; j<=ord ; ++j) + { + EXPECT_DOUBLE_EQ( exp(i,j) * binomInv, weights(i,j) ); + } + } + } + + // order 3 + { + const int ord = 3; + auto weights = primal::generateBezierCurveSectorWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(6,3); + axom::numerics::Matrix exp(ord+1,ord+1); + exp(0,0) = 0; exp(0,1) = 6; exp(0,2) = 3; exp(0,3) = 1; + exp(1,0) = -6; exp(1,1) = 0; exp(1,2) = 3; exp(1,3) = 3; + exp(2,0) = -3; exp(2,1) = -3; exp(2,2) = 0; exp(2,3) = 6; + exp(3,0) = -1; exp(3,1) = -3; exp(3,2) = -6; exp(3,3) = 0; + + for (int i=0 ; i<=ord ; ++i) + { + for (int j=0 ; j<=ord ; ++j) + { + EXPECT_DOUBLE_EQ( exp(i,j) * binomInv, weights(i,j) ); + } + } + } + + // order 4 + { + const int ord = 4; + auto weights = primal::generateBezierCurveSectorWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(8,4); + axom::numerics::Matrix exp(ord+1,ord+1); + exp(0,0) = 0; exp(0,1) = 20; exp(0,2) = 10; exp(0,3) = 4; exp(0,4) = 1; + exp(1,0) =-20; exp(1,1) = 0; exp(1,2) = 8; exp(1,3) = 8; exp(1,4) = 4; + exp(2,0) =-10; exp(2,1) = -8; exp(2,2) = 0; exp(2,3) = 8; exp(2,4) = 10; + exp(3,0) = -4; exp(3,1) = -8; exp(3,2) = -8; exp(3,3) = 0; exp(3,4) = 20; + exp(4,0) = -1; exp(4,1) = -4; exp(4,2) =-10; exp(4,3) =-20; exp(4,4) = 0; + + for (int i=0 ; i<=ord ; ++i) + { + for (int j=0 ; j<=ord ; ++j) + { + EXPECT_DOUBLE_EQ( exp(i,j) * binomInv, weights(i,j) ); + } + } + } + + + // order 5 + { + const int ord = 5; + auto weights = primal::generateBezierCurveSectorWeights(ord); + + double binomInv = 1. / axom::utilities::binomialCoefficient(10,5); + axom::numerics::Matrix exp(ord+1,ord+1); + exp(0,0) = 0; exp(0,1) = 70; exp(0,2) = 35; exp(0,3) = 15; exp(0,4) = 5; exp(0,5) = 1; + exp(1,0) =-70; exp(1,1) = 0; exp(1,2) = 25; exp(1,3) = 25; exp(1,4) = 15; exp(1,5) = 5; + exp(2,0) =-35; exp(2,1) =-25; exp(2,2) = 0; exp(2,3) = 20; exp(2,4) = 25; exp(2,5) = 15; + exp(3,0) =-15; exp(3,1) =-25; exp(3,2) =-20; exp(3,3) = 0; exp(3,4) = 25; exp(3,5) = 35; + exp(4,0) = -5; exp(4,1) =-15; exp(4,2) =-25; exp(4,3) =-25; exp(4,4) = 0; exp(4,5) = 70; + exp(5,0) = -1; exp(5,1) = -5; exp(5,2) =-15; exp(5,3) =-35; exp(5,4) =-70; exp(5,5) = 0; + + for (int i=0 ; i<=ord ; ++i) + { + for (int j=0 ; j<=ord ; ++j) + { + EXPECT_DOUBLE_EQ( exp(i,j) * binomInv, weights(i,j) ); + } + } + } +} + + //------------------------------------------------------------------------------ int main(int argc, char* argv[]) diff --git a/src/axom/primal/tests/primal_bezier_intersect.cpp b/src/axom/primal/tests/primal_bezier_intersect.cpp index 1df671987f..8061de2d74 100644 --- a/src/axom/primal/tests/primal_bezier_intersect.cpp +++ b/src/axom/primal/tests/primal_bezier_intersect.cpp @@ -199,7 +199,7 @@ TEST( primal_bezier_inter, linear_bezier) TEST( primal_bezier_inter, linear_bezier_interp_params) { - static const int DIM =2; +static const int DIM =2; using CoordType = double; using PointType = primal::Point< CoordType, DIM >; @@ -247,6 +247,99 @@ TEST( primal_bezier_inter, linear_bezier_interp_params) } } +//------------------------------------------------------------------------------ +TEST( primal_bezier_inter, no_intersections_bezier ) +{ + static const int DIM =2; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("primal: testing bezier intersection"); + SCOPED_TRACE("no intersections"); + + const int order = 3; + + // cubic line + PointType data1[order+1] = { PointType::make_point(0.0, 0.0), + PointType::make_point(1.0, 0.0), + PointType::make_point(2.0, 0.0), + PointType::make_point(3.0, 0.0)}; + + BezierCurveType curve1(data1, order); + + // Cubic curve + PointType data2[order+1] = { PointType::make_point(0.0, 0.5), + PointType::make_point(1.0,1.0), + PointType::make_point(2.0, 3.0), + PointType::make_point(3.0,1.5)}; + BezierCurveType curve2(data2, order); + + std::vector exp_intersections; + + const double eps = 1E-16; + const double eps_test = 1E-10; + + checkIntersections(curve1, curve2, + exp_intersections, exp_intersections, + eps, eps_test); +} + +//------------------------------------------------------------------------------ +TEST( primal_bezier_inter, cubic_quadratic_bezier ) +{ + static const int DIM =2; + using CoordType = double; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("primal: testing bezier intersection"); + SLIC_INFO("different orders"); + + const int order2 = 3; + + // cubic line + PointType data1 = PointType::make_point(0.0, 0.0); + + BezierCurveType curve1(0); + curve1[0]=data1; + + // Cubic curve + PointType data2[order2+1] = { PointType::make_point(0.0, 0.5), + PointType::make_point(1.0,-1.0), + PointType::make_point(2.0, 1.0), + PointType::make_point(3.0,-0.5)}; + BezierCurveType curve2(data2, order2); + + // Note: same intersection params for curve and line + std::vector exp_intersections = { 0.17267316464601146, + 0.5, + 0.827326835353989}; + + const double eps = 1E-16; + const double eps_test = 1E-10; + + for(int otherorder=1 ; otherorder<=20 ; ++otherorder) + { + curve1.setOrder(otherorder); + for (int i=0; i + +namespace primal = axom::primal; + +/** + * Helper function to compute the area and centroid of a curved polygon and to check that they match expectations, stored in \a expArea and \a expCentroid. Areas and Moments are computed within tolerance \a eps and checks use \a test_eps. + */ +template +void checkMoments(const primal::CurvedPolygon& bPolygon, + const CoordType expArea, const primal::Point& expMoment, + double eps, double test_eps) +{ + EXPECT_NEAR(expArea, bPolygon.area(eps),test_eps); + for (int i= 0; i< DIM; ++i) + { + EXPECT_NEAR(expMoment[i], bPolygon.centroid(eps)[i],test_eps); + } +} + +/* 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< CoordType, DIM>; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + 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 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, constructor ) +{ + const int DIM = 3; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + { + SLIC_INFO("Testing default CurvedPolygon constructor "); + CurvedPolygonType bPolygon; + + int expNumEdges = 0; + EXPECT_EQ(expNumEdges, bPolygon.numEdges()); + EXPECT_EQ(expNumEdges, bPolygon.getEdges().size()); + EXPECT_EQ(std::vector(), bPolygon.getEdges()); + } + + { + SLIC_INFO("Testing CurvedPolygon numEdges constructor "); + + CurvedPolygonType bPolygon(1); + int expNumEdges = 1; + EXPECT_EQ(expNumEdges, bPolygon.numEdges()); + EXPECT_EQ(expNumEdges, static_cast(bPolygon.getEdges().size())); + } +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, add_edges ) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test adding edges to empty CurvedPolygon"); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[2] = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.0, 1.6) + }; + + BezierCurveType bCurve(controlPoints,1); + + bPolygon.addEdge(bCurve); + bPolygon.addEdge(bCurve); + + EXPECT_EQ(2, bPolygon.numEdges()); + for (int p=0 ; p; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test checking if CurvedPolygon is closed."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_EQ(false,bPolygon.isClosed()); + + std::vector CP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.6, 1.2) + }; + std::vector orders = {1,1,1}; + + std::vector subCP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.3, 2.0) + }; + std::vector suborders = {1}; + CurvedPolygonType subPolygon = createPolygon(subCP,suborders); + EXPECT_EQ(false,subPolygon.isClosed()); + + bPolygon = createPolygon(CP,orders); + + EXPECT_EQ(3,bPolygon.numEdges()); + EXPECT_EQ(true,bPolygon.isClosed()); + bPolygon[2][1][0]-=2e-15; + EXPECT_EQ(false,bPolygon.isClosed(1e-15)); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, split_edge) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon edge split."); + + std::vector CP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.6, 1.2) + }; + + std::vector orders32 = {1,1,1}; + CurvedPolygonType bPolygon32 = createPolygon(CP,orders32); + std::cout << "Got here!! " << std::endl; + std::vector subCP; + + subCP.assign(CP.begin(), CP.begin()+2); + BezierCurveType bCurve(subCP,1); + bPolygon32.splitEdge(0,.5); + + BezierCurveType bCurve2; BezierCurveType bCurve3; + bCurve.split(.5,bCurve2,bCurve3); + + EXPECT_EQ(bPolygon32.numEdges(),4); + for (int i=0; i< bPolygon32[0].getOrder(); ++i) + { + for (int dimi=0; dimi< DIM; ++dimi) + { + EXPECT_EQ(bPolygon32[0][i][dimi],bCurve2[i][dimi]); + EXPECT_EQ(bPolygon32[1][i][dimi],bCurve3[i][dimi]); + } + } +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, moments_triangle_degenerate) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon degenerate triangle area computation."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + EXPECT_EQ(0.0, bPolygon.area()); + PointType origin = PointType::make_point(0.0,0.0); + + PointType controlPoints[2] = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.3, 2.0) + }; + + PointType controlPoints2[2] = { + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0, 1.6) + }; + PointType controlPoints3[2] = { + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.6, 1.2) + }; + + BezierCurveType bCurve(controlPoints,1); + bPolygon.addEdge(bCurve); + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); + + BezierCurveType bCurve2(controlPoints2,1); + bPolygon.addEdge(bCurve2); + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); + + BezierCurveType bCurve3(controlPoints3,1); + bPolygon.addEdge(bCurve3); + + bPolygon[2][1][0]-=1e-11; + checkMoments(bPolygon, 0.0, origin, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, moments_triangle_linear) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon linear triangle moment computation."); + std::vector CP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.6, 1.2) + }; + + std::vector orders = {1,1,1}; + CurvedPolygonType bPolygon = createPolygon(CP,orders); + + CoordType trueA= -.18; + PointType trueC= PointType::make_point(0.3,1.6); + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, moments_triangle_quadratic ) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon quadratic triangle area computation."); + + + SLIC_INFO("Test checking CurvedPolygon mixed order triangle moment computation."); + + std::vector CP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.4, 1.3), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.27, 1.5), + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.1 , 1.5), + PointType::make_point( 0.6, 1.2) + }; + + std::vector orders = {2,2,2}; + CurvedPolygonType bPolygon = createPolygon(CP,orders); + + CoordType trueA = -0.097333333333333; + PointType trueC= PointType::make_point(.294479452054794,1.548219178082190); + + checkMoments(bPolygon, trueA, trueC, 1e-15, 1e-14); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, moments_triangle_mixed_order) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon mixed order triangle moment computation."); + + std::vector CP = { + PointType::make_point( 0.6, 1.2), + PointType::make_point( 0.4, 1.3), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.27, 1.5), + PointType::make_point( 0.0, 1.6), + PointType::make_point( 0.6, 1.2) + }; + + std::vector orders = {2,2,1}; + CurvedPolygonType bPolygon = createPolygon(CP,orders); + + CoordType trueA = -.0906666666666666666666; + PointType trueC= PointType::make_point(.2970147058823527,1.55764705882353); + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, moments_quad_all_orders) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test checking CurvedPolygon linear triangle moment computation."); + std::vector CPorig = { + PointType::make_point( 0.0, 0.0), + PointType::make_point( 0.0, 1.0), + PointType::make_point( 1.0, 1.0), + PointType::make_point( 1.0, 0.0), + PointType::make_point( 0.0, 0.0) + }; + + std::vector orders = {1,1,1,1}; + CurvedPolygonType bPolygon = createPolygon(CPorig,orders); + + CoordType trueA= 1.0; + PointType trueC= PointType::make_point(0.5,0.5); + + checkMoments(bPolygon, trueA, trueC, 1e-14, 1e-15); + for (int p=2; p<11; ++p) + { + std::vector CP = CPorig; + for (int side=0; side<4; ++side) + { + for (int i=1; i; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test intersecting two linear triangular CurvedPolygons (single region)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.3, 2.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0 , 1.6) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point( 0.0 , 1.6), + PointType::make_point( 0.6, 1.2) + }; + + BezierCurveType bCurve(controlPoints,order); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + bPolygon.addEdge(bCurve3); + CurvedPolygonType bPolygon2=bPolygon; + for (int i=0; i bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3); + EXPECT_TRUE(didIntersect); + std::cout << bPolygons3.size() << std::endl; + // int nEd= bPolygons3.size(); + for (int i=0; i< static_cast(bPolygons3.size()); ++i) + { + std::cout << bPolygons3[i] << std::endl; + } + for (int i=0; i(bPolygons3.size()); ++idxcurve) + { + for (int k=0; k; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (single region)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.4, 1.3), + PointType::make_point(0.3, 2.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.27, 1.5), + PointType::make_point( 0.0 , 1.6) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point( 0.0 , 1.6), + PointType::make_point( 0.1 , 1.5), + PointType::make_point( 0.6, 1.2) + }; + + BezierCurveType bCurve(controlPoints,order); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + bPolygon.addEdge(bCurve3); + CurvedPolygonType bPolygon2=bPolygon; + for (int i=0; i bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3); + EXPECT_TRUE(didIntersect); + + for (int i=0; i(bPolygons3.size()); ++idxcurve) + { + for (int k=0; k; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test finding area of intersection two linear triangular CurvedPolygons (single region)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.3, 2.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.0 , 1.6) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point( 0.0 , 1.6), + PointType::make_point( 0.6, 1.2) + }; + + BezierCurveType bCurve(controlPoints,order); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + bPolygon.addEdge(bCurve3); + CurvedPolygonType bPolygon2=bPolygon; + for (int i=0; i bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3); + EXPECT_TRUE(didIntersect); + + + CoordType A=0.0; + for (int i=0; i(bPolygons3.size()); ++i) + { + A+=bPolygons3[i].area(1e-14); + } + CoordType expA=-0.0793347222222222222; + EXPECT_NEAR(A,expA,1e-10); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, area_intersection_triangle_quadratic) +{ + const int DIM = 2; + const int order = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (single region)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.4, 1.3), + PointType::make_point(0.3, 2.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.27, 1.5), + PointType::make_point( 0.0 , 1.6) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point( 0.0 , 1.6), + PointType::make_point( 0.1 , 1.5), + PointType::make_point( 0.6, 1.2) + }; + + BezierCurveType bCurve(controlPoints,order); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + bPolygon.addEdge(bCurve3); + CurvedPolygonType bPolygon2=bPolygon; + for (int i=0; i bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3,1e-14); + EXPECT_TRUE(didIntersect); + + CoordType A=0.0; + for (int i=0; i(bPolygons3.size()); ++i) + { + A+=bPolygons3[i].area(1e-8); + } + CoordType expA=-0.024649833203616; + EXPECT_NEAR(A,expA,1e-10); +} + +TEST( primal_curvedpolygon, area_intersection_triangle_quadratic_two_regions) +{ + const int DIM = 2; + const int order = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (two regions)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.6 , 1.2), + PointType::make_point(0.4 , 1.3), + PointType::make_point(0.3 , 2.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point(0.3 , 2.0), + PointType::make_point(0.27, 1.5), + PointType::make_point(0.0 , 1.6) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point(0.0 , 1.6), + PointType::make_point(0.1 , 1.5), + PointType::make_point(0.6 , 1.2) + }; + + PointType controlPoints4[order+1] = { + PointType::make_point(1.0205, 1.6699), + PointType::make_point(0.8339, 1.5467), + PointType::make_point(0.1777, 1.8101) + }; + + PointType controlPoints5[order+1] = { + PointType::make_point(0.1777, 1.8101), + PointType::make_point(0.5957, 1.5341), + PointType::make_point(0.3741, 1.3503) + }; + + PointType controlPoints6[order+1] = { + PointType::make_point(0.3741, 1.3503), + PointType::make_point(0.5107, 1.3869), + PointType::make_point(1.0205, 1.6699) + }; + CurvedPolygonType bPolygon2; + + BezierCurveType bCurve(controlPoints,order); + BezierCurveType bCurve4(controlPoints4,order); + bPolygon2.addEdge(bCurve4); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + BezierCurveType bCurve5(controlPoints5,order); + bPolygon2.addEdge(bCurve5); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + BezierCurveType bCurve6(controlPoints6,order); + bPolygon2.addEdge(bCurve6); + bPolygon.addEdge(bCurve3); + + std::vector bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3); + EXPECT_TRUE(didIntersect); + EXPECT_EQ(bPolygons3.size(),2); +} + +TEST( primal_curvedpolygon, area_intersection_triangle_inclusion) +{ + const int DIM = 2; + const int order = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + using BezierCurveType = primal::BezierCurve< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (inclusion)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + PointType controlPoints[order+1] = { + PointType::make_point(0.0 , 0.0), + PointType::make_point(0.5 , 0.0), + PointType::make_point(1.0 , 0.0) + }; + + PointType controlPoints2[order+1] = { + PointType::make_point(1.0 , 0.0), + PointType::make_point(0.5, 0.5), + PointType::make_point(0.0 , 1.0) + }; + + PointType controlPoints3[order+1] = { + PointType::make_point(0.0 , 1.0), + PointType::make_point(0.0 , 0.5), + PointType::make_point(0.0 , 0.0) + }; + + BezierCurveType bCurve(controlPoints,order); + bPolygon.addEdge(bCurve); + + BezierCurveType bCurve2(controlPoints2,order); + bPolygon.addEdge(bCurve2); + + BezierCurveType bCurve3(controlPoints3,order); + bPolygon.addEdge(bCurve3); + + CurvedPolygonType bPolygon2=bPolygon; + for (int i=0; i bPolygons3; + bool didIntersect=intersect(bPolygon,bPolygon2,bPolygons3); + bPolygons3.clear(); + bool didIntersect2=intersect(bPolygon2,bPolygon,bPolygons3); + EXPECT_TRUE(didIntersect); + EXPECT_TRUE(didIntersect2); + EXPECT_EQ(bPolygons3.size(),1); +} +//---------------------------------------------------------------------------------- +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + + return result; +} 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..6eb12e93c6 --- /dev/null +++ b/src/axom/primal/tests/primal_curved_polygon_intersect.cpp @@ -0,0 +1,300 @@ +// Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/* /file bezier_test.cpp + * /brief This file tests the BezierCurve.hpp and eval_bezier.hpp files + */ + +#include "gtest/gtest.h" + +#include "axom/primal/geometry/CurvedPolygon.hpp" +#include "axom/primal/operators/intersect.hpp" + +#include + +namespace primal = axom::primal; + +/** + * 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) +{ + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + + 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 + EXPECT_EQ(expbPolygon.size(),intersectionPolys.size()); + + //Check that expected intersection curves are found to within test_eps + for (int i=0; i +primal::CurvedPolygon createPolygon( + const std::vector> ControlPoints, + const std::vector orders) +{ + using PointType = primal::Point< CoordType, DIM>; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + 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 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, intersection_triangle_linear) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test intersecting two linear triangular CurvedPolygons (single region)."); + + std::vector CP = { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.3, 2.0), + PointType::make_point( 0.0 , 1.6), + PointType::make_point(0.6, 1.2) + }; + std::vector orders = {1,1,1}; + + CurvedPolygonType bPolygon1 = createPolygon(CP,orders); + + std::vector CP2 = { + PointType::make_point(0.71, 1.31), + PointType::make_point(0.41, 2.11), + PointType::make_point( 0.11 , 1.71), + PointType::make_point(0.71, 1.31) + }; + + + CurvedPolygonType bPolygon2 = createPolygon(CP2,orders); + + std::vector expCP = { + PointType::make_point(0.3091666666666666666666, 1.9755555555555555555), + PointType::make_point(0.11, 1.71), + PointType::make_point( 0.5083333333333333333 , 1.44444444444444444444), + PointType::make_point(0.3091666666666666666666, 1.9755555555555555555) + }; + std::vector exporders={1,1,1}; + CurvedPolygonType expbPolygon = createPolygon(expCP,exporders); + + std::vector expbPolygons {expbPolygon}; + checkIntersection(bPolygon1,bPolygon2,expbPolygons); +} + +//---------------------------------------------------------------------------------- +TEST( primal_curvedpolygon, intersection_triangle_quadratic) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (single region)."); + + CurvedPolygonType bPolygon; + EXPECT_EQ(0, bPolygon.numEdges()); + + std::vector CP= { + PointType::make_point(0.6, 1.2), + PointType::make_point(0.4, 1.3), + PointType::make_point( 0.3, 2.0), + PointType::make_point( 0.27, 1.5), + PointType::make_point( 0.0 , 1.6), + PointType::make_point( 0.1 , 1.5), + PointType::make_point( 0.6, 1.2) + }; + std::vector orders = {2,2,2}; + CurvedPolygonType bPolygon1 = createPolygon(CP,orders); + + std::vector CP2= { + PointType::make_point(0.71, 1.31), + PointType::make_point(0.51, 1.41), + PointType::make_point( 0.41, 2.11), + PointType::make_point( 0.38, 1.61), + PointType::make_point( 0.11 , 1.71), + PointType::make_point( 0.21 , 1.61), + PointType::make_point( 0.71, 1.31) + }; + CurvedPolygonType bPolygon2 = createPolygon(CP2,orders); + + std::vector expCP = { + PointType::make_point(0.335956890729522, 1.784126953773395), + PointType::make_point(0.297344765794753, 1.718171485335525), + PointType::make_point( 0.2395677533016981, 1.700128235793371), + PointType::make_point( 0.221884203146682, 1.662410644580941 ), + PointType::make_point( 0.199328465398189, 1.636873522352205 ), + PointType::make_point(0.277429214338182,1.579562422716502 ), + PointType::make_point(0.408882616650578,1.495574996394597 ), + PointType::make_point(0.368520120719339,1.616453177259694), + PointType::make_point(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< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (two regions)."); + + std::vector CP1 = { + PointType::make_point(0.6 , 1.2), + PointType::make_point(0.4 , 1.3), + PointType::make_point(0.3 , 2.0), + PointType::make_point(0.27, 1.5), + PointType::make_point(0.0 , 1.6), + PointType::make_point(0.1 , 1.5), + PointType::make_point(0.6 , 1.2) + }; + + std::vector CP2 = { + PointType::make_point(1.0205, 1.6699), + PointType::make_point(0.8339, 1.5467), + PointType::make_point(0.1777, 1.8101), + PointType::make_point(0.5957, 1.5341), + PointType::make_point(0.3741, 1.3503), + PointType::make_point(0.5107, 1.3869), + PointType::make_point(1.0205, 1.6699) + }; + std::vector orders = {2,2,2}; + + std::vector expCP1 = { + PointType::make_point(0.343364196589264, 1.747080669655736), + PointType::make_point(0.305984025190458, 1.760433098612141), + PointType::make_point(0.266743999290327, 1.775316659915674), + PointType::make_point(0.263419346128088, 1.763343410502168), + PointType::make_point(0.259796003065908, 1.752116885838515), + PointType::make_point(0.320641367919239, 1.705796408318085), + PointType::make_point(0.362111919147859, 1.662268860466508), + PointType::make_point(0.352450139541348, 1.702947255097842), + PointType::make_point(0.343364196589264, 1.747080669655736), + }; + + std::vector expCP2 = { + PointType::make_point(0.454478985809487, 1.379250566393211 ), + PointType::make_point(0.444689566319939, 1.400290430035245 ), + PointType::make_point(0.435276730907216, 1.423589798138227 ), + PointType::make_point(0.416268597450954, 1.385275578571685 ), + PointType::make_point(0.374100000000000, 1.350300000000000 ), + PointType::make_point(0.404839872482010, 1.358536305511285 ), + PointType::make_point(0.454478985809487, 1.379250566393211 ) + }; + + 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); +} + +TEST( primal_curvedpolygon, area_intersection_triangle_inclusion) +{ + const int DIM = 2; + using CoordType = double; + using CurvedPolygonType = primal::CurvedPolygon< CoordType, DIM>; + using PointType = primal::Point< CoordType, DIM >; + + SLIC_INFO("Test intersecting two quadratic triangular CurvedPolygons (inclusion)."); + + std::vector CP1 = { + PointType::make_point(0.0 , 0.0), + PointType::make_point(0.5 , 0.0), + PointType::make_point(1.0 , 0.0), + PointType::make_point(0.5, 0.5), + PointType::make_point(0.0 , 1.0), + PointType::make_point(0.0 , 0.5), + PointType::make_point(0.0 , 0.0) + }; + + std::vector CP2 = { + PointType::make_point(0.05 , 0.05), + PointType::make_point(0.30 , 0.05), + PointType::make_point(0.55 , 0.05), + PointType::make_point(0.30, 0.30), + PointType::make_point(0.05 , 0.55), + PointType::make_point(0.05 , 0.30), + PointType::make_point(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); +} +//---------------------------------------------------------------------------------- +int main(int argc, char* argv[]) +{ + int result = 0; + + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::UnitTestLogger logger; // create & initialize test logger, + + result = RUN_ALL_TESTS(); + + return result; +} diff --git a/src/axom/quest/examples/CMakeLists.txt b/src/axom/quest/examples/CMakeLists.txt index 63663c5f32..aba90ec95d 100644 --- a/src/axom/quest/examples/CMakeLists.txt +++ b/src/axom/quest/examples/CMakeLists.txt @@ -14,9 +14,9 @@ set(quest_example_depends quest ) -blt_list_append(TO quest_example_depends ELEMENTS openmp IF ENABLE_OPENMP) +blt_list_append(TO quest_example_depends ELEMENTS openmp IF ${ENABLE_OPENMP}) -blt_list_append(TO quest_example_depends ELEMENTS cuda IF ENABLE_CUDA) +blt_list_append(TO quest_example_depends ELEMENTS cuda IF ${ENABLE_CUDA}) blt_add_executable( NAME quest_containment_driver_ex @@ -110,3 +110,22 @@ 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 fmt + 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() \ No newline at end of file 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..109b0abb32 --- /dev/null +++ b/src/axom/quest/examples/quest_high_order_remap.cpp @@ -0,0 +1,720 @@ +// Copyright (c) 2017-2019, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level COPYRIGHT 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 + */ + +#include + +// Axom includes +#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 "fmt/fmt.hpp" + +#include + +//namespace mint = axom::mint; +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: + using BBox = primal::BoundingBox; + using Point = 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_ASSERT(mesh != nullptr); + m_mesh = mesh; + + bool isHighOrder = (m_mesh->GetNodalFESpace() != nullptr) + && (m_mesh->GetNE() > 0); + SLIC_ASSERT_MSG(isHighOrder, "The mesh must be high order."); + + bool isBernstein = isBernsteinBasis(); + SLIC_ASSERT_MSG(isBernstein, "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); + + //SLIC_INFO("Elem " << elemId + // << " edge " << edgeIds[e] << " w/ orient " << edgeOrients[e] + // << " -- dof inds " + // << dofIndices[0] << " " << dofIndices[1] << " " << dofIndices[2] + // << " -- points " + // << spacePointFromDof(dofIndices[0], fes, nodes) << " " + // << spacePointFromDof(dofIndices[1], fes, nodes) << " " + // << spacePointFromDof(dofIndices[2], fes, nodes) + // ); + + // 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 */ + Point spacePointFromDof(int idx, + const mfem::FiniteElementSpace* fes, + const mfem::GridFunction* nodes) + { + return Point::make_point( + (*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 +{ +private: + using GridType = spin::ImplicitGrid<2, int>; +public: + using CandidateList = std::vector; + +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(); + if (nodal_fe_space == nullptr) { std::cerr << "project_to_pos_basis(): nodal_fe_space is NULL!" << std::endl; } + + const mfem::FiniteElementCollection *nodal_fe_coll = nodal_fe_space->FEColl(); + if (nodal_fe_coll == nullptr) { std::cerr << "project_to_pos_basis(): nodal_fe_coll is NULL!" << std::endl; } + + mfem::Mesh *gf_mesh = nodal_fe_space->GetMesh(); + if (gf_mesh == nullptr) { std::cerr << "project_to_pos_basis(): gf_mesh is NULL!" << std::endl; } + + int order = nodal_fe_space->GetOrder(0); + int dim = gf_mesh->Dimension(); + 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); + + auto* pos_fe_coll = new mfem::H1_FECollection(5, dim, + mfem::BasisType::Positive); + // 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_ASSERT_MSG( + // pos_fe_coll != AXOM_NULLPTR, + // "Problem generating a positive finite element collection " + // << "corresponding to the mesh's '"<< nodal_fe_coll->Name() + // << "' finite element collection."); + + if(pos_fe_coll != nullptr) + { + //DEBUG + //std::cerr << "Good so far... pos_fe_coll is not null. Making FESpace and GridFunction." << std::endl; + 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; + } + //DEBUG + else std::cerr << "BAD... pos_fe_coll is NULL. Could not make FESpace or GridFunction." << std::endl; + + //DEBUG + if (!out_pos_gf) + { + std::cerr << "project_to_pos_basis(): Construction failed; out_pos_gf is NULL!" << std::endl; + } + + // } + + return out_pos_gf; + } + + + /*! Set up the source and target meshes */ + void setupMeshes(int res1, int res2, int order) + { + const auto quadType = mfem::Element::QUADRILATERAL; + const int dim = 2; + + // paramters for source mesh -- quad mesh covering unit square + const int src_res = res2; + const int src_ord = order; + + // paramters for target mesh -- quad mesh covering (part of) unit square + const int tgt_res = res1; + const int tgt_ord = order; + const double tgt_scale = .712378102150; + const double tgt_trans1 = .1345747181586; + const double tgt_trans2 = .1345747181586; + + // create the source mesh + { + // create mfem mesh + auto* mesh = new mfem::Mesh(src_res, src_res, quadType, true); + + // create finite element collection for nodes + auto* fec = new mfem::H1_FECollection(src_ord, dim, + mfem::BasisType::Positive); + fecMap.Register("src_fec", fec, true); + + // create finite element space for nodes + auto* fes = new mfem::FiniteElementSpace(mesh, fec, dim); + fesMap.Register("src_fes", fes, true); + mesh->SetNodalFESpace(fes); + + // SLIC_INFO("Writing to: " << axom::utilities::filesystem::getCWD()); + { + std::ofstream file; + file.open("source_mesh.mfem"); + mesh->Print(file); + } + + srcMesh.setMesh(mesh); + } + + // create the target mesh + { + auto* mesh = new mfem::Mesh(tgt_res, tgt_res, quadType, true); + xformMesh(mesh, tgt_scale, tgt_trans1, tgt_trans2); + + auto* fec = new mfem::H1_FECollection(tgt_ord, dim, + mfem::BasisType::Positive); + fecMap.Register("tgt_fec", fec, true); + + auto* fes = new mfem::FiniteElementSpace(mesh, fec, dim); + fesMap.Register("tgt_fes", fes, true); + mesh->SetNodalFESpace(fes); + + { + std::ofstream file; + file.open("target_mesh.mfem"); + mesh->Print(file); + } + + tgtMesh.setMesh(mesh); + } + } + void loadMeshes(int res2, int order) + { + // create the source mesh + const auto quadType = mfem::Element::QUADRILATERAL; + const int dim = 2; + + // paramters for target mesh -- quad mesh covering unit square + const int src_res = res2; + const int src_ord = order; + const double src_scale = 1.5; + const double src_trans1 = .01517288412347; + const double src_trans2 = .02571238506182; + + // paramters for target mesh -- quad mesh covering (part of) unit square + const int tgt_ord = 2; + + { + // NOTE (KW): For now, assume we have AXOM_DATA_DIR + namespace fs = axom::utilities::filesystem; + std::string fname = fs::joinPath(AXOM_DATA_DIR, "mfem/disc-nurbs-80.mesh"); + + auto* mesh = new mfem::Mesh(fname.c_str(),1,1); + xformMesh(mesh, src_scale, src_trans1, src_trans2); + if (mesh->NURBSext) + { + int order =src_ord; + mesh->SetCurvature(5); + } + // xformMesh(mesh, tgt_scale, tgt_trans); + { + std::ofstream file; + file.open("src_mesh_orig.mfem"); + mesh->Print(file); + } + + bool is_mesh_gf_new; + 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); + + //auto* fec = new mfem::H1_FECollection(tgt_ord, dim, + // mfem::BasisType::Positive); + //fecMap.Register("tgt_fec", fec, true); + + //auto* fes = new mfem::FiniteElementSpace(mesh, fec, dim); + //fesMap.Register("tgt_fes", fes, true); + //mesh->SetNodalFESpace(fes); + std::cout << mesh->GetNV() << std::endl; + { + std::ofstream file; + file.open("src_mesh_set.mfem"); + mesh->Print(file); + } + //std::cout << "Got here!" << std::endl; + srcMesh.setMesh(mesh); + } + + // create the target mesh + + { + // NOTE (KW): For now, assume we have AXOM_DATA_DIR + namespace fs = axom::utilities::filesystem; + std::string fname = fs::joinPath(AXOM_DATA_DIR, "mfem/disc-nurbs-80.mesh"); + + auto* mesh = new mfem::Mesh(fname.c_str(),1,1); + if (mesh->NURBSext) + { + int order =tgt_ord; + mesh->SetCurvature(5); + } + + // xformMesh(mesh, tgt_scale, tgt_trans); + const double tgt_scale = 1.0; + const double tgt_trans1 = .001237586; + const double tgt_trans2 = -.06172376; + xformMesh(mesh, tgt_scale, tgt_trans1, tgt_trans2); + + { + std::ofstream file; + file.open("target_mesh_orig.mfem"); + mesh->Print(file); + } + + bool is_mesh_gf_new; + 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); + + //auto* fec = new mfem::H1_FECollection(tgt_ord, dim, + // mfem::BasisType::Positive); + //fecMap.Register("tgt_fec", fec, true); + + //auto* fes = new mfem::FiniteElementSpace(mesh, fec, dim); + //fesMap.Register("tgt_fes", fes, true); + //mesh->SetNodalFESpace(fes); + + { + std::ofstream file; + file.open("target_mesh_set.mfem"); + mesh->Print(file); + } + std::cout << "Got here!" << std::endl; + tgtMesh.setMesh(mesh); + } + } + /*! Setup the implicit grid spatial index over the source mesh */ + void setupGrid() + { + 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(); + // SLIC_INFO("Number of Target Elements: " << nTargetElems); + 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()) + break; + + auto tgtPoly = tgtMesh.elemAsCurvedPolygon(i); + //SLIC_INFO("Target Element: " << tgtPoly); + correctArea+=tgtPoly.area(); + //SLIC_INFO("Target elem " << i + // << " -- area " << tgtPoly.area() + // //<< " -- bbox " << tgtMesh.elementBoundingBox(i) + // ); + + double A=0.0; + for (int srcElem : candidates) + { + auto srcPoly = srcMesh.elemAsCurvedPolygon(srcElem); + // SLIC_INFO("*Source Element: " << srcPoly); + //SLIC_INFO("* Source elem " << srcElem + // << " -- area " << srcPoly.area() + // //// //<< " -- 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-=pnew[i].area(); + //SLIC_INFO("** Intersection area :" << -pnew[i].area() + // ); + } + } + srcPoly.reverseOrientation(); + tgtPoly.reverseOrientation(); + //SLIC_INFO("* Calculated area: " << srcElem + // << " -- area " << A + //// //<< " -- bbox " << srcMesh.elementBoundingBox(srcElem) + // ); + } + calcE+=abs(tgtPoly.area()-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; + } + +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; + } + } +}; + + + + + + +//------------------------------------------------------------------------------ +int main( int argc, char** argv ) +{ + axom::slic::UnitTestLogger logger; // create & initialize logger + + SLIC_INFO("The application conservatively maps fields from a source \n" + <<"high order mesh to a target high order mesh!"); +// int res1=1; +// Remapper remap; +// res1=5; +// int res2=res1+1; +// remap.loadMeshes(res1,2); +// remap.setupGrid(); +// double Area = remap.computeOverlapAreas(); +// std::cout << std::endl << Area << std::endl; + +// return 0; + + int res1=1; + std::cout << "["; + for (int i=2; i<=2; ++i) + { + Remapper remap; + + std::cout.precision(16); + // res1= res1*2; + + res1=i; + int res2= res1+2; + // Setup the two meshes in the Bernstein basis + // The current implementation hard-codes the two meshes + // TODO: Read in two meshes from disk. + // In that case, we will need to convert the FEC to the Bernstein basis + // remap.setupMeshes(res1, res2, i); + remap.loadMeshes(res1, i); + + + // Set up the spatial index + remap.setupGrid(); + auto start = std::chrono::high_resolution_clock::now(); + // Computes the overlaps between elements of the target and source meshes + double Area =remap.computeOverlapAreas(); + auto finish = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = finish-start; + std::cout << i << ", " << Area << "," << elapsed.count() << std::endl; + } + std::cout << "]" << std::endl; + return 0; +} diff --git a/src/axom/quest/tests/CMakeLists.txt b/src/axom/quest/tests/CMakeLists.txt index 496c3089ab..81f4a58364 100644 --- a/src/axom/quest/tests/CMakeLists.txt +++ b/src/axom/quest/tests/CMakeLists.txt @@ -56,6 +56,7 @@ if(MFEM_FOUND) SOURCES ${test_name}.cpp OUTPUT_DIR ${TEST_OUTPUT_DIRECTORY} DEPENDS_ON ${quest_tests_depends} mfem core + FOLDER axom/quest/tests ) string(STRIP "${AXOM_DISABLE_UNUSED_PARAMETER_WARNINGS}" MFEM_COMPILE_FLAGS)