# Path: finat/quadrature.py
# def make_quadrature(ref_el, degree, scheme="default"):
#     """
#     Generate quadrature rule for given reference element
#     that will integrate an polynomial of order 'degree' exactly.
# 
#     For low-degree (<=6) polynomials on triangles and tetrahedra, this
#     uses hard-coded rules, otherwise it falls back to a collapsed
#     Gauss scheme on simplices.  On tensor-product cells, it is a
#     tensor-product quadrature rule of the subcells.
# 
#     :arg ref_el: The FIAT cell to create the quadrature for.
#     :arg degree: The degree of polynomial that the rule should
#         integrate exactly.
#     """
#     if ref_el.get_shape() == TENSORPRODUCT:
#         try:
#             degree = tuple(degree)
#         except TypeError:
#             degree = (degree,) * len(ref_el.cells)
# 
#         assert len(ref_el.cells) == len(degree)
#         quad_rules = [make_quadrature(c, d, scheme)
#                       for c, d in zip(ref_el.cells, degree)]
#         return TensorProductQuadratureRule(quad_rules)
# 
#     if ref_el.get_shape() == QUADRILATERAL:
#         return make_quadrature(ref_el.product, degree, scheme)
# 
#     if degree < 0:
#         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())


from abc import ABCMeta, abstractproperty, abstractmethod
from itertools import chain
from gem.interpreter import evaluate
from gem.optimise import delta_elimination, sum_factorise, traverse_product
from gem.utils import cached_property
from finat.quadrature import make_quadrature
import numpy
import gem
        Note that there are two entities on dimension ``0`` (vertices),
        each of which has only one possible orientation, while there is
        a single entity on dimension ``1`` (interval), which has two
        possible orientations representing non-reflected and reflected
        intervals.
        '''
        raise NotImplementedError(f"entity_permutations not yet implemented for {type(self)}")

    @cached_property
    def _entity_closure_dofs(self):
        # Compute the nodes on the closure of each sub_entity.
        entity_dofs = self.entity_dofs()
        return {dim: {e: sorted(chain(*[entity_dofs[d][se]
                                        for d, se in sub_entities]))
                      for e, sub_entities in entities.items()}
                for dim, entities in self.cell.sub_entities.items()}

    def entity_closure_dofs(self):
        '''Return the map of topological entities to degrees of
        freedom on the closure of those entities for the finite
        element.'''
        return self._entity_closure_dofs

    @cached_property
    def _entity_support_dofs(self):
        esd = {}
        for entity_dim in self.cell.sub_entities.keys():
            beta = self.get_indices()
            zeta = self.get_value_indices()


# Path: finat/quadrature.py
# def make_quadrature(ref_el, degree, scheme="default"):
#     """
#     Generate quadrature rule for given reference element
#     that will integrate an polynomial of order 'degree' exactly.
# 
#     For low-degree (<=6) polynomials on triangles and tetrahedra, this
#     uses hard-coded rules, otherwise it falls back to a collapsed
#     Gauss scheme on simplices.  On tensor-product cells, it is a
#     tensor-product quadrature rule of the subcells.
# 
#     :arg ref_el: The FIAT cell to create the quadrature for.
#     :arg degree: The degree of polynomial that the rule should
#         integrate exactly.
#     """
#     if ref_el.get_shape() == TENSORPRODUCT:
#         try:
#             degree = tuple(degree)
#         except TypeError:
#             degree = (degree,) * len(ref_el.cells)
# 
#         assert len(ref_el.cells) == len(degree)
#         quad_rules = [make_quadrature(c, d, scheme)
#                       for c, d in zip(ref_el.cells, degree)]
#         return TensorProductQuadratureRule(quad_rules)
# 
#     if ref_el.get_shape() == QUADRILATERAL:
#         return make_quadrature(ref_el.product, degree, scheme)
# 
#     if degree < 0:
#         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())


# Path: finat/finiteelementbase.py
from abc import ABCMeta, abstractproperty, abstractmethod
from itertools import chain
from gem.interpreter import evaluate
from gem.optimise import delta_elimination, sum_factorise, traverse_product
from gem.utils import cached_property
from finat.quadrature import make_quadrature
import numpy
import gem

        Note that there are two entities on dimension ``0`` (vertices),
        each of which has only one possible orientation, while there is
        a single entity on dimension ``1`` (interval), which has two
        possible orientations representing non-reflected and reflected
        intervals.
        '''
        raise NotImplementedError(f"entity_permutations not yet implemented for {type(self)}")

    @cached_property
    def _entity_closure_dofs(self):
        # Compute the nodes on the closure of each sub_entity.
        entity_dofs = self.entity_dofs()
        return {dim: {e: sorted(chain(*[entity_dofs[d][se]
                                        for d, se in sub_entities]))
                      for e, sub_entities in entities.items()}
                for dim, entities in self.cell.sub_entities.items()}

    def entity_closure_dofs(self):
        '''Return the map of topological entities to degrees of
        freedom on the closure of those entities for the finite
        element.'''
        return self._entity_closure_dofs

    @cached_property
    def _entity_support_dofs(self):
        esd = {}
        for entity_dim in self.cell.sub_entities.keys():
            beta = self.get_indices()
            zeta = self.get_value_indices()


            entity_cell = self.cell.construct_subelement(entity_dim)