class FiniteElement(object):
class CiarletElement(FiniteElement):
    def __init__(self, ref_el, dual, order, formdegree=None, mapping="affine"):
    def get_reference_element(self):
    def get_dual_set(self):
    def get_order(self):
    def dual_basis(self):
    def entity_dofs(self):
    def entity_closure_dofs(self):
    def get_formdegree(self):
    def mapping(self):
    def num_sub_elements(self):
    def space_dimension(self):
    def tabulate(self, order, points, entity=None):
    def is_nodal():
    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_el=None):
    def degree(self):
    def get_nodal_basis(self):
    def get_coeffs(self):
    def tabulate(self, order, points, entity=None):
    def value_shape(self):
    def dmats(self):
    def get_num_members(self, arg):
    def is_nodal():
def entity_support_dofs(elem, entity_dim):
            A = numpy.reshape(dualmat, (dualmat.shape[0], num_cols))
            B = numpy.reshape(old_coeffs, (old_coeffs.shape[0], num_cols))
            A = dualmat
            B = old_coeffs
        V = numpy.dot(A, numpy.transpose(B)) class DualSet(object):
    def __init__(self, nodes, ref_el, entity_ids):
    def get_nodes(self):
    def get_entity_closure_ids(self):
    def get_entity_ids(self):
    def get_reference_element(self):
    def to_riesz(self, poly_set):
def make_entity_closure_ids(ref_el, entity_ids): def index_iterator(shp):
    def __init__(self, ref_el, target_shape, pt_dict, deriv_dict,
                 functional_type):
    def evaluate(self, f):
    def __call__(self, fn):
    def get_point_dict(self):
    def get_reference_element(self):
    def get_type_tag(self):
    def to_riesz(self, poly_set):
    def tostr(self):
    def __init__(self, ref_el, x):
    def __call__(self, fn):
    def tostr(self):
    def __init__(self, ref_el, comp, shp, x):
    def tostr(self):
    def __init__(self, ref_el, x, alpha):
    def __call__(self, fn):
    def __init__(self, ref_el, facet_no, pt):
    def __init__(self, ref_el, facet_no, pt):
    def __init__(self, ref_el, Q, f_at_qpts, comp=tuple(), shp=tuple()):
    def __call__(self, fn):
    def __init__(self, ref_el, facet_no, Q, f_at_qpts):
    def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""):
    def __init__(self, cell, entity, mom_deg, comp_deg):
    def __init__(self, cell, entity, mom_deg, comp_deg):
    def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""):
    def __init__(self, cell, entity, mom_deg, comp_deg):
    def __init__(self, cell, entity, mom_deg, comp_deg):
    def __init__(self, ref_el, Q, f_at_qpts):
    def __init__(self, ref_el, Q, f_at_qpts):
    def __init__(self, ref_el, Q, f_at_qpts):
    def __init__(self, ref_el, facet_no, pt):
    def __init__(self, ref_el, edge_no, pt):
    def tostr(self):
    def __init__(self, ref_el, Q, P_at_qpts, edge):
    def __init__(self, ref_el, face_no, tno, pt):
    def tostr(self):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
    def __init__(self, ref_el, facet_no, pt):
    def tostr(self):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
    def __init__(self, ref_el, v, w, p):
    def __init__(self, ref_el, v, w, Q, f_at_qpts, comp_deg):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
    def __init__(self, ref_el, Q, P_at_qpts, facet):
class Functional(object):
class PointEvaluation(Functional):
class ComponentPointEvaluation(Functional):
class PointDerivative(Functional):
class PointNormalDerivative(Functional):
class PointNormalSecondDerivative(Functional):
class IntegralMoment(Functional):
class IntegralMomentOfNormalDerivative(Functional):
class IntegralLegendreDirectionalMoment(Functional):
class IntegralLegendreNormalMoment(IntegralLegendreDirectionalMoment):
class IntegralLegendreTangentialMoment(IntegralLegendreDirectionalMoment):
class IntegralLegendreBidirectionalMoment(Functional):
class IntegralLegendreNormalNormalMoment(IntegralLegendreBidirectionalMoment):
class IntegralLegendreNormalTangentialMoment(IntegralLegendreBidirectionalMoment):
class IntegralMomentOfDivergence(Functional):
class IntegralMomentOfTensorDivergence(Functional):
class FrobeniusIntegralMoment(Functional):
class PointNormalEvaluation(Functional):
class PointEdgeTangentEvaluation(Functional):
class IntegralMomentOfEdgeTangentEvaluation(Functional):
class PointFaceTangentEvaluation(Functional):
class IntegralMomentOfFaceTangentEvaluation(Functional):
class MonkIntegralMoment(Functional):
class PointScaledNormalEvaluation(Functional):
class IntegralMomentOfScaledNormalEvaluation(Functional):
class PointwiseInnerProductEvaluation(Functional):
class TensorBidirectionalMomentInnerProductEvaluation(Functional):
class IntegralMomentOfNormalEvaluation(Functional):
class IntegralMomentOfTangentialEvaluation(Functional):
class IntegralMomentOfNormalNormalEvaluation(Functional):
        X = sympy.DeferredVector('x')
        Q = GaussLegendreQuadratureLineRule(interval(), quadpoints)
        Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) 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 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) 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 FacetBubble(CodimBubble):
    """The facet bubble finite element: the dofs of the Lagrange FE in the interior of the facets"""

    def __init__(self, ref_el, degree):
        super(FacetBubble, self).__init__(ref_el, degree, codim=1) 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) 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]) def create_quadrature(ref_el, degree, scheme="default"):
    """
    Generate quadrature rule for given reference element
    that will integrate an polynomial of order 'degree' exactly.

    For low-degree (<=6) polynomials on triangles and tetrahedra, this
    uses hard-coded rules, otherwise it falls back to a collapsed
    Gauss scheme on simplices.  On tensor-product cells, it is a
    tensor-product quadrature rule of the subcells.

    :arg cell: The FIAT cell to create the quadrature for.
    :arg degree: The degree of polynomial that the rule should
        integrate exactly.
    """
    if ref_el.get_shape() == TENSORPRODUCT:
        try:
            degree = tuple(degree)
        except TypeError:
            degree = (degree,) * len(ref_el.cells)

        assert len(ref_el.cells) == len(degree)
        quad_rules = [create_quadrature(c, d, scheme)
                      for c, d in zip(ref_el.cells, degree)]
        return make_tensor_product_quadrature(*quad_rules)

    if ref_el.get_shape() in [QUADRILATERAL, HEXAHEDRON]:
        return create_quadrature(ref_el.product, degree, scheme)

    if degree < 0:
        raise ValueError("Need positive degree, not %d" % degree)

    if scheme == "default":
        # TODO: Point transformation to support good schemes on
        # non-UFC reference elements.
        if isinstance(ref_el, UFCTriangle):
            return _triangle_scheme(degree)
        elif isinstance(ref_el, UFCTetrahedron):
            return _tetrahedron_scheme(degree)
        else:
            return _fiat_scheme(ref_el, degree)
    elif scheme == "canonical":
        return _fiat_scheme(ref_el, degree)
    elif scheme == "KMV":  # Kong-Mulder-Veldhuizen scheme
        return _kmv_lump_scheme(ref_el, degree)
    else:
        raise ValueError("Unknown quadrature scheme: %s." % scheme)
