TrapezoidalTorus.hpp File Reference

Trapezoidal product quadrature on a torus, with the point counts fixed at compile time or run time. More...

Classes
class	mundy::TrapezoidalTorus< Scalar, N, M, AspectRatio >
	The N x M trapezoidal rule on the torus with minor radius 1 and major radius AspectRatio. More...

Namespaces
namespace	mundy

Functions
template<class Scalar>
void	mundy::trapezoidal_torus_rule (unsigned n, unsigned m, double aspect_ratio, std::vector< Scalar > &points, std::vector< Scalar > &weights)
	The n x m trapezoidal rule on the torus with minor radius 1 and major radius aspect_ratio, into std::vectors.
template<class PointsView, class WeightsView>
void	mundy::trapezoidal_torus_rule (unsigned n, unsigned m, double aspect_ratio, PointsView &points, WeightsView &weights)
	The n x m trapezoidal rule on the torus with minor radius 1 and major radius aspect_ratio, into rank-1 views.

Detailed Description

The torus has minor radius 1 and major radius a > 1, its aspect ratio: the points ((a + cos theta) cos phi, (a + cos theta) sin phi, sin theta) for the poloidal angle theta around the tube and the toroidal angle phi around the z axis, with area element (a + cos theta) dtheta dphi. The rule is the trapezoidal rule in both angles: N equispaced theta times M equispaced phi, N M points in all. Both angles are periodic, so the rule converges spectrally for smooth integrands, and it is exact for every polynomial in (x, y, z) of degree at most min(N - 2, M - 1).

Points are (x, y, z) triples, ring by ring in theta = 2 pi i / N from the outer equator (theta = 0) over the top, and by increasing phi = 2 pi j / M from +x within a ring. For a torus with minor radius r, scale the points by r and the weights by r^2; the outward unit normal at a point is (cos theta cos phi, cos theta sin phi, sin theta). Points and weights are correctly rounded to Scalar.

TrapezoidalTorus uses compile-time work that grows like N M. At the time of writing, GCC 13.3.0 compiles for N = M <= 238 whereas clang 18.1.8 compiles for N = M <= 64. Larger rules need -fconstexpr-ops-limit (GCC) or -fconstexpr-steps (clang). Its table holds 4 N M scalars, so for large rules or on the device prefer trapezoidal_torus_rule(n, m, aspect_ratio, points, weights), which builds the same rule, bit for bit, at run time.