# Path: finat/point_set.py
# class PointSet(AbstractPointSet):
#     """A basic point set with no internal structure representing a vector of
#     points."""
# 
#     def __init__(self, points):
#         """Build a PointSet from a vector of points
# 
#         :arg points: A vector of N points of shape (N, D) where D is the
#             dimension of each point."""
#         points = numpy.asarray(points)
#         assert len(points.shape) == 2
#         self.points = points
# 
#     @cached_property
#     def points(self):
#         pass  # set at initialisation
# 
#     @cached_property
#     def indices(self):
#         return (gem.Index(extent=len(self.points)),)
# 
#     @cached_property
#     def expression(self):
#         return gem.partial_indexed(gem.Literal(self.points), self.indices)
# 
#     def almost_equal(self, other, tolerance=1e-12):
#         """Approximate numerical equality of point sets"""
#         return type(self) == type(other) and \
#             self.points.shape == other.points.shape and \
#             numpy.allclose(self.points, other.points, rtol=0, atol=tolerance)
# 
# class GaussLegendrePointSet(PointSet):
#     """Gauss-Legendre quadrature points on the interval.
# 
#     This facilitates implementing discontinuous spectral elements.
#     """
#     def __init__(self, points):
#         super(GaussLegendrePointSet, self).__init__(points)
#         assert self.points.shape[1] == 1
# 
# class TensorPointSet(AbstractPointSet):
# 
#     def __init__(self, factors):
#         self.factors = tuple(factors)
# 
#     @cached_property
#     def points(self):
#         return numpy.array([list(chain(*pt_tuple))
#                             for pt_tuple in product(*[ps.points
#                                                       for ps in self.factors])])
# 
#     @cached_property
#     def indices(self):
#         return tuple(chain(*[ps.indices for ps in self.factors]))
# 
#     @cached_property
#     def expression(self):
#         result = []
#         for point_set in self.factors:
#             for i in range(point_set.dimension):
#                 result.append(gem.Indexed(point_set.expression, (i,)))
#         return gem.ListTensor(result)
# 
#     def almost_equal(self, other, tolerance=1e-12):
#         """Approximate numerical equality of point sets"""
#         return type(self) == type(other) and \
#             len(self.factors) == len(other.factors) and \
#             all(s.almost_equal(o, tolerance=tolerance)
#                 for s, o in zip(self.factors, other.factors))


from abc import ABCMeta, abstractproperty
from functools import reduce
from gem.utils import cached_property
from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT
from FIAT.quadrature import GaussLegendreQuadratureLineRule
from FIAT.quadrature_schemes import create_quadrature as fiat_scheme
from finat.point_set import PointSet, GaussLegendrePointSet, TensorPointSet
import numpy
import gem
        raise ValueError("Need positive degree, not %d" % degree)

    if ref_el.get_shape() == LINE:
        # FIAT uses Gauss-Legendre line quadature, however, since we
        # symbolically label it as such, we wish not to risk attaching
        # the wrong label in case FIAT changes.  So we explicitly ask
        # for Gauss-Legendre line quadature.
        num_points = (degree + 1 + 1) // 2  # exact integration
        fiat_rule = GaussLegendreQuadratureLineRule(ref_el, num_points)
        point_set = GaussLegendrePointSet(fiat_rule.get_points())
        return QuadratureRule(point_set, fiat_rule.get_weights())

    fiat_rule = fiat_scheme(ref_el, degree, scheme)
    return QuadratureRule(PointSet(fiat_rule.get_points()), fiat_rule.get_weights())


class AbstractQuadratureRule(metaclass=ABCMeta):
    """Abstract class representing a quadrature rule as point set and a
    corresponding set of weights."""

    @abstractproperty
    def point_set(self):
        """Point set object representing the quadrature points."""

    @abstractproperty
    def weight_expression(self):
        """GEM expression describing the weights, with the same free indices
        as the point set."""



