pbat::math Namespace Reference

Math related functionality. More...

Namespaces
namespace	linalg
	Linear Algebra related functionality.
namespace	optimization
	Namespace for optimization algorithms.
namespace	polynomial
	The namespace for the Polynomial module.

Classes
struct	DynamicQuadrature
	Quadrature rule with variable points and weights. More...
struct	FixedSizeVariableQuadrature
	Represents a quadrature scheme that can be constructed via existing quadrature schemes. More...
class	LinearOperator
	Zero-overhead composite type satisfying the CLinearOperator concept. More...
struct	OverflowChecked
	Wrapper around integer types that throws when integer overflow is detected. More...
struct	Rational
	Fixed size rational number \( \frac{a}{b} \) using std::int64_t for numerator and denominator. More...

Concepts
concept	CQuadratureRule
concept	CFixedPointQuadratureRule
concept	CPolynomialQuadratureRule
concept	CFixedPointPolynomialQuadratureRule
concept	CLinearOperator
	Concept for operator that satisfies linearity in the mathematical sense.

Typedefs
template<int Dims, int Order, common::CFloatingPoint TScalar = Scalar>
using	GaussLegendreQuadrature = typename detail::GaussLegendreQuadrature<Dims, Order, TScalar>
	Shifted Gauss-Legendre quadrature scheme over the unit box \( 0 \leq X \leq 1 \; \text{s.t.} \; X \in \mathbb{R}^d \) in dimensions \( d=1,2,3 \).
template<int Dims, int Order, common::CFloatingPoint TScalar = Scalar>
using	SymmetricSimplexPolynomialQuadratureRule
	Symmetric quadrature rule for reference simplices in dimensions \( d=1,2,3 \).

Functions
template<std::integral Integer>
bool	AddOverflows (Integer a, Integer b)
	Checks if the operation \( a+b \) is not in the range of values representable by the type Integer.
template<std::integral Integer>
bool	MultiplyOverflows (Integer a, Integer b)
	Checks if the operation \( ab \) is not in the range of values representable by the type of Integer.
template<std::integral Integer>
bool	NegationOverflows (Integer a)
	Checks if the operation \( -a \) is not in the range of values representable by the type of Integer.
template<CLinearOperator... TLinearOperators>
LinearOperator< TLinearOperators... >	ComposeLinearOperators (TLinearOperators const &... inOps)
	Compose multiple linear operators into a single linear operator.
template<polynomial::CBasis TBasis, CFixedPointPolynomialQuadratureRule TQuad>
Matrix< TBasis::kSize, TQuad::kPoints >	ReferenceMomentFittingMatrix (TBasis const &Pb, TQuad const &Q)
	Compute the moment fitting matrix in the reference simplex space.
template<polynomial::CBasis TBasis, CPolynomialQuadratureRule TQuad>
Matrix< TBasis::kSize, Eigen::Dynamic >	ReferenceMomentFittingMatrix (TBasis const &Pb, TQuad const &Q)
	Compute the moment fitting matrix in the reference simplex space.
template<polynomial::CBasis TBasis, class TDerivedXg>
Matrix< TBasis::kSize, Eigen::Dynamic >	ReferenceMomentFittingMatrix (TBasis const &Pb, Eigen::MatrixBase< TDerivedXg > const &Xg)
	Compute the moment fitting matrix in the reference simplex space.
template<class TDerivedP, class TDerivedB>
VectorX	MomentFittedWeights (Eigen::MatrixBase< TDerivedP > const &P, Eigen::DenseBase< TDerivedB > const &b, Index maxIterations=10, Scalar precision=std::numeric_limits< Scalar >::epsilon())
	Computes non-negative quadrature weights \( w \) by moment fitting.
template<polynomial::CBasis Polynomial, class TDerivedXg, class TDerivedWg>
Vector< Polynomial::kSize >	Integrate (Polynomial const &P, Eigen::MatrixBase< TDerivedXg > const &Xg, Eigen::DenseBase< TDerivedWg > const &wg)
	Computes \( \int_{\Omega} P(x) \, d\Omega \) given an existing quadrature rule \((Xg,wg) \).
