GaussTriangleEllipsoid.hpp File Reference

Three-point Gauss triangle quadrature on an ellipsoid meshed from a refined octahedron or icosahedron, with the refinement level fixed at compile time or run time. More...

Classes
class	mundy::GaussTriangleEllipsoid< Scalar, Polyhedron, L, BOverA, COverA >
	The three-point Gauss triangle rule on the Polyhedron refined L times and stretched onto the ellipsoid with semi-axes 1, BOverA, COverA, with outward unit normals. More...

Namespaces
namespace	mundy

Enumerations
enum class	mundy::SpherePolyhedron : unsigned { mundy::octahedron = 8 , mundy::icosahedron = 20 }
	The polyhedron whose refined faces mesh the sphere. Each enumerator's value is the polyhedron's face count. More...

Functions
template<class Scalar>
void	mundy::gauss_triangle_ellipsoid_rule (SpherePolyhedron polyhedron, unsigned l, double b_over_a, double c_over_a, std::vector< Scalar > &points, std::vector< Scalar > &normals, std::vector< Scalar > &weights)
	The three-point Gauss triangle rule on the polyhedron refined l times and stretched onto the ellipsoid with semi-axes 1, b_over_a, c_over_a, into std::vectors.
template<class PointsView, class NormalsView, class WeightsView>
void	mundy::gauss_triangle_ellipsoid_rule (SpherePolyhedron polyhedron, unsigned l, double b_over_a, double c_over_a, PointsView &points, NormalsView &normals, WeightsView &weights)
	The three-point Gauss triangle rule on the polyhedron refined l times and stretched onto the ellipsoid with semi-axes 1, b_over_a, c_over_a, into rank-1 views.

Detailed Description

The ellipsoid is meshed as in C. Pozrikidis's BEMLIB (trgl6_octa.f, trgl6_icos.f). Each face of the polyhedron is split into four L times with every new edge midpoint projected onto the unit sphere, and the resulting F 4^L six-node triangles (F = 8 or 20 faces) are stretched by (x, y, z) -> (x, (b/a) y, (c/a) z) onto the ellipsoid with semi-axes 1, b/a, c/a. Each element is the quadratic patch through its six nodes, with each edge node placed along its reference edge in proportion to its distances from the edge's ends (abc.f). The rule is the three-point Gauss rule on every patch (gauss_trgl.f), 3 F 4^L points in all. It integrates over the patches rather than the ellipsoid and converges at fourth order in the mesh spacing for smooth integrands. With b/a = c/a = 1 the ellipsoid is the unit sphere.

Points are (x, y, z) triples, element by element in BEMLIB's order, and within an element nearest its vertices 0, 1, 2 in turn. The points lie on the patches, within O(h^3) of the ellipsoid (O(h^4) of the sphere), so the normals are the patches' outward unit normals. For semi-axes a, b, c, scale the points by a and the weights by a^2; BEMLIB's volume-equivalent radius r gives a = r / (b/a c/a)^(1/3). Points, normals, and weights are correctly rounded to Scalar.

GaussTriangleEllipsoid uses compile-time work that grows like 4^L. At the time of writing, GCC 13.3.0 compiles for L <= 2 (octahedron) and L <= 1 (icosahedron) whereas clang 18.1.8 compiles for L = 0. Larger L need -fconstexpr-ops-limit (GCC) or -fconstexpr-steps (clang). Its tables hold 21 F 4^L scalars, so for large L or on the device prefer gauss_triangle_ellipsoid_rule(polyhedron, l, b_over_a, c_over_a, points, normals, weights), which builds the same rule, bit for bit, at run time.