POINT = 0
LINE = 1
TRIANGLE = 2
TETRAHEDRON = 3
QUADRILATERAL = 11
HEXAHEDRON = 111
TENSORPRODUCT = 99
    C = numpy.dot(numpy.transpose(qa), qb)
    U = numpy.dot(qa, y)
        A = numpy.zeros([subdim, subdim])
        B = numpy.zeros([celldim, subdim])
        C = B.dot(A)
    A = numpy.reshape(sol[:dim_x * dim_y], (dim_y, dim_x))
def lattice_iter(start, finish, depth):
def make_lattice(verts, n, interior=0):
def linalg_subspace_intersection(A, B):
    def __init__(self, shape, vertices, topology):
    def _key(self):
    def __eq__(self, other):
    def __ne__(self, other):
    def __hash__(self):
    def get_shape(self):
    def get_vertices(self):
    def get_spatial_dimension(self):
    def get_topology(self):
    def get_connectivity(self):
    def get_vertices_of_subcomplex(self, t):
    def get_dimension(self):
    def construct_subelement(self, dimension):
    def get_entity_transform(self, dim, entity_i):
    def compute_normal(self, facet_i):
    def compute_tangents(self, dim, i):
    def compute_normalized_tangents(self, dim, i):
    def compute_edge_tangent(self, edge_i):
    def compute_normalized_edge_tangent(self, edge_i):
    def compute_face_tangents(self, face_i):
    def compute_face_edge_tangents(self, dim, entity_id):
    def make_points(self, dim, entity_id, order):
    def volume(self):
    def volume_of_subcomplex(self, dim, facet_no):
    def compute_scaled_normal(self, facet_i):
    def compute_reference_normal(self, facet_dim, facet_i):
    def get_entity_transform(self, dim, entity):
    def get_dimension(self):
    def get_facet_element(self):
    def construct_subelement(self, dimension):
    def contains_point(self, point, epsilon=0):
    def __init__(self):
    def construct_subelement(self, dimension):
    def __init__(self):
    def get_facet_element(self):
    def __init__(self):
    def __init__(self):
    def get_facet_element(self):
    def __init__(self):
    def compute_normal(self, i):
    def __init__(self):
    def get_facet_element(self):
    def __init__(self):
    def get_facet_element(self):
    def __init__(self):
    def get_facet_element(self):
    def __init__(self):
    def compute_normal(self, i):
    def __init__(self, *cells):
    def _key(self):
    def _split_slices(lengths):
    def get_dimension(self):
    def construct_subelement(self, dimension):
    def get_entity_transform(self, dim, entity_i):
        def transform(point):
    def volume(self):
    def compute_reference_normal(self, facet_dim, facet_i):
    def contains_point(self, point, epsilon=0):
    def __init__(self):
    def get_dimension(self):
    def construct_subelement(self, dimension):
    def get_entity_transform(self, dim, entity_i):
    def volume(self):
    def compute_reference_normal(self, facet_dim, facet_i):
    def contains_point(self, point, epsilon=0):
    def __init__(self):
    def get_dimension(self):
    def construct_subelement(self, dimension):
    def get_entity_transform(self, dim, entity_i):
    def volume(self):
    def compute_reference_normal(self, facet_dim, facet_i):
    def contains_point(self, point, epsilon=0):