template<polynomial::CBasis Polynomial, class TDerivedXg1, class TDerivedXg2, class TDerivedWg2>
Vector< TDerivedXg1::ColsAtCompileTime >	TransferQuadrature (Polynomial const &P, Eigen::MatrixBase< TDerivedXg1 > const &Xg1, Eigen::MatrixBase< TDerivedXg2 > const &Xg2, Eigen::DenseBase< TDerivedWg2 > const &wg2, Index maxIterations=10, Scalar precision=std::numeric_limits< Scalar >::epsilon())
	Obtain weights \( w_g^1 \) by transferring an existing quadrature rule \( X_g^2, w_g^2 \) defined on a simplex onto a new quadrature rule \( X_g^1, w_g^1 \) defined on the same simplex.
template<auto Order, class TDerivedS1, class TDerivedXi1, class TDerivedS2, class TDerivedXi2, class TDerivedWg2>
std::pair< VectorX, VectorX >	TransferQuadrature (Eigen::DenseBase< TDerivedS1 > const &S1, Eigen::MatrixBase< TDerivedXi1 > const &Xi1, Eigen::DenseBase< TDerivedS2 > const &S2, Eigen::MatrixBase< TDerivedXi2 > const &Xi2, Eigen::DenseBase< TDerivedWg2 > const &wi2, Index nSimplices=-1, bool bEvaluateError=false, Index maxIterations=10, Scalar precision=std::numeric_limits< Scalar >::epsilon())
	Obtain weights \( w_g^1 \) by transferring an existing quadrature rule \( X_g^2, w_g^2 \) defined on a domain of simplices onto a new quadrature rule \( X_g^1, w_g^1 \) defined on the same domain.
template<int Order, class TDerivedS1, class TDerivedX1, class TDerivedS2, class TDerivedX2, class TDerivedW2>
std::tuple< MatrixX, MatrixX, IndexVectorX >	ReferenceMomentFittingSystems (Eigen::DenseBase< TDerivedS1 > const &S1, Eigen::MatrixBase< TDerivedX1 > const &X1, Eigen::DenseBase< TDerivedS2 > const &S2, Eigen::MatrixBase< TDerivedX2 > const &X2, Eigen::DenseBase< TDerivedW2 > const &w2, Index nSimplices=Index(-1))
	Computes reference moment fitting systems for all simplices \( S_1 \), given an existing quadrature rule \( (X_g^2, w_g^2) \) defined on a domain of simplices. The quadrature points of the new rule \( X_g^1 \) are given and fixed.
template<class TDerivedM, class TDerivedP>
CSRMatrix	BlockDiagonalReferenceMomentFittingSystem (Eigen::MatrixBase< TDerivedM > const &M, Eigen::DenseBase< TDerivedP > const &P)
	Block diagonalizes the given reference moment fitting systems (see ReferenceMomentFittingSystems()) into a large sparse matrix.
template<std::integral Integer>
Rational	operator- (Integer a, Rational const &b)
	Subtraction operation between Rational and integral type.
template<std::integral Integer>
Rational	operator+ (Integer a, Rational const &b)
	Addition operation between Rational and integral type.
template<std::integral Integer>
Rational	operator* (Integer a, Rational const &b)
	Multiplication operation between Rational and integral type.
template<std::integral Integer>
Rational	operator/ (Integer a, Rational const &b)
	Division operation between Rational and integral type.

Detailed Description

Math related functionality.

Typedef Documentation

GaussLegendreQuadrature

template<int Dims, int Order, common::CFloatingPoint TScalar = Scalar>

using pbat::math::GaussLegendreQuadrature = typename detail::GaussLegendreQuadrature<Dims, Order, TScalar>

Shifted Gauss-Legendre quadrature scheme over the unit box \( 0 \leq X \leq 1 \; \text{s.t.} \; X \in \mathbb{R}^d \) in dimensions \( d=1,2,3 \).

The points are specified as \( d+1 \) coordinate tuples in affine coordinates, i.e. the first coordinate \( X_0 = 1 - \sum_{i=1}^d X_i \).

GaussLegendreQuadrature instantiations expose static members

