Adds integral evaluation on Bezier curves/bounded regions

RELEASE-NOTES.md

@@ -90,8 +90,9 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
   in quest's signed_distance C API.
 - Adds utility function for linear interpolation (`lerp`) of two numbers
 - Adds utility function to compute binomial coefficients  
-- Adds a `CurvedPolygon` class to primal representing a polygon with `BezierCurves` as edges
+- Adds a `CurvedPolygon` class to primal representing a polygon with `BezierCurve`s as edges
 - Adds functions to compute the moments (area, centroid) of a `CurvedPolygon`
+- Adds functions to evaluate integrals over `BezierCurve` and `CurvedPolygon` objects 
 
 ###  Changed
 - Axom now requires C++14 and will default to that if not specified via `BLT_CXX_STD`.

src/axom/primal/CMakeLists.txt

@@ -68,6 +68,13 @@ set( primal_sources
      operators/detail/intersect_impl.cpp
    )
 
+blt_list_append(
+    TO       primal_headers 
+    ELEMENTS operators/evaluate_integral.hpp
+             operators/detail/evaluate_integral_impl.hpp
+    IF       MFEM_FOUND 
+               )
+
 #------------------------------------------------------------------------------
 # Build and install the library
 #------------------------------------------------------------------------------
@@ -78,6 +85,7 @@ set( primal_dependencies
 
 blt_list_append( TO primal_dependencies ELEMENTS cuda IF ENABLE_CUDA )
 blt_list_append( TO primal_dependencies ELEMENTS blt::hip_runtime IF ENABLE_HIP )
+blt_list_append( TO primal_dependencies ELEMENTS mfem IF MFEM_FOUND )
 
 blt_add_library(
     NAME                  primal

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

@@ -0,0 +1,150 @@
+// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and
+// other Axom Project Developers. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: (BSD-3-Clause)
+
+#ifndef PRIMAL_EVAL_INTEGRAL_IMPL_HPP_
+#define PRIMAL_EVAL_INTEGRAL_IMPL_HPP_
+
+// Axom includes
+#include "axom/config.hpp"  // for compile-time configuration options
+#include "axom/primal.hpp"
+
+// MFEM includes
+#ifdef AXOM_USE_MFEM
+  #include "mfem.hpp"
+#else
+  #error "Primal's integral evaluation functions require mfem library."
+#endif
+
+// C++ includes
+#include <cmath>
+
+namespace axom
+{
+namespace primal
+{
+namespace detail
+{
+/*!
+ * \brief Evaluate a scalar field line integral on a single Bezier curve.
+ *
+ * Evaluate the scalar field line integral with MFEM integration rule
+ *
+ * \param [in] c the Bezier curve object
+ * \param [in] integrand the lambda function representing the integrand. 
+ * Must accept a 2D point as input and return a double
+ * \param [in] quad the mfem integration rule containing nodes and weights
+ * \return the value of the integral
+ */
+inline double evaluate_line_integral_component(
+  const primal::BezierCurve<double, 2>& c,
+  std::function<double(Point2D)> integrand,
+  const mfem::IntegrationRule& quad)
+{
+  // Store/compute quadrature result
+  double full_quadrature = 0.0;
+  for(int q = 0; q < quad.GetNPoints(); q++)
+  {
+    // Get intermediate quadrature point
+    //  at which to evaluate tangent vector
+    auto x_q = c.evaluate(quad.IntPoint(q).x);
+    auto dx_q = c.dt(quad.IntPoint(q).x);
+
+    full_quadrature += quad.IntPoint(q).weight * integrand(x_q) * dx_q.norm();
+  }
+
+  return full_quadrature;
+}
+
+/*!
+ * \brief Evaluate a vector field line integral on a single Bezier curve.
+ *
+ * Evaluate the vector field line integral with MFEM integration rule
+ *
+ * \param [in] c the Bezier curve object
+ * \param [in] integrand the lambda function representing the integrand. 
+ * Must accept a 2D point as input and return a 2D axom vector object
+ * \param [in] quad the mfem integration rule containing nodes and weights
+ * \return the value of the integral
+ */
+inline double evaluate_line_integral_component(
+  const primal::BezierCurve<double, 2>& c,
+  std::function<Vector2D(Point2D)> vec_field,
+  const mfem::IntegrationRule& quad)
+{
+  // Store/compute quadrature result
+  double full_quadrature = 0.0;
+  for(int q = 0; q < quad.GetNPoints(); q++)
+  {
+    // Get intermediate quadrature point
+    //  on which to evaluate dot product
+    auto x_q = c.evaluate(quad.IntPoint(q).x);
+    auto dx_q = c.dt(quad.IntPoint(q).x);
+    auto func_val = vec_field(x_q);
+
+    full_quadrature +=
+      quad.IntPoint(q).weight * Vector2D::dot_product(func_val, dx_q);
+  }
+
+  return full_quadrature;
+}
+
+/*!
+ * \brief Evaluate the area integral across one component of the curved polygon.
+ *
+ * Intended to be called for each BezierCurve object in a curved polygon.
+ * Uses a Spectral Mesh-Free Quadrature derived from Green's theorem, evaluating
+ * the area integral as a line integral of the antiderivative over the curve.
+ * For algorithm details, see "Spectral Mesh-Free Quadrature for Planar 
+ * Regions Bounded by Rational Parametric Curves" by David Gunderman et al.
+ * 
+ * \param [in] cs the array of Bezier curve objects that bound the region
+ * \param [in] integrand the lambda function representing the integrand. 
+ * Must accept a 2D point as input and return a double
+ * \param [in] The lower bound of integration for the antiderivatives
+ * \param [in] quad_Q the quadrature rule for the line integral
+ * \param [in] quad_P the quadrature rule for the antiderivative
+ * \return the value of the integral, which is mathematically meaningless.
+ */
+template <class Lambda>
+double evaluate_area_integral_component(const primal::BezierCurve<double, 2>& c,
+                                        Lambda&& integrand,
+                                        double int_lb,
+                                        const mfem::IntegrationRule& quad_Q,
+                                        const mfem::IntegrationRule& quad_P)
+{
+  // Store some intermediate values
+  double antiderivative = 0.0;
+
+  // Store/compute quadrature result
+  double full_quadrature = 0.0;
+  for(int q = 0; q < quad_Q.GetNPoints(); q++)
+  {
+    // Get intermediate quadrature point
+    //  on which to evaluate antiderivative
+    auto x_q = c.evaluate(quad_Q.IntPoint(q).x);
+
+    // Evaluate the antiderivative at x_q, add it to full quadrature
+    for(int xi = 0; xi < quad_P.GetNPoints(); xi++)
+    {
+      // Define interior quadrature points
+      auto x_qxi =
+        Point2D({x_q[0], (x_q[1] - int_lb) * quad_P.IntPoint(xi).x + int_lb});
+
+      antiderivative =
+        quad_P.IntPoint(xi).weight * (x_q[1] - int_lb) * integrand(x_qxi);
+
+      full_quadrature += quad_Q.IntPoint(q).weight *
+        c.dt(quad_Q.IntPoint(q).x)[0] * -antiderivative;
+    }
+  }
+
+  return full_quadrature;
+}
+
+}  // end namespace detail
+}  // end namespace primal
+}  // end namespace axom
+
+#endif