# Path: finat/fiat_elements.py
# class ScalarFiatElement(FiatElement):
#     @property
#     def value_shape(self):
#         return ()
# 
# Path: finat/physically_mapped.py
# class NeedsCoordinateMappingElement(metaclass=ABCMeta):
# class PhysicallyMappedElement(NeedsCoordinateMappingElement):
# class DirectlyDefinedElement(NeedsCoordinateMappingElement):
# class PhysicalGeometry(metaclass=ABCMeta):
#     def __init__(self, *args, **kwargs):
#     def basis_transformation(self, coordinate_mapping):
#     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
#         def matvec(table):
#     def point_evaluation(self, order, refcoords, entity=None):
#     def cell_size(self):
#     def jacobian_at(self, point):
#     def detJ_at(self, point):
#     def reference_normals(self):
#     def physical_normals(self):
#     def physical_tangents(self):
#     def physical_edge_lengths(self):
#     def physical_points(self, point_set, entity=None):
#     def physical_vertices(self):
#         M = self.basis_transformation(coordinate_mapping)


import numpy
import FIAT
from gem import Literal, ListTensor
from finat.fiat_elements import ScalarFiatElement
from finat.physically_mapped import PhysicallyMappedElement, Citations

        V = numpy.zeros((21, 18), dtype=object)

        for multiindex in numpy.ndindex(V.shape):
            V[multiindex] = Literal(V[multiindex])

        for v in range(3):
            s = 6*v
            V[s, s] = Literal(1)
            for i in range(2):
                for j in range(2):
                    V[s+1+i, s+1+j] = J[j, i]
            V[s+3, s+3] = J[0, 0]*J[0, 0]
            V[s+3, s+4] = 2*J[0, 0]*J[1, 0]
            V[s+3, s+5] = J[1, 0]*J[1, 0]
            V[s+4, s+3] = J[0, 0]*J[0, 1]
            V[s+4, s+4] = J[0, 0]*J[1, 1] + J[1, 0]*J[0, 1]
            V[s+4, s+5] = J[1, 0]*J[1, 1]
            V[s+5, s+3] = J[0, 1]*J[0, 1]
            V[s+5, s+4] = 2*J[0, 1]*J[1, 1]
            V[s+5, s+5] = J[1, 1]*J[1, 1]

        for e in range(3):
            v0id, v1id = [i for i in range(3) if i != e]

            # nhat . J^{-T} . t
            foo = (rns[e, 0]*(J[0, 0]*pts[e, 0] + J[1, 0]*pts[e, 1])
                   + rns[e, 1]*(J[0, 1]*pts[e, 0] + J[1, 1]*pts[e, 1]))

            # vertex points

# Path: finat/fiat_elements.py
# class ScalarFiatElement(FiatElement):
#     @property
#     def value_shape(self):
#         return ()
# 
# Path: finat/physically_mapped.py
# class NeedsCoordinateMappingElement(metaclass=ABCMeta):
# class PhysicallyMappedElement(NeedsCoordinateMappingElement):
# class DirectlyDefinedElement(NeedsCoordinateMappingElement):
# class PhysicalGeometry(metaclass=ABCMeta):
#     def __init__(self, *args, **kwargs):
#     def basis_transformation(self, coordinate_mapping):
#     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
#         def matvec(table):
#     def point_evaluation(self, order, refcoords, entity=None):
#     def cell_size(self):
#     def jacobian_at(self, point):
#     def detJ_at(self, point):
#     def reference_normals(self):
#     def physical_normals(self):
#     def physical_tangents(self):
#     def physical_edge_lengths(self):
#     def physical_points(self, point_set, entity=None):
#     def physical_vertices(self):
#         M = self.basis_transformation(coordinate_mapping)


# Path: finat/bell.py
import numpy
import FIAT
from gem import Literal, ListTensor
from finat.fiat_elements import ScalarFiatElement
from finat.physically_mapped import PhysicallyMappedElement, Citations


        V = numpy.zeros((21, 18), dtype=object)

        for multiindex in numpy.ndindex(V.shape):
            V[multiindex] = Literal(V[multiindex])

        for v in range(3):
            s = 6*v
            V[s, s] = Literal(1)
            for i in range(2):
                for j in range(2):
                    V[s+1+i, s+1+j] = J[j, i]
            V[s+3, s+3] = J[0, 0]*J[0, 0]
            V[s+3, s+4] = 2*J[0, 0]*J[1, 0]
            V[s+3, s+5] = J[1, 0]*J[1, 0]
            V[s+4, s+3] = J[0, 0]*J[0, 1]
            V[s+4, s+4] = J[0, 0]*J[1, 1] + J[1, 0]*J[0, 1]
            V[s+4, s+5] = J[1, 0]*J[1, 1]
            V[s+5, s+3] = J[0, 1]*J[0, 1]
            V[s+5, s+4] = 2*J[0, 1]*J[1, 1]
            V[s+5, s+5] = J[1, 1]*J[1, 1]

        for e in range(3):
            v0id, v1id = [i for i in range(3) if i != e]

            # nhat . J^{-T} . t
            foo = (rns[e, 0]*(J[0, 0]*pts[e, 0] + J[1, 0]*pts[e, 1])
                   + rns[e, 1]*(J[0, 1]*pts[e, 0] + J[1, 1]*pts[e, 1]))

            # vertex points

            V[18+e, 6*v0id] = -1/21 * (foo / pel[e])