template <int Dims, int Order>
struct GaussLegendreQuadrature
{
    inline static std::uint8_t constexpr kDims;
    inline static std::uint8_t constexpr kOrder;
    inline static int constexpr kPoints;
    inline static std::array<TScalar, (kDims + 1) * kPoints> constexpr points;
    inline static std::array<TScalar, kPoints> constexpr weights;
};

Example usage:

using Quad = GaussLegendreQuadrature<3, 2>;
auto Xg = pbat::common::ToEigen(Quad::points).reshaped(Quad::kDims + 1, Quad::kPoints);
auto wg = pbat::common::ToEigen(Quad::weights);
for (auto g = 0; g < Quad::kPoints; ++g)
{
    auto X = Xg.col(g).segment(1, Quad::kDims);
    integrand += wg(g)*f(X);
}

Template Parameters

Dims	Spatial dimensions
Order	Polynomial quadrature order
TScalar	Floating point type, defaults to Scalar

SymmetricSimplexPolynomialQuadratureRule

template<int Dims, int Order, common::CFloatingPoint TScalar = Scalar>

using pbat::math::SymmetricSimplexPolynomialQuadratureRule

Initial value:

 
    typename detail::SymmetricSimplexPolynomialQuadratureRule<Dims, Order, TScalar>

Symmetric quadrature rule for reference simplices in dimensions \( d=1,2,3 \).

The points are specified as \( d+1 \) coordinate tuples in affine coordinates, i.e. the first coordinate \( X_0 = 1 - \sum_{i=1}^d X_i \).

Every point belongs to the reference simplex, e.g.

• the line segment \( 0,1 \) in 1D,
• the triangle \( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \end{pmatrix} \) in 2D, and
• the tetrahedron \( \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 1\end{pmatrix} \) in 3D.

SymmetricSimplexPolynomialQuadratureRule instantiations expose static members

template <int Dims, int Order>
struct SymmetricSimplexPolynomialQuadratureRule
{
    inline static std::uint8_t constexpr kDims;
    inline static std::uint8_t constexpr kOrder;
    inline static int constexpr kPoints;
    inline static std::array<TScalar, (kDims + 1) * kPoints> constexpr points;
    inline static std::array<TScalar, kPoints> constexpr weights;
};

Example usage:

using Quad = SymmetricSimplexPolynomialQuadratureRule<3, 2>;
auto Xg = pbat::common::ToEigen(Quad::points).reshaped(Quad::kDims + 1, Quad::kPoints);
auto wg = pbat::common::ToEigen(Quad::weights);
for (auto g = 0; g < Quad::kPoints; ++g)
{
    auto X = Xg.col(g).segment(1, Quad::kDims);
    integrand += wg(g)*f(X);
}

Template Parameters

Dims	Dimension of the reference simplex
Order	Polynomial order of the quadrature rule
TScalar	Scalar type of the quadrature rule, defaults to Scalar

Function Documentation

AddOverflows()

template<std::integral Integer>

bool pbat::math::AddOverflows	(	Integer	a,
		Integer	b )

Checks if the operation \( a+b \) is not in the range of values representable by the type Integer.

Template Parameters

Integer

The type of the integers to check for overflow.

Parameters

a	Left operand
b	Right operand

Returns

True if the operation overflows, false otherwise.

BlockDiagonalReferenceMomentFittingSystem()

template<class TDerivedM, class TDerivedP>

CSRMatrix pbat::math::BlockDiagonalReferenceMomentFittingSystem	(	Eigen::MatrixBase< TDerivedM > const &	M,
		Eigen::DenseBase< TDerivedP > const &	P )

Block diagonalizes the given reference moment fitting systems (see ReferenceMomentFittingSystems()) into a large sparse matrix.

Template Parameters

TDerivedM	Eigen matrix expression for moment fitting matrices
TDerivedP	Eigen vector expression for prefix into columns of moment fitting matrices

Parameters

M	Moment fitting matrices
P	Prefix into columns of moment fitting matrices

Returns

The block diagonal row sparse matrix \( \mathbf{A} \), whose diagonal blocks are the individual reference moment fitting matrices in M, such that \( \mathbf{A} \mathbf{w} = \text{vec}(\mathbf{B}) \) is the global sparse linear system to solve for quadrature weights \(\mathbf{w} \).

