Skip to content

Conversation

@jcs15c
Copy link
Contributor

@jcs15c jcs15c commented Jun 21, 2022

  • This PR is a feature
  • It adds the evaluate_integrals.hpp file to primal/operators, which contains methods to
    • Evaluate scalar line integrals on collections of Bezier curves.
    • Evaluate vector field line integrals on collections of Bezier curves.
    • Evaluate 2D area integrals for regions defined by Bezier curves.
  • It adds additional tests in primal for the above methods.
  • It adds an optional dependency for primal on mfem, used to evaluate quadrature algorithms.

edit: See https://en.wikipedia.org/wiki/Line_integral#Line_integral_of_a_scalar_field and https://en.wikipedia.org/wiki/Line_integral#Line_integral_of_a_vector_field for information about the different types of line integrals.

@jcs15c jcs15c added enhancement New feature or request Primal Issues related to Axom's 'primal component labels Jun 21, 2022
@jcs15c jcs15c requested a review from kennyweiss June 21, 2022 19:27
@jcs15c jcs15c self-assigned this Jun 21, 2022
@jcs15c jcs15c removed the enhancement New feature or request label Jun 21, 2022
Copy link
Member

@kennyweiss kennyweiss left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jcs15c -- this is great!

I added some suggestions for minor changes.
Please also add a line to the RELEASE-NOTES.

@jcs15c jcs15c force-pushed the feature/spainhour/bezier_quadrature branch from 3f5d834 to b8a0599 Compare June 22, 2022 23:50
@jcs15c jcs15c force-pushed the feature/spainhour/bezier_quadrature branch from b8a0599 to 4995da1 Compare June 23, 2022 00:14
Copy link
Member

@kennyweiss kennyweiss left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @jcs15c -- this is really nice!

Tag: @davidgunderman

*/
inline double evaluate_line_integral_component(
const primal::BezierCurve<double, 2>& c,
std::function<double(Point2D)> integrand,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the lambda should be passed by r-value ref (here and elsewhere):

Suggested change
std::function<double(Point2D)> integrand,
std::function<double(Point2D)>&& integrand,

(@kanye-quest -- could you please confirm?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm very unclear on the details, but passing the std::function's by R-value reference is causing some issues with some to-be-implemented winding number code. Specifically, if we pass to evaluate_line_integral a lambda function that is the return value for some other function (see below), then evaluate_line_integral_component can no longer fit the argument to either overloaded definition. This doesn't happen when the std::function's are passed by constant reference, or when the lambda function is not the return value for another function, so I am unsure how to keep this as general as possible.

// An example of function that returns argument to evaluate_line_integral
std::function<Vector2D(Point2D)> get_vector_field(Point2D p) 
{
  return [p](Point2D x) -> Vector2D { return Vector2D({p[1], p[0]}); };
}

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks -- ok, sounds like something to think about and try to fix in the future.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This might have something to do with this function (and the one below) accepting an std::function instead of a Lambda&& like in the other methods. The latter is a universal reference, so it can bind to both lvalues and rvalues.

For what it's worth, I think accepting Lambda&& might be better from an inlineability perspective -- it would be much more difficult for the compiler to inline through a type-erased std::function.

@kennyweiss kennyweiss requested a review from gunney1 June 23, 2022 15:46
Copy link
Contributor

@publixsubfan publixsubfan left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM, and thanks for the tests! Just a few comments to think about.

ELEMENTS operators/evaluate_integral.hpp
operators/detail/evaluate_integral_impl.hpp
IF MFEM_FOUND
)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: alignment?

*/
inline double evaluate_line_integral_component(
const primal::BezierCurve<double, 2>& c,
std::function<double(Point2D)> integrand,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This might have something to do with this function (and the one below) accepting an std::function instead of a Lambda&& like in the other methods. The latter is a universal reference, so it can bind to both lvalues and rvalues.

For what it's worth, I think accepting Lambda&& might be better from an inlineability perspective -- it would be much more difficult for the compiler to inline through a type-erased std::function.

// 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));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor: any reason to use pointers to the IntegrationRules instead of references?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only that the quest test I used for reference used them, changing it to references makes it much cleaner. Thank you!

Comment on lines 14 to 15
* Regions Bounded by Rational Parametric Curves" by David Gunderman et al.
*/
Copy link
Member

@agcapps agcapps Jun 24, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I appreciate this whole comment! I see you wrote "with ... rules supplied by MFEM," as is accurate. Please add the additional note, to make things extra clear:

Suggested change
* Regions Bounded by Rational Parametric Curves" by David Gunderman et al.
*/
* Regions Bounded by Rational Parametric Curves" by David Gunderman et al.
*
* \note This requires the MFEM third-party library.
*/

This Doxygen mirrors the compile-time checks around line 28 of this file and line 15 of evaluate_integral_impl.hpp

EXPECT_NEAR(evaluate_line_integral(connected_curve, transc_integrand, npts),
-0.574992518405,
abs_tol);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a test with disconnected line segments, as this is noted as possible in the doxygen for evaluate_line_integral. Perhaps displace the middle, linear segment by (0, 2) or some such.

EXPECT_NEAR(evaluate_line_integral(parabola_shape, winding_field, npts),
1.0,
abs_tol);
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a test for disconnected segments.

Copy link
Member

@agcapps agcapps left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for this!

@rhornung67 rhornung67 added this to the March 2022 Release milestone Jun 27, 2022
@jcs15c jcs15c merged commit 4acbe95 into develop Jun 28, 2022
@jcs15c jcs15c deleted the feature/spainhour/bezier_quadrature branch June 28, 2022 18:47
@jcs15c jcs15c restored the feature/spainhour/bezier_quadrature branch August 4, 2022 16:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Primal Issues related to Axom's 'primal component

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants