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

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>`
: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:
raise ValueError(
"Dimension mismatch between reference cell "
)
if self.points.shape != len(self.weights):
raise ValueError(
)

[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

:param cell: the :class:`~.ReferenceCell` over which this quadrature
rule is defined.
:param degree: the :ref:`degree of precision <degree-of-precision>`
"""

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, 1]

elif cell is ReferenceTriangle:
# The 2D rule is obtained using the 1D rule and the Duffy Transform.

p1 = gauss_quadrature(ReferenceInterval, degree + 1)

points = np.array([(p, q * (1 - p))
for p in p1.points
for q in q1.points])

weights = np.array([p * q * (1 - x)
for p, x in zip(p1.weights, p1.points)
for q in q1.weights])

else:
raise ValueError("Unknown reference cell")