ComposeLinearOperators()

template<CLinearOperator... TLinearOperators>

LinearOperator< TLinearOperators... > pbat::math::ComposeLinearOperators

(

TLinearOperators const &...

inOps

)

Compose multiple linear operators into a single linear operator.

Refer to LinearOperator.cpp for usage examples.

Template Parameters

TLinearOperators

Linear operator parameter types

Parameters

inOps

Linear operators to be composed

Returns

Composed linear operator

Integrate()

template<polynomial::CBasis Polynomial, class TDerivedXg, class TDerivedWg>

Vector< Polynomial::kSize > pbat::math::Integrate	(	Polynomial const &	P,
		Eigen::MatrixBase< TDerivedXg > const &	Xg,
		Eigen::DenseBase< TDerivedWg > const &	wg )

Computes \( \int_{\Omega} P(x) \, d\Omega \) given an existing quadrature rule \((Xg,wg) \).

Template Parameters

Polynomial	Polynomial basis type
TDerivedXg	Eigen matrix expression for quadrature points
TDerivedWg	Eigen vector expression for quadrature weights

Parameters

P	Polynomial basis
Xg	\( d \times n \) array of quadrature point positions defined in reference space
wg	\( n \times 1 \) array of quadrature weights

Returns

\( s \times 1 \) array of integrated polynomials

MomentFittedWeights()

template<class TDerivedP, class TDerivedB>

VectorX pbat::math::MomentFittedWeights	(	Eigen::MatrixBase< TDerivedP > const &	P,
		Eigen::DenseBase< TDerivedB > const &	b,
		Index	maxIterations = 10,
		Scalar	precision = std::numeric_limits<Scalar>::epsilon() )

Computes non-negative quadrature weights \( w \) by moment fitting.

Template Parameters

TDerivedP	Eigen matrix expression for moment fitting matrix
TDerivedB	Eigen vector expression for target integrated polynomials

Parameters

P	Moment fitting matrix \( \mathbf{P} \in \mathbb{R}^{s \times n} \), where \( s \) is the polynomial basis' size and \( n \) is the number of quadrature points
b	Target integrated polynomials \( \mathbf{b} \in \mathbb{R}^s \)
maxIterations	Maximum number of non-negative least-squares active set solver
precision	Convergence threshold

Returns

Non-negative quadrature weights \( \mathbf{w} \in \mathbb{R}^n \)

MultiplyOverflows()

template<std::integral Integer>

bool pbat::math::MultiplyOverflows	(	Integer	a,
		Integer	b )

Checks if the operation \( ab \) is not in the range of values representable by the type of Integer.

Template Parameters

Integer

The type of the integers to check for overflow.

Parameters

a	Left operand
b	Right operand

Returns

True if the operation overflows, false otherwise.

NegationOverflows()

template<std::integral Integer>

bool pbat::math::NegationOverflows

(

Integer

a

)

Checks if the operation \( -a \) is not in the range of values representable by the type of Integer.

Template Parameters

Integer

The type of the integers to check for overflow.

Parameters

a

Operand

Returns

True if the operation overflows, false otherwise.

operator*()

template<std::integral Integer>

Rational pbat::math::operator* ( Integer a, Rational const & b )

inline

Multiplication operation between Rational and integral type.

Template Parameters

Integer

Integral type

Parameters

a	Left operand
b	Right operand

Returns

Result of multiplication

operator+()

template<std::integral Integer>

Rational pbat::math::operator+ ( Integer a, Rational const & b )

inline

Addition operation between Rational and integral type.

Template Parameters

Integer

Integral type

Parameters

a	Left operand
b	Right operand

Returns

Result of addition

operator-()

template<std::integral Integer>

Rational pbat::math::operator- ( Integer a, Rational const & b )

inline

Subtraction operation between Rational and integral type.

Template Parameters

Integer

Integral type

Parameters

a	Left operand
b	Right operand

Returns

Result of subtraction

operator/()

template<std::integral Integer>

Rational pbat::math::operator/ ( Integer a, Rational const & b )

inline

Division operation between Rational and integral type.