def make_affine_mapping(xs, ys):
def default_simplex(spatial_dim):
def ufc_simplex(spatial_dim):
def ufc_cell(cell):
def volume(verts):
def tuple_sum(tree):
def is_hypercube(cell):
def flatten_reference_cube(ref_el):
def flatten_entities(topology_dict):
def compute_unflattening_map(topology_dict):
class Cell(object):
class Simplex(Cell):
class UFCSimplex(Simplex):
class Point(Simplex):
class DefaultLine(Simplex):
class UFCInterval(UFCSimplex):
class DefaultTriangle(Simplex):
class UFCTriangle(UFCSimplex):
class IntrepidTriangle(Simplex):
class DefaultTetrahedron(Simplex):
class IntrepidTetrahedron(Simplex):
class UFCTetrahedron(UFCSimplex):
class TensorProductCell(Cell):
class UFCQuadrilateral(Cell):
class UFCHexahedron(Cell): class Point(Simplex):
    """This is the reference point."""

    def __init__(self):
        verts = ((),)
        topology = {0: {0: (0,)}}
        super(Point, self).__init__(POINT, verts, topology)

    def construct_subelement(self, dimension):
        """Constructs the reference element of a cell subentity
        specified by subelement dimension.

        :arg dimension: subentity dimension (integer). Must be zero.
        """
        assert dimension == 0
        return self class UFCInterval(UFCSimplex):
    """This is the reference interval with vertices (0.0,) and (1.0,)."""

    def __init__(self):
        verts = ((0.0,), (1.0,))
        edges = {0: (0, 1)}
        topology = {0: {0: (0,), 1: (1,)},
                    1: edges}
        super(UFCInterval, self).__init__(LINE, verts, topology) class UFCTriangle(UFCSimplex):
    """This is the reference triangle with vertices (0.0,0.0),
    (1.0,0.0), and (0.0,1.0)."""

    def __init__(self):
        verts = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0))
        edges = {0: (1, 2), 1: (0, 2), 2: (0, 1)}
        faces = {0: (0, 1, 2)}
        topology = {0: {0: (0,), 1: (1,), 2: (2,)},
                    1: edges, 2: faces}
        super(UFCTriangle, self).__init__(TRIANGLE, verts, topology)

    def compute_normal(self, i):
        "UFC consistent normal"
        t = self.compute_tangents(1, i)[0]
        n = numpy.array((t[1], -t[0]))
        return n / numpy.linalg.norm(n) class UFCTetrahedron(UFCSimplex):
    """This is the reference tetrahedron with vertices (0,0,0),
    (1,0,0),(0,1,0), and (0,0,1)."""

    def __init__(self):
        verts = ((0.0, 0.0, 0.0), (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))
        vs = {0: (0, ),
              1: (1, ),
              2: (2, ),
              3: (3, )}
        edges = {0: (2, 3),
                 1: (1, 3),
                 2: (1, 2),
                 3: (0, 3),
                 4: (0, 2),
                 5: (0, 1)}
        faces = {0: (1, 2, 3),
                 1: (0, 2, 3),
                 2: (0, 1, 3),
                 3: (0, 1, 2)}
        tets = {0: (0, 1, 2, 3)}
        topology = {0: vs, 1: edges, 2: faces, 3: tets}
        super(UFCTetrahedron, self).__init__(TETRAHEDRON, verts, topology)

    def compute_normal(self, i):
        "UFC consistent normals."
        t = self.compute_tangents(2, i)
        n = numpy.cross(t[0], t[1])
        return -2.0 * n / numpy.linalg.norm(n) class Lagrange(finite_element.CiarletElement):
    """The Lagrange finite element.  It is what it is."""

    def __init__(self, ref_el, degree):
        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
        dual = LagrangeDualSet(ref_el, degree)
        formdegree = 0  # 0-form
        super(Lagrange, self).__init__(poly_set, dual, degree, formdegree) def DiscontinuousLagrange(ref_el, degree):
    if degree == 0:
        return P0.P0(ref_el)
    else:
        return HigherOrderDiscontinuousLagrange(ref_el, degree) def DiscontinuousTaylor(ref_el, degree):
    if degree == 0:
        return P0.P0(ref_el)
    else:
        return HigherOrderDiscontinuousTaylor(ref_el, degree) class P0(finite_element.CiarletElement):
    def __init__(self, ref_el):
        poly_set = polynomial_set.ONPolynomialSet(ref_el, 0)
        dual = P0Dual(ref_el)
        degree = 0
        formdegree = ref_el.get_spatial_dimension()  # n-form
        super(P0, self).__init__(poly_set, dual, degree, formdegree) class CrouzeixRaviart(finite_element.CiarletElement):
    """The Crouzeix-Raviart finite element:

    K:                 Triangle/Tetrahedron
    Polynomial space:  P_1
    Dual basis:        Evaluation at facet midpoints
    """

    def __init__(self, cell, degree):

        # Crouzeix Raviart is only defined for polynomial degree == 1
        if not (degree == 1):
            raise Exception("Crouzeix-Raviart only defined for degree 1")

        # Construct polynomial spaces, dual basis and initialize
        # FiniteElement
        space = polynomial_set.ONPolynomialSet(cell, 1)
        dual = CrouzeixRaviartDualSet(cell, 1)
        super(CrouzeixRaviart, self).__init__(space, dual, 1) class RaviartThomas(finite_element.CiarletElement):
    """
    The Raviart Thomas element

    :arg ref_el: The reference element.
    :arg k: The degree.
    :arg variant: optional variant specifying the types of nodes.

    variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
    "point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
    convergence order in the H(div)-norm
    "integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
    polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
    when you want to have (nearly) divergence-preserving interpolation.
    "integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
    """

    def __init__(self, ref_el, k, variant=None):

        degree = k - 1

        (variant, quad_deg) = check_format_variant(variant, degree, "Raviart Thomas")

        poly_set = RTSpace(ref_el, degree)
        dual = RTDualSet(ref_el, degree, variant, quad_deg)
        formdegree = ref_el.get_spatial_dimension() - 1  # (n-1)-form
        super(RaviartThomas, self).__init__(poly_set, dual, degree, formdegree,
                                            mapping="contravariant piola") class DiscontinuousRaviartThomas(finite_element.CiarletElement):
    """The discontinuous Raviart-Thomas finite element"""

    def __init__(self, ref_el, q):

        degree = q - 1
        poly_set = RTSpace(ref_el, degree)
        dual = DRTDualSet(ref_el, degree)
        super(DiscontinuousRaviartThomas, self).__init__(poly_set, dual, degree,
                                                         mapping="contravariant piola") class BrezziDouglasMarini(finite_element.CiarletElement):
    """
    The BDM element

    :arg ref_el: The reference element.
    :arg k: The degree.
    :arg variant: optional variant specifying the types of nodes.

    variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
    "point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
    convergence order in the H(div)-norm
    "integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
    polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
    when you want to have (nearly) divergence-preserving interpolation.
    "integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
    """

    def __init__(self, ref_el, k, variant=None):

        (variant, quad_deg) = check_format_variant(variant, k, "BDM")

        if k < 1:
            raise Exception("BDM_k elements only valid for k >= 1")

        sd = ref_el.get_spatial_dimension()
        poly_set = polynomial_set.ONPolynomialSet(ref_el, k, (sd, ))
        dual = BDMDualSet(ref_el, k, variant, quad_deg)
        formdegree = sd - 1  # (n-1)-form
        super(BrezziDouglasMarini, self).__init__(poly_set, dual, k, formdegree,
                                                  mapping="contravariant piola") class MixedElement(FiniteElement):
    """A FIAT-like representation of a mixed element.

    :arg elements: An iterable of FIAT elements.
    :arg ref_el: The reference element (optional).

    This object offers tabulation of the concatenated basis function
    tables along with an entity_dofs dict."""
    def __init__(self, elements, ref_el=None):
        elements = tuple(elements)

        cells = set(e.get_reference_element() for e in elements)
        if ref_el is not None:
            cells.add(ref_el)
        ref_el, = cells

        # These functionals are absolutely wrong, they all map from
        # functions of the wrong shape, and potentially of different
        # shapes.  However, they are wrong precisely as FFC hacks
        # expect them to be. :(
        nodes = [L for e in elements for L in e.dual_basis()]

        entity_dofs = concatenate_entity_dofs(ref_el, elements)

        dual = DualSet(nodes, ref_el, entity_dofs)
        super(MixedElement, self).__init__(ref_el, dual, None, mapping=None)
        self._elements = elements

    def elements(self):
        return self._elements

    def num_sub_elements(self):
        return len(self._elements)

    def value_shape(self):
        return (sum(numpy.prod(e.value_shape(), dtype=int) for e in self.elements()), )

    def mapping(self):
        return [m for e in self._elements for m in e.mapping()]

    def get_nodal_basis(self):
        raise NotImplementedError("get_nodal_basis not implemented")

    def tabulate(self, order, points, entity=None):
        """Tabulate a mixed element by appropriately splatting
        together the tabulation of the individual elements.
        """
        shape = (self.space_dimension(),) + self.value_shape() + (len(points),)

        output = {}

        sub_dims = [0] + list(e.space_dimension() for e in self.elements())
        sub_cmps = [0] + list(numpy.prod(e.value_shape(), dtype=int)
                              for e in self.elements())
        irange = numpy.cumsum(sub_dims)
        crange = numpy.cumsum(sub_cmps)

        for i, e in enumerate(self.elements()):
            table = e.tabulate(order, points, entity)

            for d, tab in table.items():
                try:
                    arr = output[d]
                except KeyError:
                    arr = numpy.zeros(shape, dtype=tab.dtype)
                    output[d] = arr

                ir = irange[i:i+2]
                cr = crange[i:i+2]
                tab = tab.reshape(ir[1] - ir[0], cr[1] - cr[0], -1)
                arr[slice(*ir), slice(*cr)] = tab

        return output

    def is_nodal(self):
        """True if primal and dual bases are orthogonal."""
        return all(e.is_nodal() for e in self._elements) class Nedelec(finite_element.CiarletElement):
    """
    Nedelec finite element

    :arg ref_el: The reference element.
    :arg k: The degree.
    :arg variant: optional variant specifying the types of nodes.

    variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
    "point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
    convergence order in the H(curl)-norm
    "integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
    polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
    when you want to have (nearly) curl-preserving interpolation.
    "integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
    """

    def __init__(self, ref_el, k, variant=None):

        degree = k - 1

        (variant, quad_deg) = check_format_variant(variant, degree, "Nedelec")

        if ref_el.get_spatial_dimension() == 3:
            poly_set = NedelecSpace3D(ref_el, degree)
            dual = NedelecDual3D(ref_el, degree, variant, quad_deg)
        elif ref_el.get_spatial_dimension() == 2:
            poly_set = NedelecSpace2D(ref_el, degree)
            dual = NedelecDual2D(ref_el, degree, variant, quad_deg)
        else:
            raise Exception("Not implemented")
        formdegree = 1  # 1-form
        super(Nedelec, self).__init__(poly_set, dual, degree, formdegree,
                                      mapping="covariant piola") class NedelecSecondKind(CiarletElement):
    """
    The H(curl) Nedelec elements of the second kind on triangles and
    tetrahedra: the polynomial space described by the full polynomials
    of degree k, with a suitable set of degrees of freedom to ensure
    H(curl) conformity.

    :arg ref_el: The reference element.
    :arg k: The degree.
    :arg variant: optional variant specifying the types of nodes.

    variant can be chosen from ["point", "integral", "integral(quadrature_degree)"]
    "point" -> dofs are evaluated by point evaluation. Note that this variant has suboptimal
    convergence order in the H(curl)-norm
    "integral" -> dofs are evaluated by quadrature rule. The quadrature degree is chosen to integrate
    polynomials of degree 5*k so that most expressions will be interpolated exactly. This is important
    when you want to have (nearly) curl-preserving interpolation.
    "integral(quadrature_degree)" -> dofs are evaluated by quadrature rule of degree quadrature_degree
    """

    def __init__(self, cell, k, variant=None):

        (variant, quad_deg) = check_format_variant(variant, k, "Nedelec Second Kind")

        # Check degree
        assert k >= 1, "Second kind Nedelecs start at 1!"

        # Get dimension
        d = cell.get_spatial_dimension()

        # Construct polynomial basis for d-vector fields
        Ps = ONPolynomialSet(cell, k, (d, ))

        # Construct dual space
        Ls = NedelecSecondKindDual(cell, k, variant, quad_deg)

        # Set form degree
        formdegree = 1  # 1-form

        # Set mapping
        mapping = "covariant piola"

        # Call init of super-class
        super(NedelecSecondKind, self).__init__(Ps, Ls, k, formdegree, mapping=mapping) class Regge(CiarletElement):
    """The generalized Regge elements for symmetric-matrix-valued functions.
       REG(r) in dimension n is the space of polynomial symmetric-matrix-valued
       functions of degree r or less with tangential-tangential continuity.
    """
    def __init__(self, cell, degree):
        assert degree >= 0, "Regge start at degree 0!"
        # shape functions
        Ps = ONSymTensorPolynomialSet(cell, degree)
        # degrees of freedom
        Ls = ReggeDual(cell, degree)
        # mapping under affine transformation
        mapping = "double covariant piola"

        super(Regge, self).__init__(Ps, Ls, degree, mapping=mapping) class HDivTrace(FiniteElement):
    """Class implementing the trace of hdiv elements. This class
    is a stand-alone element family that produces a DG-facet field.
    This element is what's produced after performing the trace
    operation on an existing H(Div) element.

    This element is also known as the discontinuous trace field that
    arises in several DG formulations.
    """

    def __init__(self, ref_el, degree):
        """Constructor for the HDivTrace element.

        :arg ref_el: A reference element, which may be a tensor product
                     cell.
        :arg degree: The degree of approximation. If on a tensor product
                     cell, then provide a tuple of degrees if you want
                     varying degrees.
        """
        sd = ref_el.get_spatial_dimension()
        if sd in (0, 1):
            raise ValueError("Cannot take the trace of a %d-dim cell." % sd)

        # Store the degrees if on a tensor product cell
        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), (
                "Number of specified degrees must be equal to the number of cells."
            )
        else:
            if ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON, QUADRILATERAL]:
                raise NotImplementedError(
                    "Trace element on a %s not implemented" % type(ref_el)
                )
            # Cannot have varying degrees for these reference cells
            if isinstance(degree, tuple):
                raise ValueError("Must have a tensor product cell if providing multiple degrees")

        # Initialize entity dofs and construct the DG elements
        # for the facets
        facet_sd = sd - 1
        dg_elements = {}
        entity_dofs = {}
        topology = ref_el.get_topology()
        for top_dim, entities in topology.items():
            cell = ref_el.construct_subelement(top_dim)
            entity_dofs[top_dim] = {}

            # We have a facet entity!
            if cell.get_spatial_dimension() == facet_sd:
                dg_elements[top_dim] = construct_dg_element(cell, degree)
            # Initialize
            for entity in entities:
                entity_dofs[top_dim][entity] = []

        # Compute the dof numbering for all facet entities
        # and extract nodes
        offset = 0
        pts = []
        for facet_dim in sorted(dg_elements):
            element = dg_elements[facet_dim]
            nf = element.space_dimension()
            num_facets = len(topology[facet_dim])

            for i in range(num_facets):
                entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
                offset += nf

                # Run over nodes and collect the points for point evaluations
                for dof in element.dual_basis():
                    facet_pt, = dof.get_point_dict()
                    transform = ref_el.get_entity_transform(facet_dim, i)
                    pts.append(tuple(transform(facet_pt)))

        # Setting up dual basis - only point evaluations
        nodes = [PointEvaluation(ref_el, pt) for pt in pts]
        dual = DualSet(nodes, ref_el, entity_dofs)

        # Degree of the element
        deg = max([e.degree() for e in dg_elements.values()])

        super(HDivTrace, self).__init__(ref_el, dual, order=deg,
                                        formdegree=facet_sd,
                                        mapping="affine")

        # Set up facet elements
        self.dg_elements = dg_elements

        # Degree for quadrature rule
        self.polydegree = deg

    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 for the trace element.")

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

    def tabulate(self, order, points, entity=None):
        """Return tabulated values of derivatives up to a 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``,
                     tabulated values are computed by geometrically
                     approximating which facet the points are on.

        .. 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 map_to_reference_facet(points, vertices, facet):
    """Given a set of points and vertices describing a facet of a simplex in n-dimensional
    coordinates (where the points lie on the facet), map the points to the reference simplex
    of dimension (n-1).

    :arg points: A set of points in n-D.
    :arg vertices: A set of vertices describing a facet of a simplex in n-D.
    :arg facet: Integer representing the facet number.
    """

    # Compute the barycentric coordinates of the points with respect to the
    # full physical simplex
    all_coords = barycentric_coordinates(points, vertices)

    # Extract vertices of the reference facet
    reference_facet_simplex = ufc_simplex(len(vertices) - 2)
    reference_vertices = reference_facet_simplex.get_vertices()

    reference_points = []
    for (i, coords) in enumerate(all_coords):
        # Extract the correct subset of barycentric coordinates since we know
        # which facet we are on
        new_coords = [coords[j] for j in range(len(coords)) if j != facet]

        # Evaluate the reference coordinate of a point in barycentric coordinates
        reference_pt = sum(np.asarray(reference_vertices[j]) * new_coords[j]
                           for j in range(len(new_coords)))

        reference_points += [reference_pt]
    return reference_points class HellanHerrmannJohnson(CiarletElement):
    """The definition of Hellan-Herrmann-Johnson element. It is defined only in
       dimension 2. It consists of piecewise polynomial symmetric-matrix-valued
       functions of degree r or less with normal-normal continuity.
    """
    def __init__(self, cell, degree):
        assert degree >= 0, "Hellan-Herrmann-Johnson starts at degree 0!"
        # shape functions
        Ps = ONSymTensorPolynomialSet(cell, degree)
        # degrees of freedom
        Ls = HellanHerrmannJohnsonDual(cell, degree)
        # mapping under affine transformation
        mapping = "double contravariant piola"

        super(HellanHerrmannJohnson, self).__init__(Ps, Ls, degree,
                                                    mapping=mapping) class BrezziDouglasFortinMarini(finite_element.CiarletElement):
    """The BDFM element"""

    def __init__(self, ref_el, degree):

        if degree != 2:
            raise Exception("BDFM_k elements only valid for k == 2")

        poly_set = BDFMSpace(ref_el, degree)
        dual = BDFMDualSet(ref_el, degree - 1)
        formdegree = ref_el.get_spatial_dimension() - 1
        super(BrezziDouglasFortinMarini, self).__init__(poly_set, dual, degree, formdegree,
                                                        mapping="contravariant piola") class GaussLegendre(finite_element.CiarletElement):
    """1D discontinuous element with nodes at the Gauss-Legendre points."""
    def __init__(self, ref_el, degree):
        if ref_el.shape != LINE:
            raise ValueError("Gauss-Legendre elements are only defined in one dimension.")
        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
        dual = GaussLegendreDualSet(ref_el, degree)
        formdegree = ref_el.get_spatial_dimension()  # n-form
        super(GaussLegendre, self).__init__(poly_set, dual, degree, formdegree)

    def tabulate(self, order, points, entity=None):
        # This overrides the default with a more numerically stable algorithm

        if entity is None:
            entity = (self.ref_el.get_dimension(), 0)

        entity_dim, entity_id = entity
        transform = self.ref_el.get_entity_transform(entity_dim, entity_id)

        xsrc = []
        for node in self.dual.nodes:
            # Assert singleton point for each node.
            (pt,), = node.get_point_dict().keys()
            xsrc.append(pt)
        xsrc = numpy.asarray(xsrc)
        xdst = numpy.array(list(map(transform, points))).flatten()
        return barycentric_interpolation(xsrc, xdst, order=order) class GaussLobattoLegendre(finite_element.CiarletElement):
    """1D continuous element with nodes at the Gauss-Lobatto points."""
    def __init__(self, ref_el, degree):
        if ref_el.shape != LINE:
            raise ValueError("Gauss-Lobatto-Legendre elements are only defined in one dimension.")
        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
        dual = GaussLobattoLegendreDualSet(ref_el, degree)
        formdegree = 0  # 0-form
        super(GaussLobattoLegendre, self).__init__(poly_set, dual, degree, formdegree)

    def tabulate(self, order, points, entity=None):
        # This overrides the default with a more numerically stable algorithm

        if entity is None:
            entity = (self.ref_el.get_dimension(), 0)

        entity_dim, entity_id = entity
        transform = self.ref_el.get_entity_transform(entity_dim, entity_id)

        xsrc = []
        for node in self.dual.nodes:
            # Assert singleton point for each node.
            (pt,), = node.get_point_dict().keys()
            xsrc.append(pt)
        xsrc = numpy.asarray(xsrc)
        xdst = numpy.array(list(map(transform, points))).flatten()
        return barycentric_interpolation(xsrc, xdst, order=order) class RestrictedElement(CiarletElement):
    """Restrict given element to specified list of dofs."""

    def __init__(self, element, indices=None, restriction_domain=None):
        '''For sake of argument, indices overrides restriction_domain'''

        if not (indices or restriction_domain):
            raise RuntimeError("Either indices or restriction_domain must be passed in")

        if not indices:
            indices = _get_indices(element, restriction_domain)

        if isinstance(indices, str):
            raise RuntimeError("variable 'indices' was a string; did you forget to use a keyword?")

        if len(indices) == 0:
            raise ValueError("No point in creating empty RestrictedElement.")

        self._element = element
        self._indices = indices

        # Fetch reference element
        ref_el = element.get_reference_element()

        # Restrict primal set
        poly_set = element.get_nodal_basis().take(indices)

        # Restrict dual set
        dof_counter = 0
        entity_ids = {}
        nodes = []
        nodes_old = element.dual_basis()
        for d, entities in element.entity_dofs().items():
            entity_ids[d] = {}
            for entity, dofs in entities.items():
                entity_ids[d][entity] = []
                for dof in dofs:
                    if dof not in indices:
                        continue
                    entity_ids[d][entity].append(dof_counter)
                    dof_counter += 1
                    nodes.append(nodes_old[dof])
        assert dof_counter == len(indices)
        dual = DualSet(nodes, ref_el, entity_ids)

        # Restrict mapping
        mapping_old = element.mapping()
        mapping_new = [mapping_old[dof] for dof in indices]
        assert all(e_mapping == mapping_new[0] for e_mapping in mapping_new)

        # Call constructor of CiarletElement
        super(RestrictedElement, self).__init__(poly_set, dual, 0, element.get_formdegree(), mapping_new[0]) 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()]) class FlattenedDimensions(FiniteElement):
    """A wrapper class that flattens entity dimensions of a FIAT element defined
    on a TensorProductCell to one with quadrilateral/hexahedron entities.
    TensorProductCell has dimension defined as a tuple of factor element dimensions
    (i, j) in 2D  and (i, j, k) in 3D.
    Flattened dimension is a sum of the tuple elements."""

    def __init__(self, element):

        nodes = element.dual.nodes
        dim = element.ref_el.get_spatial_dimension()

        if dim == 2:
            ref_el = UFCQuadrilateral()
        elif dim == 3:
            ref_el = UFCHexahedron()
        else:
            raise ValueError("Illegal element dimension %s" % dim)

        entity_ids = element.dual.entity_ids

        flat_entity_ids = flatten_entities(entity_ids)
        dual = DualSet(nodes, ref_el, flat_entity_ids)
        super(FlattenedDimensions, self).__init__(ref_el, dual, element.get_order(), element.get_formdegree(), element._mapping)
        self.element = element

        # Construct unflattening map for passing correct values to tabulate()
        self.unflattening_map = compute_unflattening_map(self.element.ref_el.get_topology())

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

    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.get_reference_element().get_spatial_dimension(), 0)

        # Entity is provided in flattened form (d, i)
        # Appropriate product entity is taken from the unflattening_map dict
        entity_dim, entity_id = entity
        product_entity = self.unflattening_map[(entity_dim, entity_id)]

        return self.element.tabulate(order, points, product_entity)

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

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

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

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

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

    def is_nodal(self):
        # This element is nodal iff unflattened element is nodal.
        return self.element.is_nodal() def Hdiv(element):
    if not isinstance(element, TensorProductElement):
        raise NotImplementedError

    if element.A.get_formdegree() is None or element.B.get_formdegree() is None:
        raise ValueError("form degree of sub-element was None (not set during initialisation), Hdiv cannot be done without this information")
    formdegree = element.A.get_formdegree() + element.B.get_formdegree()
    if formdegree != element.get_reference_element().get_spatial_dimension() - 1:
        raise ValueError("Tried to use Hdiv on a non-(n-1)-form element")

    newelement = TensorProductElement(element.A, element.B)  # make a copy to return

    # redefine value_shape()
    def value_shape(self):
        "Return the value shape of the finite element functions."
        return (self.get_reference_element().get_spatial_dimension(),)
    newelement.value_shape = types.MethodType(value_shape, newelement)

    # store old _mapping
    newelement._oldmapping = newelement._mapping

    # redefine _mapping
    newelement._mapping = "contravariant piola"

    # store formdegree
    newelement.formdegree = formdegree

    # redefine tabulate
    newelement.old_tabulate = newelement.tabulate

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

        # don't duplicate what the old function does fine...
        old_result = self.old_tabulate(order, points, entity)
        new_result = {}
        sd = self.get_reference_element().get_spatial_dimension()
        for alpha in old_result.keys():
            temp_old = old_result[alpha]

            if self._oldmapping == "affine":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[1]), dtype=temp_old.dtype)
                # both constituents affine, i.e., they were 0 forms or n-forms.
                # to sum to n-1, we must have "0-form on an interval" crossed
                # with something discontinuous.
                # look for the (continuous) 0-form, and put the value there
                if self.A.get_formdegree() == 0:
                    # first element, so (-x, 0, ...)
                    # Sign flip to ensure that a positive value of the node
                    # means a vector field having a direction "to the left"
                    # relative to direction in which the nodes are placed on an
                    # edge in case of higher-order schemes.
                    # This is required for unstructured quadrilateral meshes.
                    temp[:, 0, :] = -temp_old[:, :]
                elif self.B.get_formdegree() == 0:
                    # second element, so (..., 0, x)
                    temp[:, -1, :] = temp_old[:, :]
                else:
                    raise Exception("Hdiv affine/affine form degrees broke")

            elif self._oldmapping == "contravariant piola":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[2]), dtype=temp_old.dtype)
                Asd = self.A.get_reference_element().get_spatial_dimension()
                # one component is affine, one is contravariant piola
                # the affine one must be an n-form, hence discontinuous
                # this component/these components get zeroed out
                if element.A.mapping()[0] == "contravariant piola":
                    # first element, so (x1, ..., xn, 0, ...)
                    temp[:, :Asd, :] = temp_old[:, :, :]
                elif element.B.mapping()[0] == "contravariant piola":
                    # second element, so (..., 0, x1, ..., xn)
                    temp[:, Asd:, :] = temp_old[:, :, :]
                else:
                    raise ValueError("Hdiv contravariant piola couldn't find an existing ConPi subelement")

            elif self._oldmapping == "covariant piola":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[2]), dtype=temp_old.dtype)
                # one component is affine, one is covariant piola
                # the affine one must be an n-form, hence discontinuous
                # this component/these components get zeroed out
                # the remaining part gets perped
                if element.A.mapping()[0] == "covariant piola":
                    Asd = self.A.get_reference_element().get_spatial_dimension()
                    if not Asd == 2:
                        raise ValueError("Must be 2d shape to automatically convert covariant to contravariant")
                    temp_perp = numpy.zeros(temp_old.shape, dtype=temp_old.dtype)
                    # first element, so (x2, -x1, 0, ...)
                    temp_perp[:, 0, :] = temp_old[:, 1, :]
                    temp_perp[:, 1, :] = -temp_old[:, 0, :]
                    temp[:, :Asd, :] = temp_perp[:, :, :]
                elif element.B.mapping()[0] == "covariant piola":
                    Bsd = self.B.get_reference_element().get_spatial_dimension()
                    if not Bsd == 2:
                        raise ValueError("Must be 2d shape to automatically convert covariant to contravariant")
                    temp_perp = numpy.zeros(temp_old.shape, dtype=temp_old.dtype)
                    # second element, so (..., 0, x2, -x1)
                    temp_perp[:, 0, :] = temp_old[:, 1, :]
                    temp_perp[:, 1, :] = -temp_old[:, 0, :]
                    temp[:, Asd:, :] = temp_old[:, :, :]
                else:
                    raise ValueError("Hdiv covariant piola couldn't find an existing CovPi subelement")
            new_result[alpha] = temp
        return new_result

    newelement.tabulate = types.MethodType(tabulate, newelement)

    # splat any PointEvaluation functionals.
    # they become a nasty mix of internal and external component DOFs
    if newelement._oldmapping == "affine":
        oldnodes = newelement.dual.nodes
        newnodes = []
        for node in oldnodes:
            if isinstance(node, functional.PointEvaluation):
                newnodes.append(functional.Functional(None, None, None, {}, "Undefined"))
            else:
                newnodes.append(node)
        newelement.dual.nodes = newnodes

    return newelement def Hcurl(element):
    if not isinstance(element, TensorProductElement):
        raise NotImplementedError

    if element.A.get_formdegree() is None or element.B.get_formdegree() is None:
        raise ValueError("form degree of sub-element was None (not set during initialisation), Hcurl cannot be done without this information")
    formdegree = element.A.get_formdegree() + element.B.get_formdegree()
    if not (formdegree == 1):
        raise ValueError("Tried to use Hcurl on a non-1-form element")

    newelement = TensorProductElement(element.A, element.B)  # make a copy to return

    # redefine value_shape()
    def value_shape(self):
        "Return the value shape of the finite element functions."
        return (self.get_reference_element().get_spatial_dimension(),)
    newelement.value_shape = types.MethodType(value_shape, newelement)

    # store old _mapping
    newelement._oldmapping = newelement._mapping

    # redefine _mapping
    newelement._mapping = "covariant piola"

    # store formdegree
    newelement.formdegree = formdegree

    # redefine tabulate
    newelement.old_tabulate = newelement.tabulate

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

        # don't duplicate what the old function does fine...
        old_result = self.old_tabulate(order, points, entity)
        new_result = {}
        sd = self.get_reference_element().get_spatial_dimension()
        for alpha in old_result.keys():
            temp_old = old_result[alpha]

            if self._oldmapping == "affine":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[1]), dtype=temp_old.dtype)
                # both constituents affine, i.e., they were 0 forms or n-forms.
                # to sum to 1, we must have "1-form on an interval" crossed with
                # a bunch of 0-forms (continuous).
                # look for the 1-form, and put the value in the other place
                if self.A.get_formdegree() == 1:
                    # first element, so (x, 0, ...)
                    # No sign flip here, nor at the other branch, to ensure that
                    # a positive value of the node means a vector field having
                    # the same direction as the direction in which the nodes are
                    # placed on an edge in case of higher-order schemes.
                    # This is required for unstructured quadrilateral meshes.
                    temp[:, 0, :] = temp_old[:, :]
                elif self.B.get_formdegree() == 1:
                    # second element, so (..., 0, x)
                    temp[:, -1, :] = temp_old[:, :]
                else:
                    raise Exception("Hcurl affine/affine form degrees broke")

            elif self._oldmapping == "covariant piola":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[2]), dtype=temp_old.dtype)
                Asd = self.A.get_reference_element().get_spatial_dimension()
                # one component is affine, one is covariant piola
                # the affine one must be an 0-form, hence continuous
                # this component/these components get zeroed out
                if element.A.mapping()[0] == "covariant piola":
                    # first element, so (x1, ..., xn, 0, ...)
                    temp[:, :Asd, :] = temp_old[:, :, :]
                elif element.B.mapping()[0] == "covariant piola":
                    # second element, so (..., 0, x1, ..., xn)
                    temp[:, Asd:, :] = temp_old[:, :, :]
                else:
                    raise ValueError("Hdiv contravariant piola couldn't find an existing ConPi subelement")

            elif self._oldmapping == "contravariant piola":
                temp = numpy.zeros((temp_old.shape[0], sd, temp_old.shape[2]), dtype=temp_old.dtype)
                # one component is affine, one is contravariant piola
                # the affine one must be an 0-form, hence continuous
                # this component/these components get zeroed out
                # the remaining part gets perped
                if element.A.mapping()[0] == "contravariant piola":
                    Asd = self.A.get_reference_element().get_spatial_dimension()
                    if not Asd == 2:
                        raise ValueError("Must be 2d shape to automatically convert contravariant to covariant")
                    temp_perp = numpy.zeros(temp_old.shape, dtype=temp_old.dtype)
                    # first element, so (-x2, x1, 0, ...)
                    temp_perp[:, 0, :] = -temp_old[:, 1, :]
                    temp_perp[:, 1, :] = temp_old[:, 0, :]
                    temp[:, :Asd, :] = temp_perp[:, :, :]
                elif element.B.mapping()[0] == "contravariant piola":
                    Bsd = self.B.get_reference_element().get_spatial_dimension()
                    if not Bsd == 2:
                        raise ValueError("Must be 2d shape to automatically convert contravariant to covariant")
                    temp_perp = numpy.zeros(temp_old.shape, dtype=temp_old.dtype)
                    # second element, so (..., 0, -x2, x1)
                    temp_perp[:, 0, :] = -temp_old[:, 1, :]
                    temp_perp[:, 1, :] = temp_old[:, 0, :]
                    temp[:, Asd:, :] = temp_old[:, :, :]
                else:
                    raise ValueError("Hcurl contravariant piola couldn't find an existing CovPi subelement")
            new_result[alpha] = temp
        return new_result

    newelement.tabulate = types.MethodType(tabulate, newelement)

    # splat any PointEvaluation functionals.
    # they become a nasty mix of internal and external component DOFs
    if newelement._oldmapping == "affine":
        oldnodes = newelement.dual.nodes
        newnodes = []
        for node in oldnodes:
            if isinstance(node, functional.PointEvaluation):
                newnodes.append(functional.Functional(None, None, None, {}, "Undefined"))
            else:
                newnodes.append(node)
        newelement.dual.nodes = newnodes

    return newelement class Argyris(finite_element.CiarletElement):
    """The Argyris finite element."""

    def __init__(self, ref_el, degree):
        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
        dual = ArgyrisDualSet(ref_el, degree)
        super(Argyris, self).__init__(poly_set, dual, degree) class QuinticArgyris(finite_element.CiarletElement):
    """The Argyris finite element."""

    def __init__(self, ref_el):
        poly_set = polynomial_set.ONPolynomialSet(ref_el, 5)
        dual = QuinticArgyrisDualSet(ref_el)
        super().__init__(poly_set, dual, 5) class CubicHermite(finite_element.CiarletElement):
    """The cubic Hermite finite element.  It is what it is."""

    def __init__(self, ref_el, deg=3):
        assert deg == 3
        poly_set = polynomial_set.ONPolynomialSet(ref_el, 3)
        dual = CubicHermiteDualSet(ref_el)

        super(CubicHermite, self).__init__(poly_set, dual, 3) class Morley(finite_element.CiarletElement):
    """The Morley finite element."""

    def __init__(self, ref_el):
        poly_set = polynomial_set.ONPolynomialSet(ref_el, 2)
        dual = MorleyDualSet(ref_el)
        super().__init__(poly_set, dual, 2) class Bubble(CodimBubble):
    """The bubble finite element: the dofs of the Lagrange FE in the interior of the cell"""

    def __init__(self, ref_el, degree):
        super(Bubble, self).__init__(ref_el, degree, codim=0) class EnrichedElement(FiniteElement):
    """Class implementing a finite element that combined the degrees of freedom
    of two existing finite elements.

    This is an implementation which does not care about orthogonality of
    primal and dual basis.
    """

    def __init__(self, *elements):
        # Firstly, check it makes sense to enrich.  Elements must have:
        # - same reference element
        # - same mapping
        # - same value shape
        if len(set(e.get_reference_element() for e in elements)) > 1:
            raise ValueError("Elements must be defined on the same reference element")
        if len(set(m for e in elements for m in e.mapping())) > 1:
            raise ValueError("Elements must have same mapping")
        if len(set(e.value_shape() for e in elements)) > 1:
            raise ValueError("Elements must have the same value shape")

        # order is at least max, possibly more, though getting this
        # right isn't important AFAIK
        order = max(e.get_order() for e in elements)
        # form degree is essentially max (not true for Hdiv/Hcurl,
        # but this will raise an error above anyway).
        # E.g. an H^1 function enriched with an L^2 is now just L^2.
        if any(e.get_formdegree() is None for e in elements):
            formdegree = None
        else:
            formdegree = max(e.get_formdegree() for e in elements)

        # set up reference element and mapping, following checks above
        ref_el, = set(e.get_reference_element() for e in elements)
        mapping, = set(m for e in elements for m in e.mapping())

        # set up entity_ids - for each geometric entity, just concatenate
        # the entities of the constituent elements
        entity_ids = concatenate_entity_dofs(ref_el, elements)

        # set up dual basis - just concatenation
        nodes = list(chain.from_iterable(e.dual_basis() for e in elements))
        dual = DualSet(nodes, ref_el, entity_ids)

        super(EnrichedElement, self).__init__(ref_el, dual, order, formdegree, mapping)

        # required degree (for quadrature) is definitely max
        self.polydegree = max(e.degree() for e in elements)

        # Store subelements
        self._elements = elements

    def elements(self):
        "Return reference to original subelements"
        return self._elements

    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."""

        num_components = numpy.prod(self.value_shape())
        table_shape = (self.space_dimension(), num_components, len(points))

        table = {}
        irange = slice(0)
        for element in self._elements:

            etable = element.tabulate(order, points, entity)
            irange = slice(irange.stop, irange.stop + element.space_dimension())

            # Insert element table into table
            for dtuple in etable.keys():

                if dtuple not in table:
                    if num_components == 1:
                        table[dtuple] = numpy.zeros((self.space_dimension(), len(points)),
                                                    dtype=etable[dtuple].dtype)
                    else:
                        table[dtuple] = numpy.zeros(table_shape,
                                                    dtype=etable[dtuple].dtype)

                table[dtuple][irange][:] = etable[dtuple]

        return table

    def value_shape(self):
        """Return the value shape of the finite element functions."""
        result, = set(e.value_shape() for e in self._elements)
        return result

    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") class NodalEnrichedElement(CiarletElement):
    """NodalEnriched element is a direct sum of a sequence of
    finite elements. Dual basis is reorthogonalized to the
    primal basis for nodality.

    The following is equivalent:
        * the constructor is well-defined,
        * the resulting element is unisolvent and its basis is nodal,
        * the supplied elements are unisolvent with nodal basis and
          their primal bases are mutually linearly independent,
        * the supplied elements are unisolvent with nodal basis and
          their dual bases are mutually linearly independent.
    """

    def __init__(self, *elements):

        # Test elements are nodal
        if not all(e.is_nodal() for e in elements):
            raise ValueError("Not all elements given for construction "
                             "of NodalEnrichedElement are nodal")

        # Extract common data
        ref_el = elements[0].get_reference_element()
        expansion_set = elements[0].get_nodal_basis().get_expansion_set()
        degree = min(e.get_nodal_basis().get_degree() for e in elements)
        embedded_degree = max(e.get_nodal_basis().get_embedded_degree()
                              for e in elements)
        order = max(e.get_order() for e in elements)
        mapping = elements[0].mapping()[0]
        formdegree = None if any(e.get_formdegree() is None for e in elements) \
            else max(e.get_formdegree() for e in elements)
        value_shape = elements[0].value_shape()

        # Sanity check
        assert all(e.get_nodal_basis().get_reference_element() ==
                   ref_el for e in elements)
        assert all(type(e.get_nodal_basis().get_expansion_set()) ==
                   type(expansion_set) for e in elements)
        assert all(e_mapping == mapping for e in elements
                   for e_mapping in e.mapping())
        assert all(e.value_shape() == value_shape for e in elements)

        # Merge polynomial sets
        coeffs = _merge_coeffs([e.get_coeffs() for e in elements])
        dmats = _merge_dmats([e.dmats() for e in elements])
        poly_set = PolynomialSet(ref_el,
                                 degree,
                                 embedded_degree,
                                 expansion_set,
                                 coeffs,
                                 dmats)

        # Renumber dof numbers
        offsets = np.cumsum([0] + [e.space_dimension() for e in elements[:-1]])
        entity_ids = _merge_entity_ids((e.entity_dofs() for e in elements),
                                       offsets)

        # Merge dual bases
        nodes = [node for e in elements for node in e.dual_basis()]
        dual_set = DualSet(nodes, ref_el, entity_ids)

        # CiarletElement constructor adjusts poly_set coefficients s.t.
        # dual_set is really dual to poly_set
        super(NodalEnrichedElement, self).__init__(poly_set, dual_set, order,
                                                   formdegree=formdegree, mapping=mapping)
