from abc import ABC, abstractmethod
from typing import Optional, Tuple
import numpy as np
import numpy.polynomial.legendre as leg
[docs]
class QuadratureRule(ABC):
"""Abstract base class for a quadrature rule."""
[docs]
@abstractmethod
def get_points_weights(self) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover
"""Retrieve the quadrature points and weights."""
pass
[docs]
@staticmethod
@abstractmethod
def compute_points_weights(self) -> Tuple[np.ndarray, np.ndarray]: # pragma: no cover
"""Compute the quadrature points and weights."""
pass
@property
@abstractmethod
def dimension(self) -> int: # pragma: no cover
"""The number of dimensions of the quadrature rule."""
pass
class QuadratureRule2D(QuadratureRule):
"""Abstract base class for a 2D quadrature rule."""
@property
def dimension(self) -> int:
"""The number of dimensions of the quadrature rule."""
return 2
[docs]
class GaussLegendre(QuadratureRule):
"""
Implements Gauss-Legendre quadrature on [-1, 1]^(dimension).
Example:
>>> from ttfemesh.quadrature import GaussLegendre
>>> order = 3
>>> dim = 3
>>> qrule = GaussLegendre(order, dim)
>>> points, weights = qrule.get_points_weights()
>>> print(points)
>>> print(weights)
"""
[docs]
def __init__(self, order: int = 1, dimension: int = 1):
"""
Initialize the Gauss-Legendre quadrature rule.
Args:
order (int): The order of the quadrature rule. Must be greater than 0 (default is 1).
dimension (int): The number of dimensions (default is 1).
"""
self.order = order
self._dimension = dimension
self.points: Optional[np.ndarray] = None
self.weights: Optional[np.ndarray] = None
@property
def dimension(self) -> int:
"""The number of dimensions of the quadrature rule."""
return self._dimension
[docs]
def get_points_weights(self) -> Tuple[np.ndarray, np.ndarray]:
"""
Get the quadrature points and weights.
Returns:
Tuple[np.ndarray, np.ndarray]: Quadrature points and weights.
"""
if self.points is None or self.weights is None:
self.points, self.weights = self.compute_points_weights(self.order, self.dimension)
return self.points, self.weights
[docs]
@staticmethod
def compute_points_weights(order: int, dimension: int = 1) -> Tuple[np.ndarray, np.ndarray]:
"""
Compute the Gauss-Legendre quadrature points and weights.
Args:
order (int): The order of the quadrature rule. Must be greater than 0.
dimension (int): The number of dimensions (default is 1).
Returns:
Tuple[np.ndarray, np.ndarray]: Quadrature points and weights.
"""
points_1d, weights_1d = leg.leggauss(order)
if dimension == 1:
return points_1d[:, None], weights_1d
# Tensor product for multi-dimensional quadrature
grid = np.meshgrid(*[points_1d] * dimension, indexing="ij")
points = np.stack([g.flatten() for g in grid], axis=-1)
weights = np.prod(np.meshgrid(*[weights_1d] * dimension, indexing="ij"), axis=0)
return points, weights.flatten()
[docs]
def __repr__(self) -> str:
"""String representation of the Gauss-Legendre quadrature rule."""
return f"Gauss-Legendre Quadrature Rule (order={self.order}, dimension={self.dimension})"
[docs]
class GaussLegendre2D(GaussLegendre):
"""
Implements Gauss-Legendre quadrature on [-1, 1]^2.
Example:
>>> from ttfemesh.quadrature import GaussLegendre2D
>>> order = 3
>>> qrule2D = GaussLegendre2D(order)
>>> points2D, weights2D = qrule2D.get_points_weights()
>>> print(points2D)
>>> print(weights2D)
"""
[docs]
def __init__(self, order: int = 1):
"""
Initialize the 2D Gauss-Legendre quadrature rule.
Args:
order (int): The order of the quadrature rule. Must be greater than 0 (default is 1).
"""
super().__init__(order, dimension=2)
[docs]
def __repr__(self) -> str:
"""String representation of the 2D Gauss-Legendre quadrature rule."""
return f"2D Gauss-Legendre Quadrature Rule (order={self.order})"