class CiarletElement(FiniteElement):
    """Class implementing Ciarlet's abstraction of a finite element
    being a domain, function space, and set of nodes.

    Elements derived from this class are nodal finite elements, with a nodal
    basis generated from polynomials encoded in a `PolynomialSet`.
    """

    def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_el=None):
        ref_el = ref_el or poly_set.get_reference_element()
        super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping)

        # build generalized Vandermonde matrix
        old_coeffs = poly_set.get_coeffs()
        dualmat = dual.to_riesz(poly_set)

        shp = dualmat.shape
        if len(shp) > 2:
            num_cols = numpy.prod(shp[1:])

            A = numpy.reshape(dualmat, (dualmat.shape[0], num_cols))
            B = numpy.reshape(old_coeffs, (old_coeffs.shape[0], num_cols))
        else:
            A = dualmat
            B = old_coeffs

        V = numpy.dot(A, numpy.transpose(B))
        self.V = V

        Vinv = numpy.linalg.inv(V)

        new_coeffs_flat = numpy.dot(numpy.transpose(Vinv), B)

        new_shp = tuple([new_coeffs_flat.shape[0]] + list(shp[1:]))
        new_coeffs = numpy.reshape(new_coeffs_flat, new_shp)

        self.poly_set = PolynomialSet(ref_el,
                                      poly_set.get_degree(),
                                      poly_set.get_embedded_degree(),
                                      poly_set.get_expansion_set(),
                                      new_coeffs,
                                      poly_set.get_dmats())

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

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

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

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

        :arg order: The maximum order of derivative.
        :arg points: An iterable of points.
        :arg entity: Optional (dimension, entity number) pair
                     indicating which topological entity of the
                     reference element to tabulate on.  If ``None``,
                     default cell-wise tabulation is performed.
        """
        if entity is None:
            entity = (self.ref_el.get_spatial_dimension(), 0)

        entity_dim, entity_id = entity
        transform = self.ref_el.get_entity_transform(entity_dim, entity_id)
        return self.poly_set.tabulate(list(map(transform, points)), order)

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

    def dmats(self):
        """Return dmats: expansion coefficients for basis function
        derivatives."""
        return self.get_nodal_basis().get_dmats()

    def get_num_members(self, arg):
        "Return number of members of the expansion set."
        return self.get_nodal_basis().get_expansion_set().get_num_members(arg)

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

        All implementations/subclasses are nodal including this one.
        """
        return True class DualSet(object):
    def __init__(self, nodes, ref_el, entity_ids):
        self.nodes = nodes
        self.ref_el = ref_el
        self.entity_ids = entity_ids

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

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

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

    def get_nodes(self):
        return self.nodes

    def get_entity_closure_ids(self):
        return self.entity_closure_ids

    def get_entity_ids(self):
        return self.entity_ids

    def get_reference_element(self):
        return self.ref_el

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

        return self.mat class ONPolynomialSet(PolynomialSet):
    """Constructs an orthonormal basis out of expansion set by having an
    identity matrix of coefficients.  Can be used to specify ON bases
    for vector- and tensor-valued sets as well.

    """

    def __init__(self, ref_el, degree, shape=tuple()):

        if shape == tuple():
            num_components = 1
        else:
            flat_shape = numpy.ravel(shape)
            num_components = numpy.prod(flat_shape)
        num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
        num_members = num_components * num_exp_functions
        embedded_degree = degree
        expansion_set = expansions.get_expansion_set(ref_el)
        sd = ref_el.get_spatial_dimension()

        # set up coefficients
        coeffs_shape = tuple([num_members] + list(shape) + [num_exp_functions])
        coeffs = numpy.zeros(coeffs_shape, "d")

        # use functional's index_iterator function
        cur_bf = 0

        if shape == tuple():
            coeffs = numpy.eye(num_members)
        else:
            for idx in index_iterator(shape):
                n = expansions.polynomial_dimension(ref_el, embedded_degree)
                for exp_bf in range(n):
                    cur_idx = tuple([cur_bf] + list(idx) + [exp_bf])
                    coeffs[cur_idx] = 1.0
                    cur_bf += 1

        # construct dmats
        if degree == 0:
            dmats = [numpy.array([[0.0]], "d") for i in range(sd)]
        else:
            pts = ref_el.make_points(sd, 0, degree + sd + 1)

            v = numpy.transpose(expansion_set.tabulate(degree, pts))
            vinv = numpy.linalg.inv(v)

            dv = expansion_set.tabulate_derivatives(degree, pts)
            dtildes = [[[a[1][i] for a in dvrow] for dvrow in dv]
                       for i in range(sd)]

            dmats = [numpy.dot(vinv, numpy.transpose(dtilde))
                     for dtilde in dtildes]

        PolynomialSet.__init__(self, ref_el, degree, embedded_degree,
                               expansion_set, coeffs, dmats) class PointEdgeTangentEvaluation(Functional):
    """Implements the evaluation of the tangential component of a
    vector at a point on a facet of dimension 1."""

    def __init__(self, ref_el, edge_no, pt):
        t = ref_el.compute_edge_tangent(edge_no)
        self.t = t
        sd = ref_el.get_spatial_dimension()
        pt_dict = {pt: [(t[i], (i,)) for i in range(sd)]}
        shp = (sd,)
        super().__init__(ref_el, shp, pt_dict, {}, "PointEdgeTangent")

    def tostr(self):
        x = list(map(str, list(self.pt_dict.keys())[0]))
        return "(u.t)(%s)" % (','.join(x),) class FrobeniusIntegralMoment(Functional):

    def __init__(self, ref_el, Q, f_at_qpts):
        # f_at_qpts is (some shape) x num_qpts
        shp = tuple(f_at_qpts.shape[:-1])
        if len(Q.get_points()) != f_at_qpts.shape[-1]:
            raise Exception("Mismatch in number of quadrature points and values")

        qpts, qwts = Q.get_points(), Q.get_weights()
        pt_dict = {}

        for i, (pt_cur, wt_cur) in enumerate(zip(map(tuple, qpts), qwts)):
            pt_dict[pt_cur] = []
            for alfa in index_iterator(shp):
                qpidx = tuple(alfa + [i])
                pt_dict[pt_cur].append((wt_cur * f_at_qpts[qpidx], tuple(alfa)))

        super().__init__(ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment") 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") def make_quadrature(ref_el, m):
    """Returns the collapsed quadrature rule using m points per
    direction on the given reference element. In the tensor product
    case, m is a tuple."""

    if isinstance(m, tuple):
        min_m = min(m)
    else:
        min_m = m

    msg = "Expecting at least one (not %d) quadrature point per direction" % min_m
    assert (min_m > 0), msg

    if ref_el.get_shape() == reference_element.POINT:
        return QuadratureRule(ref_el, [()], [1])
    elif ref_el.get_shape() == reference_element.LINE:
        return GaussJacobiQuadratureLineRule(ref_el, m)
    elif ref_el.get_shape() == reference_element.TRIANGLE:
        return CollapsedQuadratureTriangleRule(ref_el, m)
    elif ref_el.get_shape() == reference_element.TETRAHEDRON:
        return CollapsedQuadratureTetrahedronRule(ref_el, m)
    elif ref_el.get_shape() == reference_element.QUADRILATERAL:
        line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
        return make_tensor_product_quadrature(line_rule, line_rule)
    elif ref_el.get_shape() == reference_element.HEXAHEDRON:
        line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m)
        return make_tensor_product_quadrature(line_rule, line_rule, line_rule)
    else:
        raise ValueError("Unable to make quadrature for cell: %s" % ref_el) class UFCTetrahedronFaceQuadratureRule(QuadratureRule):
    """Highly specialized quadrature rule for the face of a
    tetrahedron, mapped from a reference triangle, used for higher
    order Nedelecs"""

    def __init__(self, face_number, degree):

        # Create quadrature rule on reference triangle
        reference_triangle = reference_element.UFCTriangle()
        reference_rule = make_quadrature(reference_triangle, degree)
        ref_points = reference_rule.get_points()
        ref_weights = reference_rule.get_weights()

        # Get geometry information about the face of interest
        reference_tet = reference_element.UFCTetrahedron()
        face = reference_tet.get_topology()[2][face_number]
        vertices = reference_tet.get_vertices_of_subcomplex(face)

        # Use tet to map points and weights on the appropriate face
        vertices = [numpy.array(list(vertex)) for vertex in vertices]
        x0 = vertices[0]
        J = numpy.vstack([vertices[1] - x0, vertices[2] - x0]).T
        # This is just a very numpyfied way of writing J*p + x0:
        points = numpy.einsum("ij,kj->ki", J, ref_points) + x0

        # Map weights: multiply reference weights by sqrt(|J^T J|)
        detJTJ = numpy.linalg.det(numpy.dot(J.T, J))
        weights = numpy.sqrt(detJTJ) * ref_weights

        # Initialize super class with new points and weights
        QuadratureRule.__init__(self, reference_tet, points, weights)
        self._reference_rule = reference_rule
        self._J = J

    def reference_rule(self):
        return self._reference_rule

    def jacobian(self):
        return self._J 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) def check_format_variant(variant, degree, element):
    if variant is None:
        variant = "point"
        warnings.simplefilter('always', DeprecationWarning)
        warnings.warn('Variant of ' + element + ' element will change from point evaluation to integral evaluation.'
                      ' You should project into variant="integral"', DeprecationWarning)

    match = re.match(r"^integral(?:\((\d+)\))?$", variant)
    if match:
        variant = "integral"
        quad_degree, = match.groups()
        quad_degree = int(quad_degree) if quad_degree is not None else 5*(degree + 1)
        if quad_degree < degree + 1:
            raise ValueError("Warning, quadrature degree should be at least %s" % (degree + 1))
    elif variant == "point":
        quad_degree = None
    else:
        raise ValueError('Choose either variant="point" or variant="integral"'
                         'or variant="integral(Quadrature degree)"')

    return (variant, quad_degree) def mis(m, n):
    def __init__(self, ref_el, degree, embedded_degree, expansion_set, coeffs,
                 dmats):
    def tabulate_new(self, pts):
    def tabulate(self, pts, jet_order=0):
    def get_expansion_set(self):
    def get_coeffs(self):
    def get_num_members(self):
    def get_degree(self):
    def get_embedded_degree(self):
    def get_dmats(self):
    def get_reference_element(self):
    def get_shape(self):
    def take(self, items):
    def __init__(self, ref_el, degree, shape=tuple()):
