From 4ecd1b98116ac98a7ac0f3c78602a8f56c3a86a4 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 11:57:58 -0700 Subject: [PATCH 01/22] Added evaluate_integral.hpp, changed CMake file, add optional dependancy of primal on mfem --- src/axom/primal/CMakeLists.txt | 2 + .../primal/operators/evaluate_integral.hpp | 339 ++++++++++++++++++ 2 files changed, 341 insertions(+) create mode 100644 src/axom/primal/operators/evaluate_integral.hpp diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 920fd4a3b8..8aa337c77a 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -45,6 +45,7 @@ set( primal_headers operators/compute_bounding_box.hpp operators/compute_moments.hpp operators/in_sphere.hpp + operators/evaluate_integral.hpp operators/split.hpp operators/detail/clip_impl.hpp @@ -78,6 +79,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/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp new file mode 100644 index 0000000000..3556d3e3ab --- /dev/null +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -0,0 +1,339 @@ +// 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 simple 1D Gaussian Quadrature using rules generated + * by MFEM. 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar + * Regions Bounded by Rational Parametric Curves" by David Gunerman et al. + */ + +#ifndef QUEST_EVAL_INTEGRAL_INTERFACE_HPP_ +#define QUEST_EVAL_INTEGRAL_INTERFACE_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 + +// Create type definitions the different field types (vector vs scalar) +// to allow for proper overloading of line integrals +using Point2D = axom::primal::Point; +using Vector2D = axom::primal::Vector; + +typedef std::function vector_field; +typedef std::function scalar_field; + +namespace axom +{ +namespace primal +{ + +/*! + * \brief Computes the tangent vector to a Bezier curve. + * + * Uses the de Casteljau algorithm. Functionality will eventually + * be supplanted by Bezier curve class method. + * + * \param [in] c the Bezier curve + * \param [in] t on [0, 1] the parameter value on which to evaluate the derivative + * \return the tangent vector + */ +template +primal::Vector temp_dt(const primal::BezierCurve& c, T t) +{ + Vector val; + + const int ord = c.getOrder(); + std::vector dCarray(ord + 1); + + // Run de Casteljau algorithm on each dimension + for(int i = 0; i < 2; ++i) + { + for(int p = 0; p <= ord; ++p) + { + dCarray[p] = c[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 Evaluate a line integral along a collection of Bezier curves + * + * The curve need not be connected, simply adding the portion along each segment. + * Uses 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 + 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 < cs.size(); i++) + { + // Compute the line integral along each component + total_integral += evaluate_line_integral(cs[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 + mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule* quad = + &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); + + return evaluate_line_integral(c, integrand, quad); +} + +/*! + * \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(const primal::BezierCurve& c, + scalar_field integrand, + const mfem::IntegrationRule* quad) +{ + // Store intermediate values + Point2D x_q; + Vector2D dx_q; + + // 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 + x_q = c.evaluate(quad->IntPoint(q).x); + dx_q = temp_dt(c, 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(const primal::BezierCurve& c, + vector_field vec_field, + const mfem::IntegrationRule* quad) +{ + // Store intermediate values. + Point2D x_q; + Vector2D dx_q, func_val; + + // 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 + x_q = c.evaluate(quad->IntPoint(q).x); + dx_q = temp_dt(c, quad->IntPoint(q).x); + func_val = vec_field(x_q); + + full_quadrature += + quad->IntPoint(q).weight * Vector2D::dot_product(func_val, dx_q); + } + + return full_quadrature; +} + +/*! + * \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. + * Will eventually incorperate the CurverdPolygon class. + * 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 + mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + const mfem::IntegrationRule *quad_Q, *quad_P; + + // Get the quadrature for the line integral. + // Quadrature order is equal to 2*N - 1 + quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); + + // If the quadrature orders are the same (which they usually are), + // then only create one quadrature within the library and then + // the two pointers can point to the same object. Hopefully improves efficiency. + if(npts_P <= 0) + { + npts_P = npts_Q; + quad_P = quad_Q; + } + else + { + quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1)); + } + + // Define lower bound as lowest control node for all curves + double int_lb = cs[0][0][1]; // Start as y-component of first node + 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(int i = 0; i < cs.size(); i++) + { + total_integral += + evaluate_area_integral(cs[i], integrand, int_lb, quad_Q, quad_P); + } + + return total_integral; +} + +/*! + * \brief Evaluate the integral across a single component of the 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. + * + * \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(const primal::BezierCurve& c, + Lambda&& integrand, + double int_lb, + const mfem::IntegrationRule* quad_Q, + const mfem::IntegrationRule* quad_P) +{ + // Store some intermediate values + Point2D x_q, x_qxi; + 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 + 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 + x_qxi[0] = x_q[0]; + x_qxi[1] = (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); + + // TODO: Check if the antiderivative should be negated for Green's theorem + full_quadrature += quad_Q->IntPoint(q).weight * + temp_dt(c, quad_Q->IntPoint(q).x)[0] * -antiderivative; + //full_quadrature += quad_Q->IntPoint(q).weight * + // temp_dt(c, quad_Q->IntPoint(q).x)[0] * antiderivative; + } + } + + return full_quadrature; +} + +} // namespace primal +} // end namespace axom + +#endif From 3aa6a0414cb798b89947f9fecfb200cd86d2b8cd Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 11:58:56 -0700 Subject: [PATCH 02/22] Added test for evaluate_integral methods, changed CMake file --- src/axom/primal/tests/CMakeLists.txt | 3 +- src/axom/primal/tests/primal_integral.cpp | 212 ++++++++++++++++++++++ 2 files changed, 214 insertions(+), 1 deletion(-) create mode 100644 src/axom/primal/tests/primal_integral.cpp diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index 6d246266bb..c47a45061b 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -16,6 +16,7 @@ set( primal_tests primal_compute_moments.cpp primal_curved_polygon.cpp primal_in_sphere.cpp + primal_integral.cpp primal_intersect.cpp primal_intersect_impl.cpp primal_numeric_array.cpp @@ -61,4 +62,4 @@ foreach ( test ${primal_tests} ) COMMAND ${test_name}_test ) -endforeach() +endforeach() \ No newline at end of file diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp new file mode 100644 index 0000000000..6ac965b2af --- /dev/null +++ b/src/axom/primal/tests/primal_integral.cpp @@ -0,0 +1,212 @@ +// primal_integral_test.cpp : This file contains the 'main' function. Program execution begins and ends there. + +#include "axom/primal.hpp" +#include "axom/slic.hpp" +#include "axom/fmt.hpp" + +#include "gtest/gtest.h" + +#include + +namespace primal = axom::primal; +namespace slic = axom::slic; + + +/*! + * \brief Utility function to initialize the logger + */ +void initializeLogger() +{ + // Initialize Logger + slic::initialize(); + slic::setLoggingMsgLevel(axom::slic::message::Info); + + slic::LogStream* logStream; + +#ifdef AXOM_USE_MPI + std::string fmt = "[][]: \n"; + #ifdef AXOM_USE_LUMBERJACK + const int RLIMIT = 8; + logStream = new slic::LumberjackStream(&std::cout, MPI_COMM_WORLD, RLIMIT, fmt); + #else + logStream = new slic::SynchronizedStream(&std::cout, MPI_COMM_WORLD, fmt); + #endif +#else + std::string fmt = "[]: \n"; + logStream = new slic::GenericOutputStream(&std::cout, fmt); +#endif // AXOM_USE_MPI + + slic::addStreamToAllMsgLevels(logStream); +} + +/*! + * \brief Utility function to finalize the logger + */ +void finalizeLogger() +{ + if(slic::isInitialized()) + { + slic::flushStreams(); + slic::finalize(); + } +} + +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 = 15; + + // 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 + double trinodes1[] = {0.0, 1.0, 0.0, 0.0}; + Bezier tri1(trinodes1, 1); + + double trinodes2[] = {1.0, 0.0, 0.0, 1.0}; + Bezier tri2(trinodes2, 1); + + double trinodes3[] = {0.0, 0.0, 1.0, 0.0}; + Bezier tri3(trinodes3, 1); + + axom::Array triangle({tri1, tri2, tri3}); + 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) + double paranodes1[] = {2.0, 1.0, 0.0, + 0.0, 2.0, 0.0}; + Bezier para1(paranodes1, 2); + + double paranodes2[] = {0.0, 1.0, 2.0, + 0.0, -2.0, 0.0}; + Bezier para2(paranodes2, 2); + + axom::Array parabola_shape({para1, para2}); + 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); +} + +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 + double paranodes[] = {-1.0, 0.5, 2.0, + 1.0, -2.0, 4.0}; + Bezier parabola_segment(paranodes, 2); + + 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 + double segnodes1[] = {-1.0, -1.0/3.0, 1.0/3.0, 1.0, + -1.0, 1.0, -1.0, 1.0}; + Bezier cubic_segment(segnodes1, 3); + + double segnodes2[] = {1.0, -1.0, + 1.0, 0.0}; + Bezier linear_segment(segnodes2, 1); + + double segnodes3[] = {-1.0, -3.0, -1.0, + 0.0, 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(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]}); + }; + + double segnodes[] = {3.0, 0.0, 7.0, 12.0}; + Bezier linear_segment(segnodes, 1); + + 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}); + }; + + double paranodes1[] = {1.0, 0.0, -1.0, + 0.0, 2.0, 0.0}; + Bezier para1(paranodes1, 2); + + double paranodes2[] = {-1.0, 0.0, 1.0, + 0.0, -2.0, 0.0}; + Bezier para2(paranodes2, 2); + + axom::Array parabola_shape({para1, para2}); + EXPECT_NEAR(evaluate_line_integral(parabola_shape, area_field, npts), + 8.0 / 3.0, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(parabola_shape, conservative_field, npts), + 0.0, + abs_tol); + EXPECT_NEAR(evaluate_line_integral(parabola_shape, winding_field, npts), + 1.0, + abs_tol); +} + + +int main(int argc, char* argv[]) +{ + // -- Initialize logger + initializeLogger(); + + ::testing::InitGoogleTest(&argc, argv); + + int result = RUN_ALL_TESTS(); + + // -- Finalize logger + finalizeLogger(); + + return 0; +} \ No newline at end of file From 52a24a7a84084591a9b72a6c7d079a80f96f04a3 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 12:26:11 -0700 Subject: [PATCH 03/22] Added copyright info --- src/axom/primal/tests/primal_integral.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 6ac965b2af..de388fcb53 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -1,4 +1,7 @@ -// primal_integral_test.cpp : This file contains the 'main' function. Program execution begins and ends there. +// 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" From ff23aadd292b34078de6fe0f7e9b36c770054baf Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 16:12:31 -0700 Subject: [PATCH 04/22] Corrected style --- src/axom/primal/tests/primal_integral.cpp | 87 +++++++++++++---------- 1 file changed, 51 insertions(+), 36 deletions(-) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index de388fcb53..da6c1ddfe5 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -14,7 +14,6 @@ namespace primal = axom::primal; namespace slic = axom::slic; - /*! * \brief Utility function to initialize the logger */ @@ -79,25 +78,35 @@ TEST(primal_integral, evaluate_area_integral) double trinodes3[] = {0.0, 0.0, 1.0, 0.0}; Bezier tri3(trinodes3, 1); - + axom::Array triangle({tri1, tri2, tri3}); - 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); + 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) - double paranodes1[] = {2.0, 1.0, 0.0, - 0.0, 2.0, 0.0}; + double paranodes1[] = {2.0, 1.0, 0.0, 0.0, 2.0, 0.0}; Bezier para1(paranodes1, 2); - double paranodes2[] = {0.0, 1.0, 2.0, - 0.0, -2.0, 0.0}; + double paranodes2[] = {0.0, 1.0, 2.0, 0.0, -2.0, 0.0}; Bezier para2(paranodes2, 2); - axom::Array parabola_shape({para1, para2}); - 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); + axom::Array parabola_shape({para1, para2}); + 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); } TEST(primal_integral, evaluate_line_integral_scalar) @@ -117,31 +126,40 @@ TEST(primal_integral, evaluate_line_integral_scalar) }; // Test on single parabolic segment - double paranodes[] = {-1.0, 0.5, 2.0, - 1.0, -2.0, 4.0}; + double paranodes[] = {-1.0, 0.5, 2.0, 1.0, -2.0, 4.0}; Bezier parabola_segment(paranodes, 2); - 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); + 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 - double segnodes1[] = {-1.0, -1.0/3.0, 1.0/3.0, 1.0, - -1.0, 1.0, -1.0, 1.0}; + double segnodes1[] = {-1.0, -1.0 / 3.0, 1.0 / 3.0, 1.0, -1.0, 1.0, -1.0, 1.0}; Bezier cubic_segment(segnodes1, 3); - double segnodes2[] = {1.0, -1.0, - 1.0, 0.0}; + double segnodes2[] = {1.0, -1.0, 1.0, 0.0}; Bezier linear_segment(segnodes2, 1); - double segnodes3[] = {-1.0, -3.0, -1.0, - 0.0, 1.0, 2.0}; + double segnodes3[] = {-1.0, -3.0, -1.0, 0.0, 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); + 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(primal_integral, evaluate_line_integral_vector) @@ -168,7 +186,7 @@ TEST(primal_integral, evaluate_line_integral_vector) // Test on a closed curve auto area_field = [](Point2D x) -> Vector2D { - return Vector2D({-0.5*x[1], 0.5*x[0]}); + 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]}); @@ -177,13 +195,11 @@ TEST(primal_integral, evaluate_line_integral_vector) double denom = 2 * M_PI * (x[0] * x[0] + x[1] * x[1]); return Vector2D({-x[1] / denom, x[0] / denom}); }; - - double paranodes1[] = {1.0, 0.0, -1.0, - 0.0, 2.0, 0.0}; + + double paranodes1[] = {1.0, 0.0, -1.0, 0.0, 2.0, 0.0}; Bezier para1(paranodes1, 2); - double paranodes2[] = {-1.0, 0.0, 1.0, - 0.0, -2.0, 0.0}; + double paranodes2[] = {-1.0, 0.0, 1.0, 0.0, -2.0, 0.0}; Bezier para2(paranodes2, 2); axom::Array parabola_shape({para1, para2}); @@ -198,12 +214,11 @@ TEST(primal_integral, evaluate_line_integral_vector) abs_tol); } - int main(int argc, char* argv[]) { // -- Initialize logger initializeLogger(); - + ::testing::InitGoogleTest(&argc, argv); int result = RUN_ALL_TESTS(); From 2095c3e536c28f9cd9901fd4ab3eaedfb3104a7c Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 16:44:54 -0700 Subject: [PATCH 05/22] Fixed style, slic logging --- .../primal/operators/evaluate_integral.hpp | 1 - src/axom/primal/tests/primal_integral.cpp | 52 ++----------------- 2 files changed, 4 insertions(+), 49 deletions(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 3556d3e3ab..df59df60ef 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -43,7 +43,6 @@ namespace axom { namespace primal { - /*! * \brief Computes the tangent vector to a Bezier curve. * diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index da6c1ddfe5..ec77e16267 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -12,46 +12,6 @@ #include namespace primal = axom::primal; -namespace slic = axom::slic; - -/*! - * \brief Utility function to initialize the logger - */ -void initializeLogger() -{ - // Initialize Logger - slic::initialize(); - slic::setLoggingMsgLevel(axom::slic::message::Info); - - slic::LogStream* logStream; - -#ifdef AXOM_USE_MPI - std::string fmt = "[][]: \n"; - #ifdef AXOM_USE_LUMBERJACK - const int RLIMIT = 8; - logStream = new slic::LumberjackStream(&std::cout, MPI_COMM_WORLD, RLIMIT, fmt); - #else - logStream = new slic::SynchronizedStream(&std::cout, MPI_COMM_WORLD, fmt); - #endif -#else - std::string fmt = "[]: \n"; - logStream = new slic::GenericOutputStream(&std::cout, fmt); -#endif // AXOM_USE_MPI - - slic::addStreamToAllMsgLevels(logStream); -} - -/*! - * \brief Utility function to finalize the logger - */ -void finalizeLogger() -{ - if(slic::isInitialized()) - { - slic::flushStreams(); - slic::finalize(); - } -} TEST(primal_integral, evaluate_area_integral) { @@ -216,15 +176,11 @@ TEST(primal_integral, evaluate_line_integral_vector) int main(int argc, char* argv[]) { - // -- Initialize logger - initializeLogger(); - ::testing::InitGoogleTest(&argc, argv); - int result = RUN_ALL_TESTS(); + axom::slic::SimpleLogger logger; - // -- Finalize logger - finalizeLogger(); + int result = RUN_ALL_TESTS(); - return 0; -} \ No newline at end of file + return result; +} From b5496b144ed8b78120b68d20e22954a0aa2df321 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 16:49:07 -0700 Subject: [PATCH 06/22] Fixed lambda function warning --- src/axom/primal/tests/primal_integral.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index ec77e16267..0315350705 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -23,7 +23,7 @@ TEST(primal_integral, evaluate_area_integral) int npts = 15; // Define anonymous functions for testing - auto const_integrand = [](Point2D x) -> double { return 1.0; }; + 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]); @@ -79,7 +79,7 @@ TEST(primal_integral, evaluate_line_integral_scalar) int npts = 30; // Define anonymous functions for testing - auto const_integrand = [](Point2D x) -> double { return 1.0; }; + 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]); From 67212d7094452506514158d28658816ebcc5b8e6 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Tue, 21 Jun 2022 16:57:21 -0700 Subject: [PATCH 07/22] Fixed dependency of primal::evaluate_integral on mfem in build system --- src/axom/primal/CMakeLists.txt | 3 ++- src/axom/primal/operators/evaluate_integral.hpp | 2 +- src/axom/primal/tests/CMakeLists.txt | 3 ++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 8aa337c77a..c4c83120e6 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -45,7 +45,6 @@ set( primal_headers operators/compute_bounding_box.hpp operators/compute_moments.hpp operators/in_sphere.hpp - operators/evaluate_integral.hpp operators/split.hpp operators/detail/clip_impl.hpp @@ -69,6 +68,8 @@ set( primal_sources operators/detail/intersect_impl.cpp ) +blt_list_append( TO primal_headers ELEMENTS operators/evaluate_integral.hpp IF MFEM_FOUND ) + #------------------------------------------------------------------------------ # Build and install the library #------------------------------------------------------------------------------ diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index df59df60ef..fa382af9db 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -11,7 +11,7 @@ * * Line integrals are computed with simple 1D Gaussian Quadrature using rules generated * by MFEM. 2D area integrals computed with "Spectral Mesh-Free Quadrature for Planar - * Regions Bounded by Rational Parametric Curves" by David Gunerman et al. + * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. */ #ifndef QUEST_EVAL_INTEGRAL_INTERFACE_HPP_ diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index c47a45061b..aeeb4f11a8 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -16,7 +16,6 @@ set( primal_tests primal_compute_moments.cpp primal_curved_polygon.cpp primal_in_sphere.cpp - primal_integral.cpp primal_intersect.cpp primal_intersect_impl.cpp primal_numeric_array.cpp @@ -39,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 ) From 5a430556f3fd2e2a4d4a3e74df5fa1ab3dd3e015 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 12:09:41 -0700 Subject: [PATCH 08/22] Add implementation file to minimize interface, change CMake --- src/axom/primal/CMakeLists.txt | 7 +- .../detail/evaluate_integral_impl.hpp | 197 ++++++++++++++++++ 2 files changed, 203 insertions(+), 1 deletion(-) create mode 100644 src/axom/primal/operators/detail/evaluate_integral_impl.hpp diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index c4c83120e6..5b493c6686 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -68,7 +68,12 @@ set( primal_sources operators/detail/intersect_impl.cpp ) -blt_list_append( TO primal_headers ELEMENTS operators/evaluate_integral.hpp IF MFEM_FOUND ) +blt_list_append( TO primal_headers IF MFEM_FOUND ELEMENTS + + ## operators + operators/evaluate_integral.hpp + operators/detail/evaluate_integral_impl.hpp + ) #------------------------------------------------------------------------------ # Build and install the library 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..099ae87e88 --- /dev/null +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -0,0 +1,197 @@ +// 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 Computes the tangent vector to a Bezier curve. + * + * Uses the de Casteljau algorithm. Functionality will eventually + * be supplanted by Bezier curve class method. + * + * \param [in] c the Bezier curve + * \param [in] t on [0, 1] the parameter value on which to evaluate the derivative + * \return the tangent vector + */ +template +primal::Vector temp_dt(const primal::BezierCurve& c, T t) +{ + Vector val; + + const int ord = c.getOrder(); + std::vector dCarray(ord + 1); + + // Run de Casteljau algorithm on each dimension + for(int i = 0; i < 2; ++i) + { + for(int p = 0; p <= ord; ++p) + { + dCarray[p] = c[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 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 intermediate values + Point2D x_q; + Vector2D dx_q; + + // 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 + x_q = c.evaluate(quad->IntPoint(q).x); + dx_q = detail::temp_dt(c, 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 intermediate values. + Point2D x_q; + Vector2D dx_q, func_val; + + // 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 + x_q = c.evaluate(quad->IntPoint(q).x); + dx_q = detail::temp_dt(c, quad->IntPoint(q).x); + 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 integral across a single component of the 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. + * + * \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 + Point2D x_q, x_qxi; + 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 + 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 + x_qxi[0] = x_q[0]; + x_qxi[1] = (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 * + detail::temp_dt(c, 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 From 64c9dc371aa96b6d89c5f407baa777a74ae977da Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 12:10:55 -0700 Subject: [PATCH 09/22] Minimize interface, add static quadrature rules library --- .../primal/operators/evaluate_integral.hpp | 192 ++---------------- 1 file changed, 15 insertions(+), 177 deletions(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index fa382af9db..46b41c00be 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -14,13 +14,15 @@ * Regions Bounded by Rational Parametric Curves" by David Gunderman et al. */ -#ifndef QUEST_EVAL_INTEGRAL_INTERFACE_HPP_ -#define QUEST_EVAL_INTEGRAL_INTERFACE_HPP_ +#ifndef PRIMAL_EVAL_INTEGRAL_HPP_ +#define PRIMAL_EVAL_INTEGRAL_HPP_ // Axom includes #include "axom/config.hpp" // for compile-time configuration options #include "axom/primal.hpp" +#include "axom/primal/operators/detail/evaluate_integral_impl.hpp" + // MFEM includes #ifdef AXOM_USE_MFEM #include "mfem.hpp" @@ -36,54 +38,10 @@ using Point2D = axom::primal::Point; using Vector2D = axom::primal::Vector; -typedef std::function vector_field; -typedef std::function scalar_field; - namespace axom { namespace primal { -/*! - * \brief Computes the tangent vector to a Bezier curve. - * - * Uses the de Casteljau algorithm. Functionality will eventually - * be supplanted by Bezier curve class method. - * - * \param [in] c the Bezier curve - * \param [in] t on [0, 1] the parameter value on which to evaluate the derivative - * \return the tangent vector - */ -template -primal::Vector temp_dt(const primal::BezierCurve& c, T t) -{ - Vector val; - - const int ord = c.getOrder(); - std::vector dCarray(ord + 1); - - // Run de Casteljau algorithm on each dimension - for(int i = 0; i < 2; ++i) - { - for(int p = 0; p <= ord; ++p) - { - dCarray[p] = c[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 Evaluate a line integral along a collection of Bezier curves * @@ -105,7 +63,7 @@ double evaluate_line_integral(const axom::Array>& // 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 - mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); const mfem::IntegrationRule* quad = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); @@ -113,7 +71,8 @@ double evaluate_line_integral(const axom::Array>& for(int i = 0; i < cs.size(); i++) { // Compute the line integral along each component - total_integral += evaluate_line_integral(cs[i], integrand, quad); + total_integral += + detail::evaluate_line_integral_component(cs[i], integrand, quad); } return total_integral; @@ -140,81 +99,11 @@ double evaluate_line_integral(const primal::BezierCurve& c, // 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 - mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); const mfem::IntegrationRule* quad = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); - return evaluate_line_integral(c, integrand, quad); -} - -/*! - * \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(const primal::BezierCurve& c, - scalar_field integrand, - const mfem::IntegrationRule* quad) -{ - // Store intermediate values - Point2D x_q; - Vector2D dx_q; - - // 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 - x_q = c.evaluate(quad->IntPoint(q).x); - dx_q = temp_dt(c, 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(const primal::BezierCurve& c, - vector_field vec_field, - const mfem::IntegrationRule* quad) -{ - // Store intermediate values. - Point2D x_q; - Vector2D dx_q, func_val; - - // 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 - x_q = c.evaluate(quad->IntPoint(q).x); - dx_q = temp_dt(c, quad->IntPoint(q).x); - func_val = vec_field(x_q); - - full_quadrature += - quad->IntPoint(q).weight * Vector2D::dot_product(func_val, dx_q); - } - - return full_quadrature; + return detail::evaluate_line_integral_component(c, integrand, quad); } /*! @@ -241,7 +130,7 @@ double evaluate_area_integral(const axom::Array>& { // Generate quadrature library, defaulting to GaussLegendre quadrature. // Use the same one for every curve in the polygon - mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); + static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); const mfem::IntegrationRule *quad_Q, *quad_P; // Get the quadrature for the line integral. @@ -271,67 +160,16 @@ double evaluate_area_integral(const axom::Array>& double total_integral = 0.0; for(int i = 0; i < cs.size(); i++) { - total_integral += - evaluate_area_integral(cs[i], integrand, int_lb, quad_Q, quad_P); + total_integral += detail::evaluate_area_integral_component(cs[i], + integrand, + int_lb, + quad_Q, + quad_P); } return total_integral; } -/*! - * \brief Evaluate the integral across a single component of the 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. - * - * \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(const primal::BezierCurve& c, - Lambda&& integrand, - double int_lb, - const mfem::IntegrationRule* quad_Q, - const mfem::IntegrationRule* quad_P) -{ - // Store some intermediate values - Point2D x_q, x_qxi; - 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 - 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 - x_qxi[0] = x_q[0]; - x_qxi[1] = (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); - - // TODO: Check if the antiderivative should be negated for Green's theorem - full_quadrature += quad_Q->IntPoint(q).weight * - temp_dt(c, quad_Q->IntPoint(q).x)[0] * -antiderivative; - //full_quadrature += quad_Q->IntPoint(q).weight * - // temp_dt(c, quad_Q->IntPoint(q).x)[0] * antiderivative; - } - } - - return full_quadrature; -} - } // namespace primal } // end namespace axom From b9d0005ba8e42e30315028745477bda3fb1bd6ea Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 12:11:13 -0700 Subject: [PATCH 10/22] Add comments, increase readability --- src/axom/primal/tests/primal_integral.cpp | 55 +++++++++++++++++------ 1 file changed, 42 insertions(+), 13 deletions(-) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 0315350705..cbe1ebb9d3 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -20,7 +20,7 @@ TEST(primal_integral, evaluate_area_integral) double abs_tol = 1e-10; // Quadrature nodes. Should be sufficiently high to pass tests - int npts = 15; + int npts = 20; // Define anonymous functions for testing auto const_integrand = [](Point2D /*x*/) -> double { return 1.0; }; @@ -30,16 +30,18 @@ TEST(primal_integral, evaluate_area_integral) }; // Test on triangular domain - double trinodes1[] = {0.0, 1.0, 0.0, 0.0}; + Point2D trinodes1[] = {Point2D {0.0, 0.0}, Point2D {1.0, 0.0}}; Bezier tri1(trinodes1, 1); - double trinodes2[] = {1.0, 0.0, 0.0, 1.0}; + Point2D trinodes2[] = {Point2D {1.0, 0.0}, Point2D {0.0, 1.0}}; Bezier tri2(trinodes2, 1); - double trinodes3[] = {0.0, 0.0, 1.0, 0.0}; + 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); @@ -51,13 +53,19 @@ TEST(primal_integral, evaluate_area_integral) abs_tol); // Test on parabolic domain (between f(x) = 1-x^2 and g(x) = x^2-1, shifted to the right 1 unit) - double paranodes1[] = {2.0, 1.0, 0.0, 0.0, 2.0, 0.0}; + Point2D paranodes1[] = {Point2D {2.0, 0.0}, + Point2D {1.0, 2.0}, + Point2D {0.0, 0.0}}; Bezier para1(paranodes1, 2); - double paranodes2[] = {0.0, 1.0, 2.0, 0.0, -2.0, 0.0}; + 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); @@ -86,9 +94,12 @@ TEST(primal_integral, evaluate_line_integral_scalar) }; // Test on single parabolic segment - double paranodes[] = {-1.0, 0.5, 2.0, 1.0, -2.0, 4.0}; + 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 EXPECT_NEAR(evaluate_line_integral(parabola_segment, const_integrand, npts), 6.12572661998, abs_tol); @@ -100,17 +111,24 @@ TEST(primal_integral, evaluate_line_integral_scalar) abs_tol); // Test on a collection of Bezier curves - double segnodes1[] = {-1.0, -1.0 / 3.0, 1.0 / 3.0, 1.0, -1.0, 1.0, -1.0, 1.0}; + 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); - double segnodes2[] = {1.0, -1.0, 1.0, 0.0}; + Point2D segnodes2[] = {Point2D {1.0, 1.0}, Point2D {-1.0, 0.0}}; Bezier linear_segment(segnodes2, 1); - double segnodes3[] = {-1.0, -3.0, -1.0, 0.0, 1.0, 2.0}; + 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}); + + // Compare against hand computed/high-precision calculated values EXPECT_NEAR(evaluate_line_integral(connected_curve, const_integrand, npts), 8.28968500196, abs_tol); @@ -137,9 +155,10 @@ TEST(primal_integral, evaluate_line_integral_vector) return Vector2D({x[1] * x[1], 3 * x[0] - 6 * x[1]}); }; - double segnodes[] = {3.0, 0.0, 7.0, 12.0}; + 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); @@ -156,19 +175,29 @@ TEST(primal_integral, evaluate_line_integral_vector) return Vector2D({-x[1] / denom, x[0] / denom}); }; - double paranodes1[] = {1.0, 0.0, -1.0, 0.0, 2.0, 0.0}; + Point2D paranodes1[] = {Point2D {1.0, 0.0}, + Point2D {0.0, 2.0}, + Point2D {-1.0, 0.0}}; Bezier para1(paranodes1, 2); - double paranodes2[] = {-1.0, 0.0, 1.0, 0.0, -2.0, 0.0}; + 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); From 3b5715158d016eecf0f68ead2841d921b0f6cf36 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 14:22:35 -0700 Subject: [PATCH 11/22] Update primal/CMakeLists.txt formatting Co-authored-by: Kenny Weiss --- src/axom/primal/CMakeLists.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index 5b493c6686..cf607de8ed 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -68,11 +68,11 @@ set( primal_sources operators/detail/intersect_impl.cpp ) -blt_list_append( TO primal_headers IF MFEM_FOUND ELEMENTS - - ## operators - operators/evaluate_integral.hpp - operators/detail/evaluate_integral_impl.hpp +blt_list_append( + TO primal_headers + ELEMENTS operators/evaluate_integral.hpp + operators/detail/evaluate_integral_impl.hpp + IF MFEM_FOUND ) #------------------------------------------------------------------------------ From 1d61386c9e91b9be767ea3090a3fb0259e23706a Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 16:17:51 -0700 Subject: [PATCH 12/22] Added tests for CurvedPolygon interface. --- src/axom/primal/tests/CMakeLists.txt | 2 +- src/axom/primal/tests/primal_integral.cpp | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/tests/CMakeLists.txt b/src/axom/primal/tests/CMakeLists.txt index aeeb4f11a8..6b11631efd 100644 --- a/src/axom/primal/tests/CMakeLists.txt +++ b/src/axom/primal/tests/CMakeLists.txt @@ -63,4 +63,4 @@ foreach ( test ${primal_tests} ) COMMAND ${test_name}_test ) -endforeach() \ No newline at end of file +endforeach() diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index cbe1ebb9d3..65ba9b19be 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -75,6 +75,19 @@ TEST(primal_integral, evaluate_area_integral) 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) @@ -210,6 +223,6 @@ int main(int argc, char* argv[]) axom::slic::SimpleLogger logger; int result = RUN_ALL_TESTS(); - + return result; } From 4995da13e6b7186b367b45343758301caaa3a1d3 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 16:20:54 -0700 Subject: [PATCH 13/22] Updated style/formatting, compilation optimization --- .../detail/evaluate_integral_impl.hpp | 75 +++------------ .../primal/operators/evaluate_integral.hpp | 94 +++++++++++++------ src/axom/primal/tests/primal_integral.cpp | 2 +- 3 files changed, 80 insertions(+), 91 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 099ae87e88..2aa845d09c 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -26,47 +26,6 @@ namespace primal { namespace detail { -/*! - * \brief Computes the tangent vector to a Bezier curve. - * - * Uses the de Casteljau algorithm. Functionality will eventually - * be supplanted by Bezier curve class method. - * - * \param [in] c the Bezier curve - * \param [in] t on [0, 1] the parameter value on which to evaluate the derivative - * \return the tangent vector - */ -template -primal::Vector temp_dt(const primal::BezierCurve& c, T t) -{ - Vector val; - - const int ord = c.getOrder(); - std::vector dCarray(ord + 1); - - // Run de Casteljau algorithm on each dimension - for(int i = 0; i < 2; ++i) - { - for(int p = 0; p <= ord; ++p) - { - dCarray[p] = c[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 Evaluate a scalar field line integral on a single Bezier curve. * @@ -83,18 +42,14 @@ inline double evaluate_line_integral_component( std::function integrand, const mfem::IntegrationRule* quad) { - // Store intermediate values - Point2D x_q; - Vector2D dx_q; - // 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 - x_q = c.evaluate(quad->IntPoint(q).x); - dx_q = detail::temp_dt(c, quad->IntPoint(q).x); + 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(); } @@ -118,19 +73,15 @@ inline double evaluate_line_integral_component( std::function vec_field, const mfem::IntegrationRule* quad) { - // Store intermediate values. - Point2D x_q; - Vector2D dx_q, func_val; - // 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 - x_q = c.evaluate(quad->IntPoint(q).x); - dx_q = detail::temp_dt(c, quad->IntPoint(q).x); - func_val = vec_field(x_q); + 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); @@ -140,11 +91,14 @@ inline double evaluate_line_integral_component( } /*! - * \brief Evaluate the integral across a single component of the curved polygon. + * \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 @@ -161,7 +115,6 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, const mfem::IntegrationRule* quad_P) { // Store some intermediate values - Point2D x_q, x_qxi; double antiderivative = 0.0; // Store/compute quadrature result @@ -170,20 +123,20 @@ double evaluate_area_integral_component(const primal::BezierCurve& c, { // Get intermediate quadrature point // on which to evaluate antiderivative - x_q = c.evaluate(quad_Q->IntPoint(q).x); + 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 - x_qxi[0] = x_q[0]; - x_qxi[1] = (x_q[1] - int_lb) * quad_P->IntPoint(xi).x + int_lb; + 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 * - detail::temp_dt(c, quad_Q->IntPoint(q).x)[0] * -antiderivative; + c.dt(quad_Q->IntPoint(q).x)[0] * -antiderivative; } } diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 46b41c00be..cc0b218bc9 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -9,7 +9,7 @@ * \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 simple 1D Gaussian Quadrature using rules generated + * 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. */ @@ -18,7 +18,7 @@ #define PRIMAL_EVAL_INTEGRAL_HPP_ // Axom includes -#include "axom/config.hpp" // for compile-time configuration options +#include "axom/config.hpp" #include "axom/primal.hpp" #include "axom/primal/operators/detail/evaluate_integral_impl.hpp" @@ -33,11 +33,6 @@ // C++ includes #include -// Create type definitions the different field types (vector vs scalar) -// to allow for proper overloading of line integrals -using Point2D = axom::primal::Point; -using Vector2D = axom::primal::Vector; - namespace axom { namespace primal @@ -45,8 +40,9 @@ namespace primal /*! * \brief Evaluate a line integral along a collection of Bezier curves * - * The curve need not be connected, simply adding the portion along each segment. - * Uses Gaussian quadrature generated by MFEM. + * 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. @@ -68,11 +64,11 @@ double evaluate_line_integral(const axom::Array>& &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1)); double total_integral = 0.0; - for(int i = 0; i < cs.size(); i++) + for(const auto& curve : cs) { - // Compute the line integral along each component + // Compute the line integral along each component. total_integral += - detail::evaluate_line_integral_component(cs[i], integrand, quad); + detail::evaluate_line_integral_component(curve, integrand, quad); } return total_integral; @@ -133,34 +129,74 @@ double evaluate_area_integral(const axom::Array>& static mfem::IntegrationRules my_IntRules(mfem::Quadrature1D::GaussLegendre); const mfem::IntegrationRule *quad_Q, *quad_P; + if(npts_P <= 0) npts_P = npts_Q; + // Get the quadrature for the line integral. // Quadrature order is equal to 2*N - 1 quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); + quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1)); - // If the quadrature orders are the same (which they usually are), - // then only create one quadrature within the library and then - // the two pointers can point to the same object. Hopefully improves efficiency. - if(npts_P <= 0) - { - npts_P = npts_Q; - quad_P = quad_Q; - } - else - { - quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 1)); - } - - // Define lower bound as lowest control node for all curves - double int_lb = cs[0][0][1]; // Start as y-component of first node + // 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(int i = 0; i < cs.size(); i++) + 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); + const mfem::IntegrationRule *quad_Q, *quad_P; + + if(npts_P <= 0) npts_P = npts_Q; + + // Get the quadrature for the line integral. + // Quadrature order is equal to 2*N - 1 + quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); + 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(cs[i], + total_integral += detail::evaluate_area_integral_component(cpoly[i], integrand, int_lb, quad_Q, diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 65ba9b19be..aa3f872a2e 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -223,6 +223,6 @@ int main(int argc, char* argv[]) axom::slic::SimpleLogger logger; int result = RUN_ALL_TESTS(); - + return result; } From 6c02447bd34869a3f94018cc2e86ff1206075458 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Wed, 22 Jun 2022 17:33:21 -0700 Subject: [PATCH 14/22] update release notes --- RELEASE-NOTES.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md index d175c020f1..b62ffddeed 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`. From f4935d91591fe1570e14c10d1f03ef5442c0eec4 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 23 Jun 2022 08:53:19 -0700 Subject: [PATCH 15/22] Updated comments --- src/axom/primal/operators/evaluate_integral.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index cc0b218bc9..fa5c43ee0b 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -107,7 +107,6 @@ double evaluate_line_integral(const primal::BezierCurve& c, * * Assumes that the array of Bezier curves is closed and connected. Will compute * the integral regardless, but the result will be meaningless. - * Will eventually incorperate the CurverdPolygon class. * 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. * From b4bef28839a36a060f74aec9d53141e9d69da16e Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 23 Jun 2022 12:00:30 -0700 Subject: [PATCH 16/22] Added CurvedPolygon overload to evaluate_line_integral --- .../primal/operators/evaluate_integral.hpp | 35 +++++++++++++++++++ src/axom/primal/tests/primal_integral.cpp | 5 ++- 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index fa5c43ee0b..940f57ef8d 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -74,6 +74,41 @@ double evaluate_line_integral(const axom::Array>& 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. * diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index aa3f872a2e..6a14fcb561 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -112,10 +112,13 @@ TEST(primal_integral, evaluate_line_integral_scalar) Point2D {2.0, 4.0}}; Bezier parabola_segment(paranodes, 2); - // Compare against hand computed/high-precision calculated values + // 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); From 05349208890b86c67987ae23f02025a9ca8166db Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Thu, 23 Jun 2022 12:03:44 -0700 Subject: [PATCH 17/22] Updated to r-value reference for lambda functions --- src/axom/primal/operators/detail/evaluate_integral_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 2aa845d09c..17329e611a 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -39,7 +39,7 @@ namespace detail */ inline double evaluate_line_integral_component( const primal::BezierCurve& c, - std::function integrand, + std::function&& integrand, const mfem::IntegrationRule* quad) { // Store/compute quadrature result @@ -70,7 +70,7 @@ inline double evaluate_line_integral_component( */ inline double evaluate_line_integral_component( const primal::BezierCurve& c, - std::function vec_field, + std::function&& vec_field, const mfem::IntegrationRule* quad) { // Store/compute quadrature result From ddc9791926b8467df52a66be082a37a7e3e53afa Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 24 Jun 2022 09:00:59 -0700 Subject: [PATCH 18/22] Changed quadrature rule pointers to references --- src/axom/primal/CMakeLists.txt | 12 +++--- .../detail/evaluate_integral_impl.hpp | 38 +++++++++---------- .../primal/operators/evaluate_integral.hpp | 26 +++++++------ 3 files changed, 39 insertions(+), 37 deletions(-) diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index cf607de8ed..c917281e1c 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -68,12 +68,12 @@ 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 - ) +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 diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 17329e611a..37d240fb77 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -40,18 +40,18 @@ namespace detail inline double evaluate_line_integral_component( const primal::BezierCurve& c, std::function&& integrand, - const mfem::IntegrationRule* quad) + const mfem::IntegrationRule& quad) { // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad->GetNPoints(); q++) + 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); + 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(); + full_quadrature += quad.IntPoint(q).weight * integrand(x_q) * dx_q.norm(); } return full_quadrature; @@ -71,20 +71,20 @@ inline double evaluate_line_integral_component( inline double evaluate_line_integral_component( const primal::BezierCurve& c, std::function&& vec_field, - const mfem::IntegrationRule* quad) + const mfem::IntegrationRule& quad) { // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad->GetNPoints(); q++) + 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 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); + quad.IntPoint(q).weight * Vector2D::dot_product(func_val, dx_q); } return full_quadrature; @@ -111,32 +111,32 @@ 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) + 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++) + 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); + 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++) + 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}); + 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); + 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; + full_quadrature += quad_Q.IntPoint(q).weight * + c.dt(quad_Q.IntPoint(q).x)[0] * -antiderivative; } } diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 940f57ef8d..76ac0b4b2f 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -60,8 +60,8 @@ double evaluate_line_integral(const axom::Array>& // 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)); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); double total_integral = 0.0; for(const auto& curve : cs) @@ -95,8 +95,8 @@ double evaluate_line_integral(const primal::CurvedPolygon cpoly, // 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)); + 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++) @@ -131,8 +131,8 @@ double evaluate_line_integral(const primal::BezierCurve& c, // 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)); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); return detail::evaluate_line_integral_component(c, integrand, quad); } @@ -161,14 +161,15 @@ double evaluate_area_integral(const axom::Array>& // 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); - const mfem::IntegrationRule *quad_Q, *quad_P; if(npts_P <= 0) npts_P = npts_Q; // Get the quadrature for the line integral. // Quadrature order is equal to 2*N - 1 - quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); - quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 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]; @@ -211,14 +212,15 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, // 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); - const mfem::IntegrationRule *quad_Q, *quad_P; if(npts_P <= 0) npts_P = npts_Q; // Get the quadrature for the line integral. // Quadrature order is equal to 2*N - 1 - quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); - quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 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]; From 4d0d9e19b5ca6f8580a54d6868c30d545b6745b5 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 24 Jun 2022 09:00:59 -0700 Subject: [PATCH 19/22] Changed quadrature rule pointers to references --- src/axom/primal/CMakeLists.txt | 12 +++--- .../detail/evaluate_integral_impl.hpp | 38 +++++++++---------- .../primal/operators/evaluate_integral.hpp | 28 ++++++++------ 3 files changed, 41 insertions(+), 37 deletions(-) diff --git a/src/axom/primal/CMakeLists.txt b/src/axom/primal/CMakeLists.txt index cf607de8ed..c917281e1c 100644 --- a/src/axom/primal/CMakeLists.txt +++ b/src/axom/primal/CMakeLists.txt @@ -68,12 +68,12 @@ 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 - ) +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 diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 17329e611a..37d240fb77 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -40,18 +40,18 @@ namespace detail inline double evaluate_line_integral_component( const primal::BezierCurve& c, std::function&& integrand, - const mfem::IntegrationRule* quad) + const mfem::IntegrationRule& quad) { // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad->GetNPoints(); q++) + 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); + 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(); + full_quadrature += quad.IntPoint(q).weight * integrand(x_q) * dx_q.norm(); } return full_quadrature; @@ -71,20 +71,20 @@ inline double evaluate_line_integral_component( inline double evaluate_line_integral_component( const primal::BezierCurve& c, std::function&& vec_field, - const mfem::IntegrationRule* quad) + const mfem::IntegrationRule& quad) { // Store/compute quadrature result double full_quadrature = 0.0; - for(int q = 0; q < quad->GetNPoints(); q++) + 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 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); + quad.IntPoint(q).weight * Vector2D::dot_product(func_val, dx_q); } return full_quadrature; @@ -111,32 +111,32 @@ 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) + 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++) + 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); + 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++) + 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}); + 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); + 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; + full_quadrature += quad_Q.IntPoint(q).weight * + c.dt(quad_Q.IntPoint(q).x)[0] * -antiderivative; } } diff --git a/src/axom/primal/operators/evaluate_integral.hpp b/src/axom/primal/operators/evaluate_integral.hpp index 940f57ef8d..12d22a8a94 100644 --- a/src/axom/primal/operators/evaluate_integral.hpp +++ b/src/axom/primal/operators/evaluate_integral.hpp @@ -12,6 +12,8 @@ * 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_ @@ -60,8 +62,8 @@ double evaluate_line_integral(const axom::Array>& // 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)); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); double total_integral = 0.0; for(const auto& curve : cs) @@ -95,8 +97,8 @@ double evaluate_line_integral(const primal::CurvedPolygon cpoly, // 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)); + 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++) @@ -131,8 +133,8 @@ double evaluate_line_integral(const primal::BezierCurve& c, // 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)); + const mfem::IntegrationRule& quad = + my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts - 1); return detail::evaluate_line_integral_component(c, integrand, quad); } @@ -161,14 +163,15 @@ double evaluate_area_integral(const axom::Array>& // 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); - const mfem::IntegrationRule *quad_Q, *quad_P; if(npts_P <= 0) npts_P = npts_Q; // Get the quadrature for the line integral. // Quadrature order is equal to 2*N - 1 - quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); - quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 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]; @@ -211,14 +214,15 @@ double evaluate_area_integral(const primal::CurvedPolygon cpoly, // 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); - const mfem::IntegrationRule *quad_Q, *quad_P; if(npts_P <= 0) npts_P = npts_Q; // Get the quadrature for the line integral. // Quadrature order is equal to 2*N - 1 - quad_Q = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_Q - 1)); - quad_P = &(my_IntRules.Get(mfem::Geometry::SEGMENT, 2 * npts_P - 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]; From 702de762d5a6a3a52f5f2e8917598c447267a313 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 24 Jun 2022 09:57:28 -0700 Subject: [PATCH 20/22] Temporarily removed R-value references --- src/axom/primal/operators/detail/evaluate_integral_impl.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp index 37d240fb77..c5934e6adc 100644 --- a/src/axom/primal/operators/detail/evaluate_integral_impl.hpp +++ b/src/axom/primal/operators/detail/evaluate_integral_impl.hpp @@ -39,7 +39,7 @@ namespace detail */ inline double evaluate_line_integral_component( const primal::BezierCurve& c, - std::function&& integrand, + std::function integrand, const mfem::IntegrationRule& quad) { // Store/compute quadrature result @@ -70,7 +70,7 @@ inline double evaluate_line_integral_component( */ inline double evaluate_line_integral_component( const primal::BezierCurve& c, - std::function&& vec_field, + std::function vec_field, const mfem::IntegrationRule& quad) { // Store/compute quadrature result From de51f54ac67ecfd1acb800b68325801adb4b282f Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Fri, 24 Jun 2022 10:37:19 -0700 Subject: [PATCH 21/22] Add tests for line integrals on disconnected curves --- src/axom/primal/tests/primal_integral.cpp | 34 +++++++++++++++++++++-- 1 file changed, 32 insertions(+), 2 deletions(-) diff --git a/src/axom/primal/tests/primal_integral.cpp b/src/axom/primal/tests/primal_integral.cpp index 6a14fcb561..cb6bd7924a 100644 --- a/src/axom/primal/tests/primal_integral.cpp +++ b/src/axom/primal/tests/primal_integral.cpp @@ -143,8 +143,6 @@ TEST(primal_integral, evaluate_line_integral_scalar) axom::Array connected_curve( {cubic_segment, linear_segment, quadratic_segment}); - - // Compare against hand computed/high-precision calculated values EXPECT_NEAR(evaluate_line_integral(connected_curve, const_integrand, npts), 8.28968500196, abs_tol); @@ -154,6 +152,18 @@ TEST(primal_integral, evaluate_line_integral_scalar) 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) @@ -217,6 +227,26 @@ TEST(primal_integral, evaluate_line_integral_vector) 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[]) From 8f75254819b2ef9cb399f8e767924c7f48297bc6 Mon Sep 17 00:00:00 2001 From: Jacob Spainhour Date: Mon, 27 Jun 2022 12:59:01 -0700 Subject: [PATCH 22/22] Added extra dependencies for Primal --- src/docs/sphinx/quickstart_guide/config_build.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docs/sphinx/quickstart_guide/config_build.rst b/src/docs/sphinx/quickstart_guide/config_build.rst index 502c9e23f1..387c9002fc 100644 --- a/src/docs/sphinx/quickstart_guide/config_build.rst +++ b/src/docs/sphinx/quickstart_guide/config_build.rst @@ -59,7 +59,7 @@ Library Dependencies `c2c`_ Optional C2C_DIR `HDF5`_ Optional: Sidre HDF5_DIR `Lua`_ Optional: Inlet LUA_DIR - `MFEM`_ Optional: Quest 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