Source code for fe_utils.quadrature

from numpy.polynomial.legendre import leggauss
from .reference_elements import ReferenceInterval, ReferenceTriangle
import numpy as np


[docs]class QuadratureRule(object): def __init__(self, cell, degree, points, weights): """A quadrature rule implementing integration of the reference cell provided. :param cell: the :class:`~.ReferenceCell` over which this quadrature rule is defined. :param degree: the :ref:`degree of precision <degree-of-precision>` of this quadrature rule. :points: a list of the position vectors of the quadrature points. :weights: the corresponding vector of quadrature weights. """ #: The :class:`~.ReferenceCell` over which this quadrature #: rule is defined. self.cell = cell #: Two dimensional array, the rows of which form the position #: vectors of the quadrature points. self.points = np.array(points, dtype=np.double) #: The corresponding array of quadrature weights. self.weights = np.array(weights, dtype=np.double) #: The degree of precision of the quadrature rule. self.degree = degree if self.cell.dim != self.points.shape[1]: raise ValueError( "Dimension mismatch between reference cell " "and quadrature points" ) if self.points.shape[0] != len(self.weights): raise ValueError( "Number of quadrature points and quadrature weights must match" )
[docs] def integrate(self, function): """Integrate the function provided using this quadrature rule. :param function: A Python function taking a position vector as its single argument and returning a scalar value. The implementation of this method is left as an :ref:`exercise <ex-integrate>`. """ raise NotImplementedError
[docs]def gauss_quadrature(cell, degree): """Return a Gauss-Legendre :class:`QuadratureRule`. :param cell: the :class:`~.ReferenceCell` over which this quadrature rule is defined. :param degree: the :ref:`degree of precision <degree-of-precision>` of this quadrature rule. """ if cell is ReferenceInterval: # We can obtain the 1D gauss-legendre rule from numpy # and change coordinates. # Gauss-legendre quadrature has degree = 2 * npoints - 1 # The extra + 1 deals with truncation. npoints = int((degree + 1 + 1) / 2) points, weights = leggauss(npoints) # Numpy uses the [-1, 1] interval, we need to change to [0, 1] points = (points + 1.) / 2. # Rescale the weights to sum to 1 instead of 2. weights = weights / 2. # We expect points to be an n x dim array. points.shape = [points.shape[0], 1] elif cell is ReferenceTriangle: # The 2D rule is obtained using the 1D rule and the Duffy Transform. p1 = gauss_quadrature(ReferenceInterval, degree + 1) q1 = gauss_quadrature(ReferenceInterval, degree) points = np.array([(p[0], q[0] * (1 - p[0])) for p in p1.points for q in q1.points]) weights = np.array([p * q * (1 - x[0]) for p, x in zip(p1.weights, p1.points) for q in q1.weights]) else: raise ValueError("Unknown reference cell") return QuadratureRule(cell, degree, points, weights)