Merge pull request #1161 from LLNL/feature/spainhour/bezier_patch_der…

…ivatives

Add additional derivatives to Bezier curve/patches

src/axom/primal/operators/detail/winding_number_impl.hpp

@@ -311,30 +311,56 @@ double stokes_winding_number(const Point<T, 3>& query,
 
     // Compute one of three vector field line integrals depending on
     //  the orientation of the original surface, indicated through ax.
-    if(ax == SingularityAxis::x)
+    switch(ax)
     {
+    case(SingularityAxis::x):
       quadrature += quad_rule.IntPoint(q).weight *
         (node[2] * node[0] * node_dt[1] - node[1] * node[0] * node_dt[2]) /
         (node[1] * node[1] + node[2] * node[2]) / node_norm;
-    }
-    else if(ax == SingularityAxis::y)
-    {
+      break;
+    case(SingularityAxis::y):
       quadrature += quad_rule.IntPoint(q).weight *
         (node[0] * node[1] * node_dt[2] - node[2] * node[1] * node_dt[0]) /
         (node[0] * node[0] + node[2] * node[2]) / node_norm;
-    }
-    else  // ax == SingularityAxis::z || ax == SingularityAxis::rotated
-    {
+      break;
+    case(SingularityAxis::z):
+    case(SingularityAxis::rotated):
       quadrature += quad_rule.IntPoint(q).weight *
         (node[1] * node[2] * node_dt[0] - node[0] * node[2] * node_dt[1]) /
         (node[0] * node[0] + node[1] * node[1]) / node_norm;
+      break;
     }
   }
 
-  // Adaptively refine quadrature over curves if query is not far enough away.
+  // Adaptively refine quadrature over curves if query is not far enough away
+  //  from the singularity axis. If rotated, assume you need to adapt.
+  bool needs_adapt = false;
   BoundingBox<T, 3> cBox(curve.boundingBox());
-  if(squared_distance(query, cBox.getCentroid()) <=
-     2 * cBox.range().squared_norm())
+  Point<T, 3> centroid = cBox.getCentroid();
+
+  switch(ax)
+  {
+  case(SingularityAxis::x):
+    needs_adapt = (query[1] - centroid[1]) * (query[1] - centroid[1]) +
+        (query[2] - centroid[2]) * (query[2] - centroid[2]) <=
+      cBox.range().squared_norm();
+    break;
+  case(SingularityAxis::y):
+    needs_adapt = (query[0] - centroid[0]) * (query[0] - centroid[0]) +
+        (query[2] - centroid[2]) * (query[2] - centroid[2]) <=
+      cBox.range().squared_norm();
+    break;
+  case(SingularityAxis::z):
+    needs_adapt = (query[0] - centroid[0]) * (query[0] - centroid[0]) *
+        (query[1] - centroid[1]) * (query[1] - centroid[1]) <=
+      cBox.range().squared_norm();
+    break;
+  case(SingularityAxis::rotated):
+    needs_adapt = true;
+    break;
+  }
+
+  if(needs_adapt)
   {
     return stokes_winding_number_adaptive(query,
                                           curve,
@@ -358,6 +384,7 @@ double stokes_winding_number(const Point<T, 3>& query,
  * \param [in] quad_rule The mfem quadrature rule object
  * \param [in] quad_coarse The integral evaluated on the original curve
  * \param [in] quad_tol The maximum relative error allowed in each quadrature
+ * \param [in] depth The current recursive depth
  * 
  * Recursively apply quadrature for one of three possible integrals along two halfs 
  * of a curve. The sum of this integral along the subcurves should be equal to to
@@ -394,23 +421,24 @@ double stokes_winding_number_adaptive(const Point<T, 3>& query,
 
       // Compute one of three vector field line integrals depending on
       //  the orientation of the original surface, indicated through ax.
-      if(ax == SingularityAxis::x)
+      switch(ax)
       {
+      case(SingularityAxis::x):
         quad_fine[i] += quad_rule.IntPoint(q).weight *
           (node[2] * node[0] * node_dt[1] - node[1] * node[0] * node_dt[2]) /
           (node[1] * node[1] + node[2] * node[2]) / node_norm;
-      }
-      else if(ax == SingularityAxis::y)
-      {
+        break;
+      case(SingularityAxis::y):
         quad_fine[i] += quad_rule.IntPoint(q).weight *
           (node[0] * node[1] * node_dt[2] - node[2] * node[1] * node_dt[0]) /
           (node[0] * node[0] + node[2] * node[2]) / node_norm;
-      }
-      else  // ax == SingularityAxis::z || ax == SingularityAxis::rotated
-      {
+        break;
+      case(SingularityAxis::z):
+      case(SingularityAxis::rotated):
         quad_fine[i] += quad_rule.IntPoint(q).weight *
           (node[1] * node[2] * node_dt[0] - node[0] * node[2] * node_dt[1]) /
           (node[0] * node[0] + node[1] * node[1]) / node_norm;
+        break;
       }
     }
   }
@@ -420,7 +448,7 @@ double stokes_winding_number_adaptive(const Point<T, 3>& query,
      axom::utilities::isNearlyEqualRelative(quad_fine[0] + quad_fine[1],
                                             quad_coarse,
                                             quad_tol,
-                                            0.0))
+                                            1e-10))
   {
     return 0.25 * M_1_PI * (quad_fine[0] + quad_fine[1]);
   }