# Path: finat/point_set.py
# class PointSet(AbstractPointSet):
#     """A basic point set with no internal structure representing a vector of
#     points."""
# 
#     def __init__(self, points):
#         """Build a PointSet from a vector of points
# 
#         :arg points: A vector of N points of shape (N, D) where D is the
#             dimension of each point."""
#         points = numpy.asarray(points)
#         assert len(points.shape) == 2
#         self.points = points
# 
#     @cached_property
#     def points(self):
#         pass  # set at initialisation
# 
#     @cached_property
#     def indices(self):
#         return (gem.Index(extent=len(self.points)),)
# 
#     @cached_property
#     def expression(self):
#         return gem.partial_indexed(gem.Literal(self.points), self.indices)
# 
#     def almost_equal(self, other, tolerance=1e-12):
#         """Approximate numerical equality of point sets"""
#         return type(self) == type(other) and \
#             self.points.shape == other.points.shape and \
#             numpy.allclose(self.points, other.points, rtol=0, atol=tolerance)
# 
# class GaussLegendrePointSet(PointSet):
#     """Gauss-Legendre quadrature points on the interval.
# 
#     This facilitates implementing discontinuous spectral elements.
#     """
#     def __init__(self, points):
#         super(GaussLegendrePointSet, self).__init__(points)
#         assert self.points.shape[1] == 1
# 
# class TensorPointSet(AbstractPointSet):
# 
#     def __init__(self, factors):
#         self.factors = tuple(factors)
# 
#     @cached_property
#     def points(self):
#         return numpy.array([list(chain(*pt_tuple))
#                             for pt_tuple in product(*[ps.points
#                                                       for ps in self.factors])])
# 
#     @cached_property
#     def indices(self):
#         return tuple(chain(*[ps.indices for ps in self.factors]))
# 
#     @cached_property
#     def expression(self):
#         result = []
#         for point_set in self.factors:
#             for i in range(point_set.dimension):
#                 result.append(gem.Indexed(point_set.expression, (i,)))
#         return gem.ListTensor(result)
# 
#     def almost_equal(self, other, tolerance=1e-12):
#         """Approximate numerical equality of point sets"""
#         return type(self) == type(other) and \
#             len(self.factors) == len(other.factors) and \
#             all(s.almost_equal(o, tolerance=tolerance)
#                 for s, o in zip(self.factors, other.factors))


# Path: finat/quadrature.py
from abc import ABCMeta, abstractproperty
from functools import reduce
from gem.utils import cached_property
from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT
from FIAT.quadrature import GaussLegendreQuadratureLineRule
from FIAT.quadrature_schemes import create_quadrature as fiat_scheme
from finat.point_set import PointSet, GaussLegendrePointSet, TensorPointSet
import numpy
import gem

        raise ValueError("Need positive degree, not %d" % degree)

    if ref_el.get_shape() == LINE:
        # FIAT uses Gauss-Legendre line quadature, however, since we
        # symbolically label it as such, we wish not to risk attaching
        # the wrong label in case FIAT changes.  So we explicitly ask
        # for Gauss-Legendre line quadature.
        num_points = (degree + 1 + 1) // 2  # exact integration
        fiat_rule = GaussLegendreQuadratureLineRule(ref_el, num_points)
        point_set = GaussLegendrePointSet(fiat_rule.get_points())
        return QuadratureRule(point_set, fiat_rule.get_weights())

    fiat_rule = fiat_scheme(ref_el, degree, scheme)
    return QuadratureRule(PointSet(fiat_rule.get_points()), fiat_rule.get_weights())


class AbstractQuadratureRule(metaclass=ABCMeta):
    """Abstract class representing a quadrature rule as point set and a
    corresponding set of weights."""

    @abstractproperty
    def point_set(self):
        """Point set object representing the quadrature points."""

    @abstractproperty
    def weight_expression(self):
        """GEM expression describing the weights, with the same free indices
        as the point set."""



class QuadratureRule(AbstractQuadratureRule):