import random
import numpy as np
import pytest
    import os
from FIAT.reference_element import LINE, ReferenceElement
from FIAT.reference_element import Point, UFCInterval, UFCTriangle, UFCTetrahedron
from FIAT.lagrange import Lagrange
from FIAT.discontinuous_lagrange import DiscontinuousLagrange   # noqa: F401
from FIAT.discontinuous_taylor import DiscontinuousTaylor       # noqa: F401
from FIAT.P0 import P0                                          # noqa: F401
from FIAT.crouzeix_raviart import CrouzeixRaviart               # noqa: F401
from FIAT.raviart_thomas import RaviartThomas                   # noqa: F401
from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas  # noqa: F401
from FIAT.brezzi_douglas_marini import BrezziDouglasMarini      # noqa: F401
from FIAT.mixed import MixedElement
from FIAT.nedelec import Nedelec                                # noqa: F401
from FIAT.nedelec_second_kind import NedelecSecondKind          # noqa: F401
from FIAT.regge import Regge                                    # noqa: F401
from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet   # noqa: F401
from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson  # noqa: F401
from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini  # noqa: F401
from FIAT.gauss_legendre import GaussLegendre                   # noqa: F401
from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre    # noqa: F401
from FIAT.restricted import RestrictedElement                   # noqa: F401
from FIAT.tensor_product import TensorProductElement            # noqa: F401
from FIAT.tensor_product import FlattenedDimensions             # noqa: F401
from FIAT.hdivcurl import Hdiv, Hcurl                           # noqa: F401
from FIAT.argyris import Argyris, QuinticArgyris                # noqa: F401
from FIAT.hermite import CubicHermite                           # noqa: F401
from FIAT.morley import Morley                                  # noqa: F401
from FIAT.bubble import Bubble
from FIAT.enriched import EnrichedElement                       # noqa: F401
from FIAT.nodal_enriched import NodalEnrichedElement
    xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"),
    xfail_impl("HDivTrace(T, 1)"),
    xfail_impl("EnrichedElement("
               "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
               "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
               ")"),
    xfail_impl("EnrichedElement("
               "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
               "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
               ")"),
    # Following elements are checked using tabulate
    xfail_impl("HDivTrace(T, 0)"),
    xfail_impl("HDivTrace(T, 1)"),
    xfail_impl("HDivTrace(T, 2)"),
    xfail_impl("HDivTrace(T, 3)"),
    xfail_impl("HDivTrace(S, 0)"),
    xfail_impl("HDivTrace(S, 1)"),
    xfail_impl("HDivTrace(S, 2)"),
    xfail_impl("HDivTrace(S, 3)"),
    xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"),
    xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"),
    xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"),
    xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))"),
    xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))"),
    xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))"),
    xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))"),
    xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))"),
]