src/axom/primal/operators/is_convex.hpp

@@ -9,6 +9,7 @@
 #include "axom/primal/geometry/Segment.hpp"
 #include "axom/primal/geometry/Polygon.hpp"
 #include "axom/primal/operators/orientation.hpp"
+#include "axom/primal/operators/squared_distance.hpp"
 
 namespace axom
 {
@@ -19,9 +20,10 @@ namespace primal
  * 
  * \param [in] poly The polygon
  * 
- * Uses dot products to detect whether vertices extend in the "convex" direction.
- * Uses the edge P[0]P[N] as a reference, meaning its adjacent edges are
- * always considered to be oriented correctly.
+ * A 2D polygon is convex when every line that does not contain an edge
+ * intersects the shape at most twice.
+ * Checks whether for each pair of vertices P[i-1]P[i+1], the point
+ * P[i] and (P[0] or P[n]) lie on the same side of the line connecting them.
  * 
  * Algorithm adapted from:
  * Y. L. Ma, W.T. Hewitt. "Point inversion and projection for NURBS curve 
@@ -54,7 +56,7 @@ bool is_convex(const Polygon<T, 2>& poly, double EPS = 1e-8)
     }
 
     // Ensure other point to check against isn't adjacent
-    if(res1 == orientation(poly[(i < n / 2) ? n : 0], seg, EPS))
+    if(res1 == orientation(poly[(i <= n / 2) ? n : 0], seg, EPS))
     {
       return false;
     }

src/axom/primal/operators/winding_number.hpp

@@ -454,7 +454,8 @@ int winding_number(const Point<T, 3>& query,
 }
 
 #ifdef AXOM_USE_MFEM
-/*!
+
+/*
  * \brief Computes the solid angle winding number for a Bezier patch
  *
  * \param [in] query The query point to test
@@ -463,6 +464,7 @@ int winding_number(const Point<T, 3>& query,
  *                      considered indistinguishable
  * \param [in] quad_tol The maximum relative error allowed in the quadrature
  * \param [in] EPS Miscellaneous numerical tolerance level for nonphysical distances
+ * \param [in] depth The current recursive depth
  * 
  * Computes the generalized winding number for a Bezier patch using Stokes theorem.
  *
@@ -473,7 +475,7 @@ int winding_number(const Point<T, 3>& query,
  */
 template <typename T>
 double winding_number(const Point<T, 3>& query,
-                      const BezierPatch<T>& bPatch,
+                      const BezierPatch<T, 3>& bPatch,
                       const double edge_tol = 1e-8,
                       const double quad_tol = 1e-8,
                       const double EPS = 1e-8,
@@ -507,7 +509,7 @@ double winding_number(const Point<T, 3>& query,
   constexpr double edge_offset = 0.01;
   if(squared_distance(query, bPatch(0, 0)) <= edge_tol_sq)
   {
-    BezierPatch<T> p1, p2, p3, p4;
+    BezierPatch<T, 3> p1, p2, p3, p4;
     bPatch.split(0.0 + edge_offset, 0.0 + edge_offset, p1, p2, p3, p4);
     double new_edge_tol = 0.5 *
       sqrt(axom::utilities::min(
@@ -521,7 +523,7 @@ double winding_number(const Point<T, 3>& query,
   }
   if(squared_distance(query, bPatch(ord_u, 0)) <= edge_tol_sq)
   {
-    BezierPatch<T> p1, p2, p3, p4;
+    BezierPatch<T, 3> p1, p2, p3, p4;
     bPatch.split(1.0 - edge_offset, 0.0 + edge_offset, p1, p2, p3, p4);
     double new_edge_tol = 0.5 *
       sqrt(axom::utilities::min(
@@ -535,7 +537,7 @@ double winding_number(const Point<T, 3>& query,
   }
   if(squared_distance(query, bPatch(0, ord_v)) <= edge_tol_sq)
   {
-    BezierPatch<T> p1, p2, p3, p4;
+    BezierPatch<T, 3> p1, p2, p3, p4;
     bPatch.split(0.0 + edge_offset, 1.0 - edge_offset, p1, p2, p3, p4);
     double new_edge_tol = 0.5 *
       sqrt(axom::utilities::min(
@@ -549,7 +551,7 @@ double winding_number(const Point<T, 3>& query,
   }
   if(squared_distance(query, bPatch(ord_u, ord_v)) <= edge_tol_sq)
   {
-    BezierPatch<T> p1, p2, p3, p4;
+    BezierPatch<T, 3> p1, p2, p3, p4;
     bPatch.split(1.0 - edge_offset, 1.0 - edge_offset, p1, p2, p3, p4);
     double new_edge_tol = 0.5 *
       sqrt(axom::utilities::min(
@@ -602,7 +604,7 @@ double winding_number(const Point<T, 3>& query,
     OrientedBoundingBox<T, 3> oBox(bPatch.orientedBoundingBox().expand(edge_tol));
     if(oBox.contains(query))
     {
-      BezierPatch<T> p1, p2, p3, p4;
+      BezierPatch<T, 3> p1, p2, p3, p4;
       bPatch.split(0.5, 0.5, p1, p2, p3, p4);
       return winding_number(query, p1, edge_tol, quad_tol, EPS, depth + 1) +
         winding_number(query, p2, edge_tol, quad_tol, EPS, depth + 1) +