Template Parameters

Integer

Integral type

Parameters

a	Left operand
b	Right operand

Returns

Result of division

ReferenceMomentFittingMatrix() [1/3]

template<polynomial::CBasis TBasis, class TDerivedXg>

Matrix< TBasis::kSize, Eigen::Dynamic > pbat::math::ReferenceMomentFittingMatrix	(	TBasis const &	Pb,
		Eigen::MatrixBase< TDerivedXg > const &	Xg )

Compute the moment fitting matrix in the reference simplex space.

Template Parameters

TBasis	Polynomial basis type
TDerivedXg	Eigen expression for quadrature points

Parameters

Pb	Polynomial basis
Xg	Quadrature points

Returns

Moment fitting matrix

ReferenceMomentFittingMatrix() [2/3]

template<polynomial::CBasis TBasis, CFixedPointPolynomialQuadratureRule TQuad>

Matrix< TBasis::kSize, TQuad::kPoints > pbat::math::ReferenceMomentFittingMatrix	(	TBasis const &	Pb,
		TQuad const &	Q )

Compute the moment fitting matrix in the reference simplex space.

Template Parameters

TBasis	Polynomial basis type
TQuad	Quadrature rule type

Parameters

Pb	Polynomial basis
Q	Quadrature rule

Returns

Moment fitting matrix

ReferenceMomentFittingMatrix() [3/3]

template<polynomial::CBasis TBasis, CPolynomialQuadratureRule TQuad>

Matrix< TBasis::kSize, Eigen::Dynamic > pbat::math::ReferenceMomentFittingMatrix	(	TBasis const &	Pb,
		TQuad const &	Q )

Compute the moment fitting matrix in the reference simplex space.

Template Parameters

TBasis	Polynomial basis type
TQuad	Quadrature rule type

Parameters

Pb	Polynomial basis
Q	Quadrature rule

Returns

Moment fitting matrix

ReferenceMomentFittingSystems()

template<int Order, class TDerivedS1, class TDerivedX1, class TDerivedS2, class TDerivedX2, class TDerivedW2>

std::tuple< MatrixX, MatrixX, IndexVectorX > pbat::math::ReferenceMomentFittingSystems	(	Eigen::DenseBase< TDerivedS1 > const &	S1,
		Eigen::MatrixBase< TDerivedX1 > const &	X1,
		Eigen::DenseBase< TDerivedS2 > const &	S2,
		Eigen::MatrixBase< TDerivedX2 > const &	X2,
		Eigen::DenseBase< TDerivedW2 > const &	w2,
		Index	nSimplices = Index(-1) )

Computes reference moment fitting systems for all simplices \( S_1 \), given an existing quadrature rule \( (X_g^2, w_g^2) \) defined on a domain of simplices. The quadrature points of the new rule \( X_g^1 \) are given and fixed.

Template Parameters

Order	Order of the quadrature rules
TDerivedS1	Eigen vector expression for simplex indices
TDerivedX1	Eigen matrix expression for quadrature points
TDerivedS2	Eigen vector expression for simplex indices
TDerivedX2	Eigen matrix expression for existing quadrature points
TDerivedW2	Eigen vector expression for existing quadrature weights

Parameters

S1	\( |S_1| \times 1 \) index array giving the simplex containing the corresponding quadrature point in columns of X1, i.e. S1(g) is the simplex containing X1.col(g)
X1	\( d \times n_1 \) array \( X_g^1 \) of quadrature point positions defined in reference space
S2	\( |S_2| \times 1 \) index array giving the simplex containing the corresponding quadrature point in columns of X2, i.e. S2(g) is the simplex containing X2.col(g)
X2	\( d \times n_2 \) array of quadrature point positions defined in reference simplex space
w2	\( n_2 \times 1 \) array \( w_g^2 \) of existing quadrature weights
nSimplices	Number of simplices in the domain. If nSimplices < 1, the number of simplices is inferred from S1 and S2.

Returns

(P, B, prefix) where P is the moment fitting matrix, B is the target integrated polynomials, and prefix is the prefix into columns of P for each simplex, i.e. the block P[:,prefix[s]:prefix[s+1]] = B[:,s] represents the moment fitting system for simplex s.

