Skip to content

Conversation

@jcs15c
Copy link
Contributor

@jcs15c jcs15c commented Dec 15, 2025

Summary

This PR is an enhancement, adding an additional primal::evaluate_line_integral method to evaluate_integral.hpp which permits integral evaluation over any lambda function which has an "integrable" return type, i.e. one for which addition and scalar multiplication are defined. Most prominently, this permits evaluation over integrands which return primal::Vector objects. The same functionality is added to primal::evaluate_area_integral. These methods are supported by additional tests in primal_integral.cpp.

More generally, this addition is meant to simplify the process of evaluating physical quantities for geometric objects. For example, each component of the average surface normal of a 3D NURBS surface is computed by integrating each component of the surface normal vector over the surface. Evaluating this as three separate integrals results in a significant amount of redundant computation which can now be avoided.

A note on names

Implementing this feature necessarily introduces some amount of naming ambiguity between methods in evaluate_integral.hpp, as now there are two types of "vector field integrals" with distinct meanings. Currently defined in develop are methods evaluate_vector_line_integral, which evaluates integrals of the form

$$\int_C \vec{F}(\vec{x}) \cdot d\vec{r},$$

which assume that $\vec{F}$ : $\mathbb{R}^n \to \mathbb{R}^n$. Typically $n$ is equal to 2 or 3, but technically other values are possible mathematically. Either way, this method always returns a scalar value, i.e. a double, as would be suggested by the dot product in the formulation.

The proposed changes introduces methods of evaluating integrals of the form

$$\int_C \vec{F}(\vec{x}) dr,$$

which assumes that $\vec{F}$ : $\mathbb{R}^n \to \mathbb{R}^m$, with an output in $\mathbb{R}^m$. While $n$ will most often be 2 or 3, $m$ can be significantly larger. Specifically, I'm targeting a future use with $m = 39$, corresponding to a function with a return type primal::Vector<double, 39>. These could also be characterized as "vector field integrals," even though they behave more like "component-wise scalar field integrals." This ambiguity is essentially unavoidable, because even the input for "scalar field integrals" ($\mathbb{R}^n$) are themselves a "vector space over a scalar field" in a mathematical sense.

To make sure any ambiguity is inherently resolved by templates, the functionality of evaluate_scalar_line_integral methods in evaluate_integral.hpp is merged into an evaluate_line_integral method which performs this more general component-wise integration. Conveniently, this mirrors the evaluate_area_integral method.

@jcs15c jcs15c self-assigned this Dec 15, 2025
@jcs15c jcs15c added enhancement New feature or request Primal Issues related to Axom's 'primal component labels Dec 15, 2025
@jcs15c jcs15c requested review from gunney1 and rhornung67 December 15, 2025 22:05
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 -- nice contribution!

namespace internal
{
///@{
/// \name Boilerplate to support integrals of functions with general return types,
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
/// \name Boilerplate to support integrals of functions with general return types,
/// \name Type traits to support integrals of functions with general return types,

int_lb,
npts_Q,
npts_P);
total_integral += detail::evaluate_area_integral_component(beziers[j],
Copy link
Member

Choose a reason for hiding this comment

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

Looks like this needs to be updated:

Suggested change
total_integral += detail::evaluate_area_integral_component(beziers[j],
total_integral += detail::evaluate_area_integral_component(bez,

Copy link
Member

Choose a reason for hiding this comment

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

Thanks for the other changes!

// Integrate the surface normal over the patches
ret_vec += evaluate_area_integral(
nPatch.getTrimmingCurves(),
[&nPatch](Point2D x) -> Vector<T, 3> { return nPatch.normal(x[0], x[1]); } npts);
Copy link
Member

Choose a reason for hiding this comment

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

Where did the comma in front of "npts" go? Isn't it still needed?

ret_vec += evaluate_area_integral(
nPatch.getTrimmingCurves(),
[&nPatch](Point2D x) -> Vector<T, 3> { return nPatch.normal(x[0], x[1]); } npts);
};
Copy link
Member

Choose a reason for hiding this comment

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

Extra semicolon. Delete it.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

This whole function call is a mess, apparently. Sorry about that, I'll fix it and double check the rest.

Copy link
Member

@BradWhitlock BradWhitlock left a comment

Choose a reason for hiding this comment

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

Thanks!

@kennyweiss kennyweiss self-requested a review December 23, 2025 17:52
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 for the updates @jcs15c -- looks great!

EXPECT_NEAR(evaluate_scalar_line_integral(ellipse_arc, transc_integrand, npts),
1.38837959326,
abs_tol);
EXPECT_NEAR(evaluate_line_integral(ellipse_arc, const_integrand, npts), 2.42211205514, 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 update the RELEASE-NOTES w/ the changed function name

template <typename Lambda,
typename T,
int NDIMS,
typename LambdaRetTyp = std::invoke_result_t<Lambda, typename BezierCurve<T, NDIMS>::PointType>>
Copy link
Contributor

Choose a reason for hiding this comment

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

Throughout this file, sometimes I read LambdaRetType (as in (new) line 87 in the Doxygen), sometimes LambdaRetTyp (as here), and sometimes LambdaRetTYpe (as at (new) line 300). Please make them all consistent. If they have to be disambiguated, please make it a larger string or conceptual distance.

Copy link
Contributor

@Arlie-Capps Arlie-Capps 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 !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request Primal Issues related to Axom's 'primal component

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants