# Path: FIAT/expansions.py
# def jrc(a, b, n):
# def _tabulate_dpts(tabulator, D, n, order, pts):
#     def form_derivative(F):
#     def numpy_lambdify(X, F):
#     def evaluate_lambda(lmbd, x):
# def xi_triangle(eta):
# def xi_tetrahedron(eta):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def tabulate_derivatives(self, n, pts):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def tabulate_derivatives(self, n, pts):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def _tabulate(self, n, pts):
#         def idx(p, q):
#     def tabulate_derivatives(self, n, pts):
#     def tabulate_jet(self, n, pts, order=1):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def _tabulate(self, n, pts):
#         def idx(p, q, r):
#     def tabulate_derivatives(self, n, pts):
#     def tabulate_jet(self, n, pts, order=1):
# def get_expansion_set(ref_el):
# def polynomial_dimension(ref_el, degree):
#     X = sympy.DeferredVector('x')
#         D = 3
# class PointExpansionSet(object):
# class LineExpansionSet(object):
# class TriangleExpansionSet(object):
# class TetrahedronExpansionSet(object):
# 
# Path: FIAT/functional.py
# def index_iterator(shp):
#     """Constructs a generator iterating over all indices in
#     shp in generalized column-major order  So if shp = (2,2), then we
#     construct the sequence (0,0),(0,1),(1,0),(1,1)"""
#     if len(shp) == 0:
#         return
#     elif len(shp) == 1:
#         for i in range(shp[0]):
#             yield [i]
#     else:
#         shp_foo = shp[1:]
#         for i in range(shp[0]):
#             for foo in index_iterator(shp_foo):
#                 yield [i] + foo


import numpy
from FIAT import expansions
from FIAT.functional import index_iterator
                         A.get_embedded_degree(),
                         A.get_expansion_set(),
                         coeffs,
                         A.get_dmats())


class ONSymTensorPolynomialSet(PolynomialSet):
    """Constructs an orthonormal basis for symmetric-tensor-valued
    polynomials on a reference element.

    """

    def __init__(self, ref_el, degree, size=None):

        sd = ref_el.get_spatial_dimension()
        if size is None:
            size = sd

        shape = (size, size)
        num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
        num_components = size * (size + 1) // 2
        num_members = num_components * num_exp_functions
        embedded_degree = degree
        expansion_set = expansions.get_expansion_set(ref_el)

        # set up coefficients for symmetric tensors
        coeffs_shape = tuple([num_members] + list(shape) + [num_exp_functions])
        coeffs = numpy.zeros(coeffs_shape, "d")
        cur_bf = 0
        for [i, j] in index_iterator(shape):

# Path: FIAT/expansions.py
# def jrc(a, b, n):
# def _tabulate_dpts(tabulator, D, n, order, pts):
#     def form_derivative(F):
#     def numpy_lambdify(X, F):
#     def evaluate_lambda(lmbd, x):
# def xi_triangle(eta):
# def xi_tetrahedron(eta):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def tabulate_derivatives(self, n, pts):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def tabulate_derivatives(self, n, pts):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def _tabulate(self, n, pts):
#         def idx(p, q):
#     def tabulate_derivatives(self, n, pts):
#     def tabulate_jet(self, n, pts, order=1):
#     def __init__(self, ref_el):
#     def get_num_members(self, n):
#     def tabulate(self, n, pts):
#     def _tabulate(self, n, pts):
#         def idx(p, q, r):
#     def tabulate_derivatives(self, n, pts):
#     def tabulate_jet(self, n, pts, order=1):
# def get_expansion_set(ref_el):
# def polynomial_dimension(ref_el, degree):
#     X = sympy.DeferredVector('x')
#         D = 3
# class PointExpansionSet(object):
# class LineExpansionSet(object):
# class TriangleExpansionSet(object):
# class TetrahedronExpansionSet(object):
# 
# Path: FIAT/functional.py
# def index_iterator(shp):
#     """Constructs a generator iterating over all indices in
#     shp in generalized column-major order  So if shp = (2,2), then we
#     construct the sequence (0,0),(0,1),(1,0),(1,1)"""
#     if len(shp) == 0:
#         return
#     elif len(shp) == 1:
#         for i in range(shp[0]):
#             yield [i]
#     else:
#         shp_foo = shp[1:]
#         for i in range(shp[0]):
#             for foo in index_iterator(shp_foo):
#                 yield [i] + foo


# Path: FIAT/polynomial_set.py
import numpy
from FIAT import expansions
from FIAT.functional import index_iterator

                         A.get_embedded_degree(),
                         A.get_expansion_set(),
                         coeffs,
                         A.get_dmats())


class ONSymTensorPolynomialSet(PolynomialSet):
    """Constructs an orthonormal basis for symmetric-tensor-valued
    polynomials on a reference element.

    """

    def __init__(self, ref_el, degree, size=None):

        sd = ref_el.get_spatial_dimension()
        if size is None:
            size = sd

        shape = (size, size)
        num_exp_functions = expansions.polynomial_dimension(ref_el, degree)
        num_components = size * (size + 1) // 2
        num_members = num_components * num_exp_functions
        embedded_degree = degree
        expansion_set = expansions.get_expansion_set(ref_el)

        # set up coefficients for symmetric tensors
        coeffs_shape = tuple([num_members] + list(shape) + [num_exp_functions])
        coeffs = numpy.zeros(coeffs_shape, "d")
        cur_bf = 0
        for [i, j] in index_iterator(shape):

            n = expansions.polynomial_dimension(ref_el, embedded_degree)