TransferQuadrature() [1/2]

template<auto Order, class TDerivedS1, class TDerivedXi1, class TDerivedS2, class TDerivedXi2, class TDerivedWg2>

std::pair< VectorX, VectorX > pbat::math::TransferQuadrature	(	Eigen::DenseBase< TDerivedS1 > const &	S1,
		Eigen::MatrixBase< TDerivedXi1 > const &	Xi1,
		Eigen::DenseBase< TDerivedS2 > const &	S2,
		Eigen::MatrixBase< TDerivedXi2 > const &	Xi2,
		Eigen::DenseBase< TDerivedWg2 > const &	wi2,
		Index	nSimplices = -1,
		bool	bEvaluateError = false,
		Index	maxIterations = 10,
		Scalar	precision = std::numeric_limits<Scalar>::epsilon() )

Obtain weights \( w_g^1 \) by transferring an existing quadrature rule \( X_g^2, w_g^2 \) defined on a domain of simplices onto a new quadrature rule \( X_g^1, w_g^1 \) defined on the same domain.

Template Parameters

Order	Order of the quadrature rules
TDerivedS1	Eigen vector expression for simplex indices
TDerivedXi1	Eigen matrix expression for quadrature points
TDerivedS2	Eigen vector expression for simplex indices
TDerivedXi2	Eigen matrix expression for existing quadrature points
TDerivedWg2	Eigen vector expression for existing quadrature weights

Parameters

S1	\( |S_1| \times 1 \) index array giving the simplex containing the corresponding quadrature point in columns of Xi1, i.e. S1(g) is the simplex containing Xi1.col(g)
Xi1	\( d \times n_1 \) array \( \xi_g^1 \) of quadrature point positions defined in simplex space
S2	\( |S_2| \times 1 \) index array giving the simplex containing the corresponding quadrature point in columns of Xi2, i.e. S2(g) is the simplex containing Xi2.col(g)
Xi2	\( d \times n_2 \) array of quadrature point positions defined in reference simplex space
wi2	\( n_2 \times 1 \) array \( w_g^2 \) of existing quadrature weights
nSimplices	Number of simplices in the domain. If nSimplices < 1, the number of simplices is inferred from S1 and S2.
bEvaluateError	Whether to compute the integration error on the new quadrature rule
maxIterations	Maximum number of non-negative least-squares active set solver
precision	Convergence threshold

Returns

(w, e) where w are the quadrature weights associated with points Xi1, and e is the integration error in each simplex (zeros if not bEvaluateError).

TransferQuadrature() [2/2]

template<polynomial::CBasis Polynomial, class TDerivedXg1, class TDerivedXg2, class TDerivedWg2>

Vector< TDerivedXg1::ColsAtCompileTime > pbat::math::TransferQuadrature	(	Polynomial const &	P,
		Eigen::MatrixBase< TDerivedXg1 > const &	Xg1,
		Eigen::MatrixBase< TDerivedXg2 > const &	Xg2,
		Eigen::DenseBase< TDerivedWg2 > const &	wg2,
		Index	maxIterations = 10,
		Scalar	precision = std::numeric_limits<Scalar>::epsilon() )

Obtain weights \( w_g^1 \) by transferring an existing quadrature rule \( X_g^2, w_g^2 \) defined on a simplex onto a new quadrature rule \( X_g^1, w_g^1 \) defined on the same simplex.

Template Parameters

TDerivedXg1	Eigen matrix expression for quadrature points
TDerivedXg2	Eigen matrix expression for existing quadrature points
TDerivedWg2	Eigen vector expression for existing quadrature weights
Polynomial	Polynomial basis type

Parameters

P	Polynomial basis
Xg1	\( d \times n_1 \) array \( X_g^1 \) of quadrature point positions defined in reference space
Xg2	\( d \times n_2 \) array \( X_g^2 \) of quadrature point positions defined in reference space
wg2	\( n_2 \times 1 \) array \( w_g^2 \) of existing quadrature weights
maxIterations	Maximum number of non-negative least-squares active set solver
precision	Convergence threshold

Returns

\( n_1 \times 1 \) array of quadrature weights