def project(f, U, Q):
def form_matrix_product(mats, alpha):
def polynomial_set_union_normalized(A, B):
    def __init__(self, ref_el, degree, size=None):
class PolynomialSet(object):
class ONPolynomialSet(PolynomialSet):
class ONSymTensorPolynomialSet(PolynomialSet):
                    D = form_matrix_product(self.dmats, alpha)
                    D = numpy.eye(1) class QuadratureRule(object):
class GaussJacobiQuadratureLineRule(QuadratureRule):
class GaussLobattoLegendreQuadratureLineRule(QuadratureRule):
class GaussLegendreQuadratureLineRule(QuadratureRule):
class RadauQuadratureLineRule(QuadratureRule):
class CollapsedQuadratureTriangleRule(QuadratureRule):
class CollapsedQuadratureTetrahedronRule(QuadratureRule):
class UFCTetrahedronFaceQuadratureRule(QuadratureRule):
    def __init__(self, ref_el, pts, wts):
    def get_points(self):
    def get_weights(self):
    def integrate(self, f):
    def __init__(self, ref_el, m):
    def __init__(self, ref_el, m):
    def __init__(self, ref_el, m):
    def __init__(self, ref_el, m, right=True):
    def __init__(self, ref_el, m):
    def __init__(self, ref_el, m):
    def __init__(self, face_number, degree):
    def reference_rule(self):
    def jacobian(self):