@pytest.mark.parametrize('element', elements)
def test_nodality(element):
    """Check that generated elements are nodal, i.e. nodes evaluated
    on basis functions give Kronecker delta
    """
    # Instantiate element lazily
    element = eval(element)

    # Fetch primal and dual basis
    poly_set = element.get_nodal_basis()
    dual_set = element.get_dual_set()
    assert poly_set.get_reference_element() == dual_set.get_reference_element()

    # Get coeffs of primal and dual bases w.r.t. expansion set
    coeffs_poly = poly_set.get_coeffs()
    coeffs_dual = dual_set.to_riesz(poly_set)
    assert coeffs_poly.shape == coeffs_dual.shape

    # Check nodality
    for i in range(coeffs_dual.shape[0]):
        for j in range(coeffs_poly.shape[0]):
            assert np.isclose(
                coeffs_dual[i].flatten().dot(coeffs_poly[j].flatten()),
                1.0 if i == j else 0.0
            )


@pytest.mark.parametrize('elements', [
    (Lagrange(I, 2), Bubble(I, 2)),
    (Lagrange(T, 3), Bubble(T, 3)),
    (Lagrange(S, 4), Bubble(S, 4)),
    (Lagrange(I, 1), Lagrange(I, 1)),
    (Lagrange(I, 1), Bubble(I, 2), Bubble(I, 2)),
])
def test_illposed_nodal_enriched(elements):
    """Check that nodal enriched element fails on ill-posed
    (non-unisolvent) case
    """
    with pytest.raises(np.linalg.LinAlgError):
        NodalEnrichedElement(*elements)


