def DiscontinuousLagrange(ref_el, degree):
    if degree == 0:
        return P0.P0(ref_el)
    else:
        return HigherOrderDiscontinuousLagrange(ref_el, degree) class DualSet(object):
    def __init__(self, nodes, ref_el, entity_ids):
        self.nodes = nodes
        self.ref_el = ref_el
        self.entity_ids = entity_ids

        # Compute the nodes on the closure of each sub_entity.
        self.entity_closure_ids = {}
        for dim, entities in ref_el.sub_entities.items():
            self.entity_closure_ids[dim] = {}

            for e, sub_entities in entities.items():
                ids = []

                for d, se in sub_entities:
                    ids += self.entity_ids[d][se]
                ids.sort()
                self.entity_closure_ids[d][e] = ids

    def get_nodes(self):
        return self.nodes

    def get_entity_closure_ids(self):
        return self.entity_closure_ids

    def get_entity_ids(self):
        return self.entity_ids

    def get_reference_element(self):
        return self.ref_el

    def to_riesz(self, poly_set):
        r"""This method gives the action of the entire dual set
        on each member of the expansion set underlying poly_set.
        Then, applying the linear functionals of the dual set to an
        arbitrary polynomial in poly_set is accomplished by (generalized)
        matrix multiplication.

        For scalar-valued spaces, this produces a matrix
        :\math:`R_{i, j}` such that
        :\math:`\ell_i(f) = \sum_{j} a_j \ell_i(\phi_j)`
        for :\math:`f=\sum_{j} a_j \phi_j`.

        More generally, it will have shape concatenating
        the number of functionals in the dual set, the value shape
        of functions it takes, and the number of members of the
        expansion set.
        """

        # This rather technical code queries the low-level information
        # in pt_dict and deriv_dict
        # for each functional to find out where it evaluates its
        # inputs and/or their derivatives.  Then, it tabulates the
        # expansion set one time for all the function values and
        # another for all of the derivatives.  This circumvents
        # needing to call the to_riesz method of each functional and
        # also limits the number of different calls to tabulate.

        tshape = self.nodes[0].target_shape
        num_nodes = len(self.nodes)
        es = poly_set.get_expansion_set()
        ed = poly_set.get_embedded_degree()
        num_exp = es.get_num_members(poly_set.get_embedded_degree())

        riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp])

        self.mat = numpy.zeros(riesz_shape, "d")

        # Dictionaries mapping pts to which functionals they come from
        pts_to_ells = collections.OrderedDict()
        dpts_to_ells = collections.OrderedDict()

        for i, ell in enumerate(self.nodes):
            for pt in ell.pt_dict:
                if pt in pts_to_ells:
                    pts_to_ells[pt].append(i)
                else:
                    pts_to_ells[pt] = [i]

            for pt in ell.deriv_dict:
                if pt in dpts_to_ells:
                    dpts_to_ells[pt].append(i)
                else:
                    dpts_to_ells[pt] = [i]

        # Now tabulate the function values
        pts = list(pts_to_ells.keys())
        expansion_values = es.tabulate(ed, pts)

        for j, pt in enumerate(pts):
            which_ells = pts_to_ells[pt]

            for k in which_ells:
                pt_dict = self.nodes[k].pt_dict
                wc_list = pt_dict[pt]

                for i in range(num_exp):
                    for (w, c) in wc_list:
                        self.mat[k][c][i] += w*expansion_values[i, j]

        # Tabulate the derivative values that are needed
        max_deriv_order = max([ell.max_deriv_order for ell in self.nodes])
        if max_deriv_order > 0:
            dpts = list(dpts_to_ells.keys())
            # It's easiest/most efficient to get derivatives of the
            # expansion set through the polynomial set interface.
            # This is creating a short-lived set to do just this.
            expansion = polynomial_set.ONPolynomialSet(self.ref_el, ed)
            dexpansion_values = expansion.tabulate(dpts, max_deriv_order)

            for j, pt in enumerate(dpts):
                which_ells = dpts_to_ells[pt]

                for k in which_ells:
                    dpt_dict = self.nodes[k].deriv_dict
                    wac_list = dpt_dict[pt]

                    for i in range(num_exp):
                        for (w, alpha, c) in wac_list:
                            self.mat[k][c][i] += w*dexpansion_values[alpha][i, j]

        return self.mat class FiniteElement(object):
    """Class implementing a basic abstraction template for general
    finite element families. Finite elements which inherit from
    this class are non-nodal unless they are CiarletElement subclasses.
    """

    def __init__(self, ref_el, dual, order, formdegree=None, mapping="affine"):
        # Relevant attributes that do not necessarily depend on a PolynomialSet object:
        # The order (degree) of the polynomial basis
        self.order = order
        self.formdegree = formdegree

        # The reference element and the appropriate dual
        self.ref_el = ref_el
        self.dual = dual

        # The appropriate mapping for the finite element space
        self._mapping = mapping

    def get_reference_element(self):
        """Return the reference element for the finite element."""
        return self.ref_el

    def get_dual_set(self):
        """Return the dual for the finite element."""
        return self.dual

    def get_order(self):
        """Return the order of the element (may be different from the degree)."""
        return self.order

    def dual_basis(self):
        """Return the dual basis (list of functionals) for the finite
        element."""
        return self.dual.get_nodes()

    def entity_dofs(self):
        """Return the map of topological entities to degrees of
        freedom for the finite element."""
        return self.dual.get_entity_ids()

    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.dual.get_entity_closure_ids()

    def get_formdegree(self):
        """Return the degree of the associated form (FEEC)"""
        return self.formdegree

    def mapping(self):
        """Return a list of appropriate mappings from the reference
        element to a physical element for each basis function of the
        finite element."""
        return [self._mapping] * self.space_dimension()

    def num_sub_elements(self):
        """Return the number of sub-elements."""
        return 1

    def space_dimension(self):
        """Return the dimension of the finite element space."""
        return len(self.dual_basis())

    def tabulate(self, order, points, entity=None):
        """Return tabulated values of derivatives up to given order of
        basis functions at given points.

        :arg order: The maximum order of derivative.
        :arg points: An iterable of points.
        :arg entity: Optional (dimension, entity number) pair
                     indicating which topological entity of the
                     reference element to tabulate on.  If ``None``,
                     default cell-wise tabulation is performed.
        """
        raise NotImplementedError("Must be specified in the element subclass of FiniteElement.")

    @staticmethod
    def is_nodal():
        """True if primal and dual bases are orthogonal. If false,
        dual basis is not implemented or is undefined.

        Subclasses may not necessarily be nodal, unless it is a CiarletElement.
        """
        return False class PointEvaluation(Functional):
    """Class representing point evaluation of scalar functions at a
    particular point x."""

    def __init__(self, ref_el, x):
        pt_dict = {x: [(1.0, tuple())]}
        Functional.__init__(self, ref_el, tuple(), pt_dict, {}, "PointEval")

    def __call__(self, fn):
        """Evaluate the functional on the function fn."""
        return fn(tuple(self.pt_dict.keys())[0])

    def tostr(self):
        x = list(map(str, list(self.pt_dict.keys())[0]))
        return "u(%s)" % (','.join(x),) def mis(m, n):
    """Returns all m-tuples of nonnegative integers that sum up to n."""
    if m == 1:
        return [(n,)]
    elif n == 0:
        return [tuple([0] * m)]
    else:
        return [tuple([n - i] + list(foo))
                for i in range(n + 1)
                for foo in mis(m - 1, i)] def ufc_simplex(spatial_dim):
    """Factory function that maps spatial dimension to an instance of
    the UFC reference simplex of that dimension."""
    if spatial_dim == 0:
        return Point()
    elif spatial_dim == 1:
        return UFCInterval()
    elif spatial_dim == 2:
        return UFCTriangle()
    elif spatial_dim == 3:
        return UFCTetrahedron()
    else:
        raise RuntimeError("Can't create UFC simplex of dimension %s." % str(spatial_dim)) POINT = 0 LINE = 1 QUADRILATERAL = 11 TRIANGLE = 2 TETRAHEDRON = 3 TENSORPRODUCT = 99 class TensorProductElement(FiniteElement):
    """Class implementing a finite element that is the tensor product
    of two existing finite elements."""

    def __init__(self, A, B):
        # set up simple things
        order = min(A.get_order(), B.get_order())
        if A.get_formdegree() is None or B.get_formdegree() is None:
            formdegree = None
        else:
            formdegree = A.get_formdegree() + B.get_formdegree()

        # set up reference element
        ref_el = TensorProductCell(A.get_reference_element(),
                                   B.get_reference_element())

        if A.mapping()[0] != "affine" and B.mapping()[0] == "affine":
            mapping = A.mapping()[0]
        elif B.mapping()[0] != "affine" and A.mapping()[0] == "affine":
            mapping = B.mapping()[0]
        elif A.mapping()[0] == "affine" and B.mapping()[0] == "affine":
            mapping = "affine"
        else:
            raise ValueError("check tensor product mappings - at least one must be affine")

        # set up entity_ids
        Adofs = A.entity_dofs()
        Bdofs = B.entity_dofs()
        Bsdim = B.space_dimension()
        entity_ids = {}

        for curAdim in Adofs:
            for curBdim in Bdofs:
                entity_ids[(curAdim, curBdim)] = {}
                dim_cur = 0
                for entityA in Adofs[curAdim]:
                    for entityB in Bdofs[curBdim]:
                        entity_ids[(curAdim, curBdim)][dim_cur] = \
                            [x*Bsdim + y for x in Adofs[curAdim][entityA]
                                for y in Bdofs[curBdim][entityB]]
                        dim_cur += 1

        # set up dual basis
        Anodes = A.dual_basis()
        Bnodes = B.dual_basis()

        # build the dual set by inspecting the current dual
        # sets item by item.
        # Currently supported cases:
        # PointEval x PointEval = PointEval [scalar x scalar = scalar]
        # PointScaledNormalEval x PointEval = PointScaledNormalEval [vector x scalar = vector]
        # ComponentPointEvaluation x PointEval [vector x scalar = vector]
        nodes = []
        for Anode in Anodes:
            if isinstance(Anode, functional.PointEvaluation):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: PointEval x PointEval
                        # the PointEval functional just requires the
                        # coordinates. these are currently stored as
                        # the key of a one-item dictionary. we retrieve
                        # these by calling get_point_dict(), and
                        # use the concatenation to make a new PointEval
                        nodes.append(functional.PointEvaluation(ref_el, _first_point(Anode) + _first_point(Bnode)))
                    elif isinstance(Bnode, functional.IntegralMoment):
                        # dummy functional for product with integral moments
                        nodes.append(functional.Functional(None, None, None,
                                                           {}, "Undefined"))
                    elif isinstance(Bnode, functional.PointDerivative):
                        # dummy functional for product with point derivative
                        nodes.append(functional.Functional(None, None, None,
                                                           {}, "Undefined"))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.PointScaledNormalEvaluation):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: PointScaledNormalEval x PointEval
                        # this could be wrong if the second shape
                        # has spatial dimension >1, since we are not
                        # explicitly scaling by facet size
                        if len(_first_point(Bnode)) > 1:
                            # TODO: support this case one day
                            raise NotImplementedError("PointScaledNormalEval x PointEval is not yet supported if the second shape has dimension > 1")
                        # We cannot make a new functional.PSNEval in
                        # the natural way, since it tries to compute
                        # the normal vector by itself.
                        # Instead, we create things manually, and
                        # call Functional() with these arguments
                        sd = ref_el.get_spatial_dimension()
                        # The pt_dict is a one-item dictionary containing
                        # the details of the functional.
                        # The key is the spatial coordinate, which
                        # is just a concatenation of the two parts.
                        # The value is a list of tuples, representing
                        # the normal vector (scaled by the volume of
                        # the facet) at that point.
                        # Each tuple looks like (foo, (i,)); the i'th
                        # component of the scaled normal is foo.

                        # The following line is only valid when the second
                        # shape has spatial dimension 1 (enforced above)
                        Apoint, Avalue = _first_point_pair(Anode)
                        pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}

                        # The following line should be used in the
                        # general case
                        # pt_dict = {Anode.get_point_dict().keys()[0] + Bnode.get_point_dict().keys()[0]: Anode.get_point_dict().values()[0] + [(0.0, (ii,)) for ii in range(len(Anode.get_point_dict().keys()[0]), len(Anode.get_point_dict().keys()[0]) + len(Bnode.get_point_dict().keys()[0]))]}

                        # THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
                        shp = (sd,)
                        nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointScaledNormalEval"))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.PointEdgeTangentEvaluation):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: PointEdgeTangentEval x PointEval
                        # this is very similar to the case above, so comments omitted
                        if len(_first_point(Bnode)) > 1:
                            raise NotImplementedError("PointEdgeTangentEval x PointEval is not yet supported if the second shape has dimension > 1")
                        sd = ref_el.get_spatial_dimension()
                        Apoint, Avalue = _first_point_pair(Anode)
                        pt_dict = {Apoint + _first_point(Bnode): Avalue + [(0.0, (len(Apoint),))]}

                        # THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
                        shp = (sd,)
                        nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "PointEdgeTangent"))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.ComponentPointEvaluation):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: ComponentPointEval x PointEval
                        # the CptPointEval functional requires the component
                        # and the coordinates. very similar to PE x PE case.
                        sd = ref_el.get_spatial_dimension()
                        nodes.append(functional.ComponentPointEvaluation(ref_el, Anode.comp, (sd,), _first_point(Anode) + _first_point(Bnode)))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.FrobeniusIntegralMoment):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: FroIntMom x PointEval
                        sd = ref_el.get_spatial_dimension()
                        pt_dict = {}
                        pt_old = Anode.get_point_dict()
                        for pt in pt_old:
                            pt_dict[pt+_first_point(Bnode)] = pt_old[pt] + [(0.0, sd-1)]
                        # THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
                        shp = (sd,)
                        nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment"))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.IntegralMoment):
                for Bnode in Bnodes:
                    if isinstance(Bnode, functional.PointEvaluation):
                        # case: IntMom x PointEval
                        sd = ref_el.get_spatial_dimension()
                        pt_dict = {}
                        pt_old = Anode.get_point_dict()
                        for pt in pt_old:
                            pt_dict[pt+_first_point(Bnode)] = pt_old[pt]
                        # THE FOLLOWING IS PROBABLY CORRECT BUT UNTESTED
                        shp = (sd,)
                        nodes.append(functional.Functional(ref_el, shp, pt_dict, {}, "IntegralMoment"))
                    else:
                        raise NotImplementedError("unsupported functional type")

            elif isinstance(Anode, functional.Functional):
                # this should catch everything else
                for Bnode in Bnodes:
                    nodes.append(functional.Functional(None, None, None, {}, "Undefined"))
            else:
                raise NotImplementedError("unsupported functional type")

        dual = dual_set.DualSet(nodes, ref_el, entity_ids)

        super(TensorProductElement, self).__init__(ref_el, dual, order, formdegree, mapping)
        # Set up constituent elements
        self.A = A
        self.B = B

        # degree for quadrature rule
        self.polydegree = max(A.degree(), B.degree())

    def degree(self):
        """Return the degree of the (embedding) polynomial space."""
        return self.polydegree

    def get_nodal_basis(self):
        """Return the nodal basis, encoded as a PolynomialSet object,
        for the finite element."""
        raise NotImplementedError("get_nodal_basis not implemented")

    def get_coeffs(self):
        """Return the expansion coefficients for the basis of the
        finite element."""
        raise NotImplementedError("get_coeffs not implemented")

    def tabulate(self, order, points, entity=None):
        """Return tabulated values of derivatives up to given order of
        basis functions at given points."""
        if entity is None:
            entity = (self.ref_el.get_dimension(), 0)
        entity_dim, entity_id = entity

        shape = tuple(len(c.get_topology()[d])
                      for c, d in zip(self.ref_el.cells, entity_dim))
        idA, idB = numpy.unravel_index(entity_id, shape)

        # Factor the entity argument to get entities of the component elements
        entityA_dim, entityB_dim = entity_dim
        entityA = (entityA_dim, idA)
        entityB = (entityB_dim, idB)

        pointsAdim, pointsBdim = [c.get_spatial_dimension()
                                  for c in self.ref_el.construct_subelement(entity_dim).cells]
        pointsA = [point[:pointsAdim] for point in points]
        pointsB = [point[pointsAdim:pointsAdim + pointsBdim] for point in points]

        Asdim = self.A.ref_el.get_spatial_dimension()
        Bsdim = self.B.ref_el.get_spatial_dimension()
        # Note that for entities other than cells, the following
        # tabulations are already appropriately zero-padded so no
        # additional zero padding is required.
        Atab = self.A.tabulate(order, pointsA, entityA)
        Btab = self.B.tabulate(order, pointsB, entityB)
        npoints = len(points)

        # allow 2 scalar-valued FE spaces, or 1 scalar-valued,
        # 1 vector-valued. Combining 2 vector-valued spaces
        # into a tensor-valued space via an outer-product
        # seems to be a sensible general option, but I don't
        # know how to handle the nestedness of the arrays
        # if someone then tries to make a new "tensor finite
        # element" where one component is already a
        # tensor-valued space!
        A_valuedim = len(self.A.value_shape())  # scalar: 0, vector: 1
        B_valuedim = len(self.B.value_shape())  # scalar: 0, vector: 1
        if A_valuedim + B_valuedim > 1:
            raise NotImplementedError("tabulate does not support two vector-valued inputs")
        result = {}
        for i in range(order + 1):
            alphas = mis(Asdim+Bsdim, i)  # thanks, Rob!
            for alpha in alphas:
                if A_valuedim == 0 and B_valuedim == 0:
                    # for each point, get outer product of (A's basis
                    # functions f1, f2, ... evaluated at that point)
                    # with (B's basis functions g1, g2, ... evaluated
                    # at that point). This gives temp[point][f_i][g_j].
                    # Flatten this, so bfs are
                    # in the order f1g1, f1g2, ..., f2g1, f2g2, ...
                    # which is compatible with the entity_dofs order.
                    # We now have temp[point][full basis function]
                    # Transpose this to get temp[bf][point],
                    # and we are done.
                    temp = numpy.array([numpy.outer(
                                       Atab[alpha[0:Asdim]][..., j],
                                       Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
                        .ravel() for j in range(npoints)])
                    result[alpha] = temp.transpose()
                elif A_valuedim == 1 and B_valuedim == 0:
                    # similar to above, except A's basis functions
                    # are now vector-valued. numpy.outer flattens the
                    # array, so it's like taking the OP of
                    # f1_x, f1_y, f2_x, f2_y, ... with g1, g2, ...
                    # this gives us
                    # temp[point][f1x, f1y, f2x, f2y, ...][g_j].
                    # reshape once to get temp[point][f_i][x/y][g_j]
                    # transpose to get temp[point][x/y][f_i][g_j]
                    # reshape to flatten the last two indices, this
                    # gives us temp[point][x/y][full bf_i]
                    # finally, transpose the first and last indices
                    # to get temp[bf_i][x/y][point], and we are done.
                    temp = numpy.array([numpy.outer(
                                       Atab[alpha[0:Asdim]][..., j],
                                       Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
                        for j in range(npoints)])
                    assert temp.shape[1] % 2 == 0
                    temp2 = temp.reshape((temp.shape[0],
                                          temp.shape[1]//2,
                                          2,
                                          temp.shape[2]))\
                        .transpose(0, 2, 1, 3)\
                        .reshape((temp.shape[0], 2, -1))\
                        .transpose(2, 1, 0)
                    result[alpha] = temp2
                elif A_valuedim == 0 and B_valuedim == 1:
                    # as above, with B's functions now vector-valued.
                    # we now do... [numpy.outer ... for ...] gives
                    # temp[point][f_i][g1x,g1y,g2x,g2y,...].
                    # reshape to temp[point][f_i][g_j][x/y]
                    # flatten middle: temp[point][full bf_i][x/y]
                    # transpose to temp[bf_i][x/y][point]
                    temp = numpy.array([numpy.outer(
                        Atab[alpha[0:Asdim]][..., j],
                        Btab[alpha[Asdim:Asdim+Bsdim]][..., j])
                        for j in range(len(Atab[alpha[0:Asdim]][0]))])
                    assert temp.shape[2] % 2 == 0
                    temp2 = temp.reshape((temp.shape[0], temp.shape[1],
                                          temp.shape[2]//2, 2))\
                        .reshape((temp.shape[0], -1, 2))\
                        .transpose(1, 2, 0)
                    result[alpha] = temp2
        return result

    def value_shape(self):
        """Return the value shape of the finite element functions."""
        if len(self.A.value_shape()) == 0 and len(self.B.value_shape()) == 0:
            return ()
        elif len(self.A.value_shape()) == 1 and len(self.B.value_shape()) == 0:
            return (self.A.value_shape()[0],)
        elif len(self.A.value_shape()) == 0 and len(self.B.value_shape()) == 1:
            return (self.B.value_shape()[0],)
        else:
            raise NotImplementedError("value_shape not implemented")

    def dmats(self):
        """Return dmats: expansion coefficients for basis function
        derivatives."""
        raise NotImplementedError("dmats not implemented")

    def get_num_members(self, arg):
        """Return number of members of the expansion set."""
        raise NotImplementedError("get_num_members not implemented")

    def is_nodal(self):
        # This element is nodal iff all factor elements are nodal.
        return all([self.A.is_nodal(), self.B.is_nodal()])
import numpy as np
from FIAT.discontinuous_lagrange import DiscontinuousLagrange
from FIAT.dual_set import DualSet
from FIAT.finite_element import FiniteElement
from FIAT.functional import PointEvaluation
from FIAT.polynomial_set import mis
from FIAT.reference_element import (ufc_simplex, POINT,
                                    LINE, QUADRILATERAL,
                                    TRIANGLE, TETRAHEDRON,
                                    TENSORPRODUCT)
from FIAT.tensor_product import TensorProductElement

        .. note ::

           Performing illegal tabulations on this element will result in either
           a tabulation table of `numpy.nan` arrays (`entity=None` case), or
           insertions of the `TraceError` exception class. This is due to the
           fact that performing cell-wise tabulations, or asking for any order
           of derivative evaluations, are not mathematically well-defined.
        """
        sd = self.ref_el.get_spatial_dimension()
        facet_sd = sd - 1

        # Initializing dictionary with zeros
        phivals = {}
        for i in range(order + 1):
            alphas = mis(sd, i)
            for alpha in alphas:
                phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points)))

        evalkey = (0,) * sd

        # If entity is None, identify facet using numerical tolerance and
        # return the tabulated values
        if entity is None:
            # NOTE: Numerical approximation of the facet id is currently only
            # implemented for simplex reference cells.
            if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]:
                raise NotImplementedError(
                    "Tabulating this element on a %s cell without providing "
                    "an entity is not currently supported." % type(self.ref_el)
                )

            # Attempt to identify which facet (if any) the given points are on
            vertices = self.ref_el.vertices
            coordinates = barycentric_coordinates(points, vertices)
            unique_facet, success = extract_unique_facet(coordinates)

            # If not successful, return NaNs
            if not success:
                for key in phivals:
                    phivals[key] = np.full(shape=(self.space_dimension(), len(points)), fill_value=np.nan)

                return phivals

            # Otherwise, extract non-zero values and insertion indices
            else:
                # Map points to the reference facet
                new_points = map_to_reference_facet(points, vertices, unique_facet)

                # Retrieve values by tabulating the DG element
                element = self.dg_elements[facet_sd]
                nf = element.space_dimension()
                nonzerovals, = element.tabulate(order, new_points).values()
                indices = slice(nf * unique_facet, nf * (unique_facet + 1))

        else:
            entity_dim, _ = entity

            # If the user is directly specifying cell-wise tabulation, return
            # TraceErrors in dict for appropriate handling in the form compiler
            if entity_dim not in self.dg_elements:
                for key in phivals:
                    msg = "The HDivTrace element can only be tabulated on facets."
                    phivals[key] = TraceError(msg)

                return phivals

            else:
                # Retrieve function evaluations (order = 0 case)
                offset = 0
                for facet_dim in sorted(self.dg_elements):
                    element = self.dg_elements[facet_dim]
                    nf = element.space_dimension()
                    num_facets = len(self.ref_el.get_topology()[facet_dim])

                    # Loop over the number of facets until we find a facet
                    # with matching dimension and id
                    for i in range(num_facets):
                        # Found it! Grab insertion indices
                        if (facet_dim, i) == entity:
                            nonzerovals, = element.tabulate(0, points).values()
                            indices = slice(offset, offset + nf)

                        offset += nf

        # If asking for gradient evaluations, insert TraceError in
        # gradient slots
        if order > 0:
            msg = "Gradients on trace elements are not well-defined."
            for key in phivals:
                if key != evalkey:
                    phivals[key] = TraceError(msg)

        # Insert non-zero values in appropriate place
        phivals[evalkey][indices, :] = nonzerovals

        return phivals

    def value_shape(self):
        """Return the value shape of the finite element functions."""
        return ()

    def dmats(self):
        """Return dmats: expansion coefficients for basis function
        derivatives."""
        raise NotImplementedError("dmats not implemented for the trace element.")

    def get_num_members(self, arg):
        """Return number of members of the expansion set."""
        raise NotImplementedError("get_num_members not implemented for the trace element.")

    @staticmethod
    def is_nodal():
        return True


def construct_dg_element(ref_el, degree):
    """Constructs a discontinuous galerkin element of a given degree
    on a particular reference cell.
    """

    if ref_el.get_shape() in [LINE, TRIANGLE]: