# Path: gusto/kernels.py
# class GaussianElimination(object):
# class Average(object):
# class AverageWeightings(object):
# class PhysicsRecoveryTop():
# class PhysicsRecoveryBottom():
#     def __init__(self, DG1):
#     def apply(self, v_DG1_old, v_DG1, act_coords, eff_coords, num_ext):
#     def __init__(self, V):
#     def apply(self, v_out, weighting, v_in):
#     def __init__(self, V):
#     def apply(self, w):
#     def __init__(self):
#     def apply(self, v_DG1, v_CG1):
#     def __init__(self):
#     def apply(self, v_DG1, v_CG1):


from enum import Enum
from firedrake import (BrokenElement, Constant, DirichletBC, FiniteElement,
                       Function, FunctionSpace, Interpolator, Projector,
                       SpatialCoordinate, TensorProductElement,
                       VectorFunctionSpace, as_vector, function, interval)
from firedrake.utils import cached_property
from gusto import kernels
import ufl
                raise NotImplementedError('The physics boundary method only works on extruded meshes')
            # base spaces
            cell = mesh._base_mesh.ufl_cell().cellname()
            w_hori = FiniteElement("DG", cell, 0, variant="equispaced")
            w_vert = FiniteElement("CG", interval, 1, variant="equispaced")
            # build element
            theta_element = TensorProductElement(w_hori, w_vert)
            # spaces
            Vtheta = FunctionSpace(mesh, theta_element)
            Vtheta_broken = FunctionSpace(mesh, BrokenElement(theta_element))
            if v_CG1.function_space() != Vtheta:
                raise ValueError("This boundary recovery method requires v_CG1 to be in DG0xCG1 TensorProductSpace.")
            if v_DG1.function_space() != Vtheta_broken:
                raise ValueError("This boundary recovery method requires v_DG1 to be in the broken DG0xCG1 TensorProductSpace.")
        else:
            raise ValueError("Boundary method should be a Boundary Method Enum object.")

        vec_DG1 = VectorFunctionSpace(DG0.mesh(), DG1_element)
        x = SpatialCoordinate(DG0.mesh())
        self.interpolator = Interpolator(self.v_CG1, self.v_DG1)

        if self.method == Boundary_Method.dynamics:

            # STRATEGY
            # obtain a coordinate field for all the nodes
            self.act_coords = Function(vec_DG1).project(x)  # actual coordinates
            self.eff_coords = eff_coords  # effective coordinates
            self.output = Function(DG1)
            self.on_exterior = find_domain_boundaries(mesh)


# Path: gusto/kernels.py
# class GaussianElimination(object):
# class Average(object):
# class AverageWeightings(object):
# class PhysicsRecoveryTop():
# class PhysicsRecoveryBottom():
#     def __init__(self, DG1):
#     def apply(self, v_DG1_old, v_DG1, act_coords, eff_coords, num_ext):
#     def __init__(self, V):
#     def apply(self, v_out, weighting, v_in):
#     def __init__(self, V):
#     def apply(self, w):
#     def __init__(self):
#     def apply(self, v_DG1, v_CG1):
#     def __init__(self):
#     def apply(self, v_DG1, v_CG1):


# Path: gusto/recovery.py
from enum import Enum
from firedrake import (BrokenElement, Constant, DirichletBC, FiniteElement,
                       Function, FunctionSpace, Interpolator, Projector,
                       SpatialCoordinate, TensorProductElement,
                       VectorFunctionSpace, as_vector, function, interval)
from firedrake.utils import cached_property
from gusto import kernels
import ufl

                raise NotImplementedError('The physics boundary method only works on extruded meshes')
            # base spaces
            cell = mesh._base_mesh.ufl_cell().cellname()
            w_hori = FiniteElement("DG", cell, 0, variant="equispaced")
            w_vert = FiniteElement("CG", interval, 1, variant="equispaced")
            # build element
            theta_element = TensorProductElement(w_hori, w_vert)
            # spaces
            Vtheta = FunctionSpace(mesh, theta_element)
            Vtheta_broken = FunctionSpace(mesh, BrokenElement(theta_element))
            if v_CG1.function_space() != Vtheta:
                raise ValueError("This boundary recovery method requires v_CG1 to be in DG0xCG1 TensorProductSpace.")
            if v_DG1.function_space() != Vtheta_broken:
                raise ValueError("This boundary recovery method requires v_DG1 to be in the broken DG0xCG1 TensorProductSpace.")
        else:
            raise ValueError("Boundary method should be a Boundary Method Enum object.")

        vec_DG1 = VectorFunctionSpace(DG0.mesh(), DG1_element)
        x = SpatialCoordinate(DG0.mesh())
        self.interpolator = Interpolator(self.v_CG1, self.v_DG1)

        if self.method == Boundary_Method.dynamics:

            # STRATEGY
            # obtain a coordinate field for all the nodes
            self.act_coords = Function(vec_DG1).project(x)  # actual coordinates
            self.eff_coords = eff_coords  # effective coordinates
            self.output = Function(DG1)
            self.on_exterior = find_domain_boundaries(mesh)


            self.gaussian_elimination_kernel = kernels.GaussianElimination(DG1)