def test_empty_bubble():
    "Check that bubble of too low degree fails"
    with pytest.raises(RuntimeError):
        Bubble(I, 1)
    with pytest.raises(RuntimeError):
        Bubble(T, 2)
    with pytest.raises(RuntimeError):
        Bubble(S, 3)


def test_nodal_enriched_implementation():
    """Following element pair should be the same.
    This might be fragile to dof reordering but works now.
    """

    e0 = RaviartThomas(T, 2)

    e1 = NodalEnrichedElement(
        RestrictedElement(RaviartThomas(T, 2), restriction_domain='facet'),
        RestrictedElement(RaviartThomas(T, 2), restriction_domain='interior')
    )

    for attr in ["degree",
                 "get_reference_element",
                 "entity_dofs",
                 "entity_closure_dofs",
                 "get_formdegree",
                 "mapping",
                 "num_sub_elements",
                 "space_dimension",
                 "value_shape",
                 "is_nodal",
                 ]:
        assert getattr(e0, attr)() == getattr(e1, attr)()
    assert np.allclose(e0.get_coeffs(), e1.get_coeffs())
    assert np.allclose(e0.dmats(), e1.dmats())
    assert np.allclose(e0.get_dual_set().to_riesz(e0.get_nodal_basis()),
                       e1.get_dual_set().to_riesz(e1.get_nodal_basis()))


def test_mixed_is_nodal():
    element = MixedElement([DiscontinuousLagrange(T, 1), RaviartThomas(T, 2)])

    assert element.is_nodal()


def test_mixed_is_not_nodal():
    element = MixedElement([

        EnrichedElement(