diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index 3849715299..acef35281e 100644 --- a/RELEASE-NOTES.md +++ b/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`. diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 920fd4a3b8..c917281e1c 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/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 diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp new file mode 100644 index 0000000000..c5934e6adc --- /dev/null +++ b/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 + +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& c, + std::function 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& c, + std::function 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 +double evaluate_area_integral_component(const primal::BezierCurve& 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 \ No newline at end of file diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp new file mode 100644 index 0000000000..12d22a8a94 --- /dev/null +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -0,0 +1,250 @@ +// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and +// other Axom Project Developers. See the top-level LICENSE file for details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +/*! + * \file evaluate_integral.hpp + * + * \brief Consists of methods that evaluate integrals over regions defined + * by Bezier curves, such as 2D area integrals and scalar/vector field line integrals + * + * Line integrals are computed with 1D quadrature rules supplied + * by MFEM. 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar + * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. + * + * \note This requires the MFEM third-party library + */ + +#ifndef PRIMAL_EVAL_INTEGRAL_HPP_ +#define PRIMAL_EVAL_INTEGRAL_HPP_ + +// Axom includes +#include "axom/config.hpp" +#include "axom/primal.hpp" + +#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" + +// MFEM includes +#ifdef AXOM_USE_MFEM + #include "mfem.hpp" +#else + #error "Primal's integral evaluation functions require mfem library." +#endif + +// C++ includes +#include + +namespace axom +{ +namespace primal +{ +/*! + * \brief Evaluate a line integral along a collection of Bezier curves + * + * The line integral is evaluated on each curve in the array, and added + * together to represent the total integral. The curves need not be connected. + * Uses 1D Gaussian quadrature generated by MFEM. + * + * \param [in] cs the array of Bezier curve objects + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input, and return a double (scalar field) or + * 2D vector object (vector field). + * \param [in] npts the number of Gaussian quadrature nodes for each component + * \return the value of the integral + */ +template +double evaluate_line_integral(const axom::Array>& cs, + Lambda&& integrand, + int npts) +{ + // Generate quadrature library, defaulting to GaussLegendre quadrature. + // Use the same one for every curve in the polygon + // Quadrature order is equal to 2*N - 1 + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + double total_integral = 0.0; + for(const auto& curve : cs) + { + // Compute the line integral along each component. + total_integral += + detail::evaluate_line_integral_component(curve, integrand, quad); + } + + return total_integral; +} + +/*! + * \brief Evaluate a line integral along the boundary of a CurvedPolygon object. + * + * See above definition for details. + * + * \param [in] cpoly the CurvedPolygon object + * \param [in] integrand the lambda function representing the integrand. + * Must accept a 2D point as input and return a double + * \param [in] npts_Q the number of quadrature points to evaluate the line integral + * \param [in] npts_P the number of quadrature points to evaluate the antiderivative + * \return the value of the integral + */ +template +double evaluate_line_integral(const primal::CurvedPolygon cpoly, + Lambda&& integrand, + int npts) +{ + // Generate quadrature library, defaulting to GaussLegendre quadrature. + // Use the same one for every curve in the polygon + // Quadrature order is equal to 2*N - 1 + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + double total_integral = 0.0; + for(int i = 0; i < cpoly.numEdges(); i++) + { + // Compute the line integral along each component. + total_integral += + detail::evaluate_line_integral_component(cpoly[i], integrand, quad); + } + + return total_integral; +} + +/*! + * \brief Evaluate a line integral on a single Bezier curve. + * + * Evaluate the line integral with a given number of Gaussian + * quadrature nodes generated by MFEM. + * + * \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 (scalar field) or + * 2D vector object (vector field). + * \param [in] npts the number of quadrature nodes + * \return the value of the integral + */ +template +double evaluate_line_integral(const primal::BezierCurve& c, + Lambda&& integrand, + int npts) +{ + // Generate quadrature library, defaulting to GaussLegendre quadrature. + // Use the same one for every curve in the polygon + // Quadrature order is equal to 2*N - 1 + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); + + return detail::evaluate_line_integral_component(c, integrand, quad); +} + +/*! + * \brief Evaluate an integral across a 2D domain bounded by Bezier curves. + * + * Assumes that the array of Bezier curves is closed and connected. Will compute + * the integral regardless, but the result will be meaningless. + * 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. + * + * \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] npts_Q the number of quadrature points to evaluate the line integral + * \param [in] npts_P the number of quadrature points to evaluate the antiderivative + * \return the value of the integral + */ +template +double evaluate_area_integral(const axom::Array>& cs, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + // Generate quadrature library, defaulting to GaussLegendre quadrature. + // Use the same one for every curve in the polygon + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + + if(npts_P <= 0) npts_P = npts_Q; + + // Get the quadrature for the line integral. + // Quadrature order is equal to 2*N - 1 + const mfem::IntegrationRule& quad_Q = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1); + const mfem::IntegrationRule& quad_P = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1); + + // Use minimum y-coord of control nodes as lower bound for integration + double int_lb = cs[0][0][1]; + for(int i = 0; i < cs.size(); i++) + for(int j = 1; j < cs[i].getOrder() + 1; j++) + int_lb = std::min(int_lb, cs[i][j][1]); + + // Evaluate the antiderivative line integral along each component + double total_integral = 0.0; + for(const auto& curve : cs) + { + total_integral += detail::evaluate_area_integral_component(curve, + integrand, + int_lb, + quad_Q, + quad_P); + } + + return total_integral; +} + +/*! + * \brief Evaluate an integral on the interior of a CurvedPolygon object. + * + * See above definition for details. + * + * \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] npts_Q the number of quadrature points to evaluate the line integral + * \param [in] npts_P the number of quadrature points to evaluate the antiderivative + * \return the value of the integral + */ +template +double evaluate_area_integral(const primal::CurvedPolygon cpoly, + Lambda&& integrand, + int npts_Q, + int npts_P = 0) +{ + // Generate quadrature library, defaulting to GaussLegendre quadrature. + // Use the same one for every curve in the polygon + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + + if(npts_P <= 0) npts_P = npts_Q; + + // Get the quadrature for the line integral. + // Quadrature order is equal to 2*N - 1 + const mfem::IntegrationRule& quad_Q = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1); + const mfem::IntegrationRule& quad_P = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1); + + // Use minimum y-coord of control nodes as lower bound for integration + double int_lb = cpoly[0][0][1]; + for(int i = 0; i < cpoly.numEdges(); i++) + for(int j = 1; j < cpoly[i].getOrder() + 1; j++) + int_lb = std::min(int_lb, cpoly[i][j][1]); + + // Evaluate the antiderivative line integral along each component + double total_integral = 0.0; + for(int i = 0; i < cpoly.numEdges(); i++) + { + total_integral += detail::evaluate_area_integral_component(cpoly[i], + integrand, + int_lb, + quad_Q, + quad_P); + } + + return total_integral; +} + +} // namespace primal +} // end namespace axom + +#endif diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index 6d246266bb..6b11631efd 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -38,6 +38,8 @@ set( primal_tests primal_zip.cpp ) +blt_list_append( TO primal_tests ELEMENTS primal_integral.cpp IF MFEM_FOUND ) + set(primal_test_depends axom fmt gtest) blt_list_append( TO primal_test_depends ELEMENTS cuda IF ENABLE_CUDA ) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp new file mode 100644 index 0000000000..cb6bd7924a --- /dev/null +++ b/src/axom/primal/tests/primal_integral.cpp @@ -0,0 +1,261 @@ +// 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) + +#include "axom/primal.hpp" +#include "axom/slic.hpp" +#include "axom/fmt.hpp" + +#include "gtest/gtest.h" + +#include + +namespace primal = axom::primal; + +TEST(primal_integral, evaluate_area_integral) +{ + using Point2D = primal::Point; + using Bezier = primal::BezierCurve; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 20; + + // Define anonymous functions for testing + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; + auto transc_integrand = [](Point2D x) -> double { + return std::sin(x[0] * x[1]); + }; + + // Test on triangular domain + Point2D trinodes1[] = {Point2D {0.0, 0.0}, Point2D {1.0, 0.0}}; + Bezier tri1(trinodes1, 1); + + Point2D trinodes2[] = {Point2D {1.0, 0.0}, Point2D {0.0, 1.0}}; + Bezier tri2(trinodes2, 1); + + Point2D trinodes3[] = {Point2D {0.0, 1.0}, Point2D {0.0, 0.0}}; + Bezier tri3(trinodes3, 1); + + axom::Array triangle({tri1, tri2, tri3}); + + // Compare against hand computed/high-precision calculated values + EXPECT_NEAR(evaluate_area_integral(triangle, const_integrand, npts), + 0.5, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(triangle, poly_integrand, npts), + 1.0 / 60.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(triangle, transc_integrand, npts), + 0.0415181074232, + abs_tol); + + // Test on parabolic domain (between f(x) = 1-x^2 and g(x) = x^2-1, shifted to the right 1 unit) + Point2D paranodes1[] = {Point2D {2.0, 0.0}, + Point2D {1.0, 2.0}, + Point2D {0.0, 0.0}}; + Bezier para1(paranodes1, 2); + + Point2D paranodes2[] = {Point2D {0.0, 0.0}, + Point2D {1.0, -2.0}, + Point2D {2.0, 0.0}}; + Bezier para2(paranodes2, 2); + + axom::Array parabola_shape({para1, para2}); + + // Compare against hand computed/high-precision calculated values + EXPECT_NEAR(evaluate_area_integral(parabola_shape, const_integrand, npts), + 8.0 / 3.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_shape, poly_integrand, npts), + 64.0 / 105.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_shape, transc_integrand, npts), + 0.0, + abs_tol); + + // Ensure compatibility with curved polygons + Bezier pedges[2] = {para1, para2}; + primal::CurvedPolygon parabola_polygon(pedges, 2); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, const_integrand, npts), + 8.0 / 3.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, poly_integrand, npts), + 64.0 / 105.0, + abs_tol); + EXPECT_NEAR(evaluate_area_integral(parabola_polygon, transc_integrand, npts), + 0.0, + abs_tol); +} + +TEST(primal_integral, evaluate_line_integral_scalar) +{ + using Point2D = primal::Point; + using Bezier = primal::BezierCurve; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + // Define anonymous functions for testing + auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; + auto poly_integrand = [](Point2D x) -> double { return x[0] * x[1] * x[1]; }; + auto transc_integrand = [](Point2D x) -> double { + return std::sin(x[0] * x[1]); + }; + + // Test on single parabolic segment + Point2D paranodes[] = {Point2D {-1.0, 1.0}, + Point2D {0.5, -2.0}, + Point2D {2.0, 4.0}}; + Bezier parabola_segment(paranodes, 2); + + // Compare against hand computed/high-precision calculated values. + + // Constant integrand line integral is equivalent to arc-length calculation + EXPECT_NEAR(evaluate_line_integral(parabola_segment, const_integrand, npts), + 6.12572661998, + abs_tol); + + EXPECT_NEAR(evaluate_line_integral(parabola_segment, poly_integrand, npts), + 37.8010703669, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(parabola_segment, transc_integrand, npts), + 0.495907795678, + abs_tol); + + // Test on a collection of Bezier curves + Point2D segnodes1[] = {Point2D {-1.0, -1.0}, + Point2D {-1.0 / 3.0, 1.0}, + Point2D {1.0 / 3.0, -1.0}, + Point2D {1.0, 1.0}}; + Bezier cubic_segment(segnodes1, 3); + + Point2D segnodes2[] = {Point2D {1.0, 1.0}, Point2D {-1.0, 0.0}}; + Bezier linear_segment(segnodes2, 1); + + Point2D segnodes3[] = {Point2D {-1.0, 0.0}, + Point2D {-3.0, 1.0}, + Point2D {-1.0, 2.0}}; + Bezier quadratic_segment(segnodes3, 2); + + axom::Array connected_curve( + {cubic_segment, linear_segment, quadratic_segment}); + EXPECT_NEAR(evaluate_line_integral(connected_curve, const_integrand, npts), + 8.28968500196, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(connected_curve, poly_integrand, npts), + -5.97565740064, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(connected_curve, transc_integrand, npts), + -0.574992518405, + abs_tol); + + // Test algorithm on disconnected curves + axom::Array disconnected_curve({cubic_segment, quadratic_segment}); + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, const_integrand, npts), + 6.05361702446, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, poly_integrand, npts), + -6.34833539689, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(disconnected_curve, transc_integrand, npts), + -0.914161242161, + abs_tol); +} + +TEST(primal_integral, evaluate_line_integral_vector) +{ + using Point2D = primal::Point; + using Vector2D = primal::Vector; + using Bezier = primal::BezierCurve; + double abs_tol = 1e-10; + + // Quadrature nodes. Should be sufficiently high to pass tests + int npts = 30; + + // Test on a single line segment + auto vec_field = [](Point2D x) -> Vector2D { + return Vector2D({x[1] * x[1], 3 * x[0] - 6 * x[1]}); + }; + + Point2D segnodes[] = {Point2D {3.0, 7.0}, Point2D {0.0, 12.0}}; + Bezier linear_segment(segnodes, 1); + + // Compare against hand computed values + EXPECT_NEAR(evaluate_line_integral(linear_segment, vec_field, npts), + -1079.0 / 2.0, + abs_tol); + + // Test on a closed curve + auto area_field = [](Point2D x) -> Vector2D { + return Vector2D({-0.5 * x[1], 0.5 * x[0]}); + }; + auto conservative_field = [](Point2D x) -> Vector2D { + return Vector2D({2 * x[0] * x[1] * x[1], 2 * x[0] * x[0] * x[1]}); + }; + auto winding_field = [](Point2D x) -> Vector2D { + double denom = 2 * M_PI * (x[0] * x[0] + x[1] * x[1]); + return Vector2D({-x[1] / denom, x[0] / denom}); + }; + + Point2D paranodes1[] = {Point2D {1.0, 0.0}, + Point2D {0.0, 2.0}, + Point2D {-1.0, 0.0}}; + Bezier para1(paranodes1, 2); + + Point2D paranodes2[] = {Point2D {-1.0, 0.0}, + Point2D {0.0, -2.0}, + Point2D {1.0, 0.0}}; + Bezier para2(paranodes2, 2); + + axom::Array parabola_shape({para1, para2}); + + // This vector field calculates the area of the region + EXPECT_NEAR(evaluate_line_integral(parabola_shape, area_field, npts), + 8.0 / 3.0, + abs_tol); + + // This vector field is conservative, so it should evaluate to zero + EXPECT_NEAR(evaluate_line_integral(parabola_shape, conservative_field, npts), + 0.0, + abs_tol); + + // This vector field is generated by a in/out query, should return 1 (inside) + EXPECT_NEAR(evaluate_line_integral(parabola_shape, winding_field, npts), + 1.0, + abs_tol); + + // Test algorithm on disconnected curves + Point2D paranodes2_shifted[] = {Point2D {-1.0, -1.0}, + Point2D {0.0, -3.0}, + Point2D {1.0, -1.0}}; + Bezier para2_shift(paranodes2_shifted, 2); + axom::Array disconnected_parabola_shape({para1, para2_shift}); + + EXPECT_NEAR( + evaluate_line_integral(disconnected_parabola_shape, area_field, npts), + 11.0 / 3.0, + abs_tol); + EXPECT_NEAR( + evaluate_line_integral(disconnected_parabola_shape, conservative_field, npts), + 0.0, + abs_tol); + EXPECT_NEAR( + evaluate_line_integral(disconnected_parabola_shape, winding_field, npts), + 0.75, + abs_tol); +} + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + + axom::slic::SimpleLogger logger; + + int result = RUN_ALL_TESTS(); + + return result; +} diff --git a/src/docs/sphinx/quickstart_guide/config_build.rst b/src/docs/sphinx/quickstart_guide/config_build.rst index 2a236a0a34..2a4cd601c9 100644 --- a/src/docs/sphinx/quickstart_guide/config_build.rst +++ b/src/docs/sphinx/quickstart_guide/config_build.rst @@ -76,7 +76,7 @@ The following table lists: `c2c`_ Optional: Quest C2C_DIR `HDF5`_ Optional: Sidre HDF5_DIR `Lua`_ Optional: Inlet LUA_DIR - `MFEM`_ Optional: Quest, Sidre MFEM_DIR + `MFEM`_ Optional: Primal, Quest, Sidre MFEM_DIR `RAJA`_ Optional: Mint, Spin, Quest RAJA_DIR `SCR`_ Optional: Sidre SCR_DIR `Umpire`_ Optional: Core, Spin, Quest UMPIRE_DIR