from FIAT import (
    finite_element,
    dual_set,
    functional,
    Bubble,
    FacetBubble,
    Lagrange,
    NodalEnrichedElement,
    RestrictedElement,
    reference_element,
)
from FIAT.quadrature_schemes import create_quadrature
# Copyright (C) 2020 Robert C. Kirby (Baylor University)
#
# contributions by Keith Roberts (University of São Paulo)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later

TRIANGLE = reference_element.UFCTriangle()
TETRAHEDRON = reference_element.UFCTetrahedron()


def _get_entity_ids(ref_el, degree):
    """The topological association in a dictionary"""
    T = ref_el.topology
    sd = ref_el.get_spatial_dimension()
    if degree == 1:  # works for any spatial dimension.
        entity_ids = {0: dict((i, [i]) for i in range(len(T[0])))}
        for d in range(1, sd + 1):
            entity_ids[d] = dict((i, []) for i in range(len(T[d])))
    elif degree == 2:
        if sd == 2:
            entity_ids = {
                0: dict((i, [i]) for i in range(3)),
                1: dict((i, [i + 3]) for i in range(3)),
                2: {0: [6]},
            }
        elif sd == 3:
            entity_ids = {
                0: dict((i, [i]) for i in range(4)),
                1: dict((i, [i + 4]) for i in range(6)),
                2: dict((i, [i + 10]) for i in range(4)),
                3: {0: [14]},
            }
    elif degree == 3:
        if sd == 2:
            etop = [[3, 4], [6, 5], [7, 8]]
            entity_ids = {
                0: dict((i, [i]) for i in range(3)),
                1: dict((i, etop[i]) for i in range(3)),
                2: {0: [9, 10, 11]},
            }
        elif sd == 3:
            etop = [[4, 5], [7, 6], [8, 9], [11, 10], [12, 13], [14, 15]]
            ftop = [[16, 17, 18], [19, 20, 21], [22, 23, 24], [25, 26, 27]]
            entity_ids = {
                0: dict((i, [i]) for i in range(4)),
                1: dict((i, etop[i]) for i in range(6)),
                2: dict((i, ftop[i]) for i in range(4)),
                3: {0: [28, 29, 30, 31]},
            }
    elif degree == 4:
        if sd == 2:
            etop = [[6, 3, 7], [9, 4, 8], [10, 5, 11]]
            entity_ids = {
                0: dict((i, [i]) for i in range(3)),
                1: dict((i, etop[i]) for i in range(3)),
                2: {0: [i for i in range(12, 18)]},
            }
    elif degree == 5:
        if sd == 2:
            etop = [[9, 3, 4, 10], [12, 6, 5, 11], [13, 7, 8, 14]]
            entity_ids = {
                0: dict((i, [i]) for i in range(3)),
                1: dict((i, etop[i]) for i in range(3)),
                2: {0: [i for i in range(15, 30)]},
            }
    return entity_ids


def bump(T, deg):
    """Increase degree of polynomial along face/edges"""
    sd = T.get_spatial_dimension()
    if deg == 1:
        return (0, 0)
    else:
        if sd == 2:
            if deg < 5:
                return (1, 1)
            elif deg == 5:
                return (2, 2)
            else:
                raise ValueError("Degree not supported")
        elif sd == 3:
            if deg < 4:
                return (1, 2)
            else:
                raise ValueError("Degree not supported")
        else:
            raise ValueError("Dimension of element is not supported")


def KongMulderVeldhuizenSpace(T, deg):
    sd = T.get_spatial_dimension()
    if deg == 1:
        return Lagrange(T, 1).poly_set
    else:

        L = Lagrange(T, deg)