# 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
class Morley(PhysicallyMappedElement, ScalarFiatElement):
    def __init__(self, cell, degree):
        if degree != 2:
            raise ValueError("Degree must be 2 for Morley element")
        if Citations is not None:
            Citations().register("Morley1971")
        super().__init__(FIAT.Morley(cell))

    def basis_transformation(self, coordinate_mapping):
        # Jacobians at edge midpoints
        J = coordinate_mapping.jacobian_at([1/3, 1/3])

        rns = coordinate_mapping.reference_normals()
        pns = coordinate_mapping.physical_normals()

        pts = coordinate_mapping.physical_tangents()

        pel = coordinate_mapping.physical_edge_lengths()

        V = numpy.eye(6, dtype=object)
        for multiindex in numpy.ndindex(V.shape):
            V[multiindex] = Literal(V[multiindex])

        for i in range(3):
            V[i+3, i+3] = (rns[i, 0]*(pns[i, 0]*J[0, 0] + pns[i, 1]*J[1, 0])
                           + rns[i, 1]*(pns[i, 0]*J[0, 1] + pns[i, 1]*J[1, 1]))

        for i, c in enumerate([(1, 2), (0, 2), (0, 1)]):
            B12 = (rns[i, 0]*(pts[i, 0]*J[0, 0] + pts[i, 1]*J[1, 0])
                   + rns[i, 1]*(pts[i, 0]*J[0, 1] + pts[i, 1]*J[1, 1]))

# 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/morley.py
import numpy
import FIAT
from gem import Literal, ListTensor
from finat.fiat_elements import ScalarFiatElement
from finat.physically_mapped import PhysicallyMappedElement, Citations

class Morley(PhysicallyMappedElement, ScalarFiatElement):
    def __init__(self, cell, degree):
        if degree != 2:
            raise ValueError("Degree must be 2 for Morley element")
        if Citations is not None:
            Citations().register("Morley1971")
        super().__init__(FIAT.Morley(cell))

    def basis_transformation(self, coordinate_mapping):
        # Jacobians at edge midpoints
        J = coordinate_mapping.jacobian_at([1/3, 1/3])

        rns = coordinate_mapping.reference_normals()
        pns = coordinate_mapping.physical_normals()

        pts = coordinate_mapping.physical_tangents()

        pel = coordinate_mapping.physical_edge_lengths()

        V = numpy.eye(6, dtype=object)
        for multiindex in numpy.ndindex(V.shape):
            V[multiindex] = Literal(V[multiindex])

        for i in range(3):
            V[i+3, i+3] = (rns[i, 0]*(pns[i, 0]*J[0, 0] + pns[i, 1]*J[1, 0])
                           + rns[i, 1]*(pns[i, 0]*J[0, 1] + pns[i, 1]*J[1, 1]))

        for i, c in enumerate([(1, 2), (0, 2), (0, 1)]):
            B12 = (rns[i, 0]*(pts[i, 0]*J[0, 0] + pts[i, 1]*J[1, 0])
                   + rns[i, 1]*(pts[i, 0]*J[0, 1] + pts[i, 1]*J[1, 1]))

            V[3+i, c[0]] = -1*B12 / pel[i]