def make_quadrature(ref_el, m):
def make_tensor_product_quadrature(*quad_rules):
def compute_gauss_jacobi_points(a, b, m):
def compute_gauss_jacobi_rule(a, b, m):
        N = m - 1
        P = numpy.zeros((N + 1, N + 2))
            P[0, :] = (-1) ** numpy.arange(0, N + 2)
        P = P[:, :-1]
        J = numpy.vstack([vertices[1] - x0, vertices[2] - x0]).T 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)
import numpy
from FIAT.finite_element import CiarletElement
from FIAT.dual_set import DualSet
from FIAT.polynomial_set import ONPolynomialSet
from FIAT.functional import PointEdgeTangentEvaluation as Tangent
from FIAT.functional import FrobeniusIntegralMoment as IntegralMoment
from FIAT.raviart_thomas import RaviartThomas
from FIAT.quadrature import make_quadrature, UFCTetrahedronFaceQuadratureRule
from FIAT.reference_element import UFCTetrahedron
from FIAT.check_format_variant import check_format_variant
from FIAT import polynomial_set, quadrature, functional
# Copyright (C) 2010-2012 Marie E. Rognes
# Modified by Andrew T. T. McRae (Imperial College London)
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later





class NedelecSecondKindDual(DualSet):
    r"""
    This class represents the dual basis for the Nedelec H(curl)
    elements of the second kind. The degrees of freedom (L) for the
    elements of the k'th degree are

    d = 2:

      vertices: None

      edges:    L(f) = f (x_i) * t       for (k+1) points x_i on each edge

      cell:     L(f) = \int f * g * dx   for g in RT_{k-1}


    d = 3:

      vertices: None

      edges:    L(f)  = f(x_i) * t         for (k+1) points x_i on each edge

      faces:    L(f) = \int_F f * g * ds   for g in RT_{k-1}(F) for each face F

      cell:     L(f) = \int f * g * dx     for g in RT_{k-2}

    Higher spatial dimensions are not yet implemented. (For d = 1,
    these elements coincide with the CG_k elements.)
    """

    def __init__(self, cell, degree, variant, quad_deg):

        # Define degrees of freedom
        (dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, quad_deg)
        # Call init of super-class
        super(NedelecSecondKindDual, self).__init__(dofs, cell, ids)

    def generate_degrees_of_freedom(self, cell, degree, variant, quad_deg):
        "Generate dofs and geometry-to-dof maps (ids)."

        dofs = []
        ids = {}

        # Extract spatial dimension and topology
        d = cell.get_spatial_dimension()
        assert (d in (2, 3)), "Second kind Nedelecs only implemented in 2/3D."

        # Zero vertex-based degrees of freedom (d+1 of these)
        ids[0] = dict(list(zip(list(range(d + 1)), ([] for i in range(d + 1)))))

        # (d+1) degrees of freedom per entity of codimension 1 (edges)
        (edge_dofs, edge_ids) = self._generate_edge_dofs(cell, degree, 0, variant, quad_deg)
        dofs.extend(edge_dofs)
        ids[1] = edge_ids

        # Include face degrees of freedom if 3D
        if d == 3:
            (face_dofs, face_ids) = self._generate_face_dofs(cell, degree,
                                                             len(dofs), variant, quad_deg)
            dofs.extend(face_dofs)
            ids[2] = face_ids

        # Varying degrees of freedom (possibly zero) per cell
        (cell_dofs, cell_ids) = self._generate_cell_dofs(cell, degree, len(dofs), variant, quad_deg)
        dofs.extend(cell_dofs)
        ids[d] = cell_ids

        return (dofs, ids)

    def _generate_edge_dofs(self, cell, degree, offset, variant, quad_deg):
        """Generate degrees of freedoms (dofs) for entities of
        codimension 1 (edges)."""

        # (degree+1) tangential component point evaluation degrees of
        # freedom per entity of codimension 1 (edges)
        dofs = []
        ids = {}

        if variant == "integral":
            edge = cell.construct_subelement(1)
            Q = quadrature.make_quadrature(edge, quad_deg)

            Pq = polynomial_set.ONPolynomialSet(edge, degree)