<commit_before>"""
The recovery operators used for lowest-order advection schemes.
"""
from firedrake import (expression, function, Function, FunctionSpace, Projector,
                       VectorFunctionSpace, SpatialCoordinate, as_vector,
                       dx, Interpolator, BrokenElement, interval, Constant,
                       TensorProductElement, FiniteElement, DirichletBC,
                       VectorElement)
from firedrake.utils import cached_property
from firedrake.parloops import par_loop, READ, INC, WRITE
from gusto.transport_equation import is_dg
from pyop2 import ON_TOP, ON_BOTTOM
import ufl
import numpy as np
from enum import Enum

__all__ = ["Averager", "Boundary_Method", "Boundary_Recoverer", "Recoverer"]


class Averager(object):
    """
    An object that 'recovers' a low order field (e.g. in DG0)
    into a higher order field (e.g. in CG1).
    The code is essentially that of the Firedrake Projector
    object, using the "average" method, and could possibly
    be replaced by it if it comes into the master branch.

    :arg v: the :class:`ufl.Expr` or
         :class:`.Function` to project.
    :arg v_out: :class:`.Function` to put the result in.
    """

    def __init__(self, v, v_out):

        if isinstance(v, expression.Expression) or not isinstance(v, (ufl.core.expr.Expr, function.Function)):
            raise ValueError("Can only recover UFL expression or Functions not '%s'" % type(v))

        # Check shape values
        if v.ufl_shape != v_out.ufl_shape:
            raise RuntimeError('Shape mismatch between source %s and target function spaces %s in project' % (v.ufl_shape, v_out.ufl_shape))

        self._same_fspace = (isinstance(v, function.Function) and v.function_space() == v_out.function_space())
        self.v = v
        self.v_out = v_out
        self.V = v_out.function_space()

        # Check the number of local dofs
        if self.v_out.function_space().finat_element.space_dimension() != self.v.function_space().finat_element.space_dimension():
            raise RuntimeError("Number of local dofs for each field must be equal.")

        # NOTE: Any bcs on the function self.v should just work.
        # Loop over node extent and dof extent
        self.shapes = {"nDOFs": self.V.finat_element.space_dimension(),
                       "dim": np.prod(self.V.shape, dtype=int)}
        # Averaging kernel
        average_domain = "{{[i, j]: 0 <= i < {nDOFs} and 0 <= j < {dim}}}".format(**self.shapes)
        average_instructions = ("""
                                for i
                                    for j
                                        vo[i,j] = vo[i,j] + v[i,j] / w[i,j]
                                    end
                                end
                                """)
        self._average_kernel = (average_domain, average_instructions)

    @cached_property
    def _weighting(self):
        """
        Generates a weight function for computing a projection via averaging.
        """
        w = Function(self.V)
        weight_domain = "{{[i, j]: 0 <= i < {nDOFs} and 0 <= j < {dim}}}".format(**self.shapes)
        weight_instructions = ("""
                               for i
                                   for j
                                      w[i,j] = w[i,j] + 1.0
                                   end
                               end
                               """)
        _weight_kernel = (weight_domain, weight_instructions)

        par_loop(_weight_kernel, dx, {"w": (w, INC)}, is_loopy_kernel=True)
        return w

    def project(self):
        """
        Apply the recovery.
        """

        # Ensure that the function being populated is zeroed out
        self.v_out.dat.zero()
        par_loop(self._average_kernel, dx, {"vo": (self.v_out, INC),
                                            "w": (self._weighting, READ),
                                            "v": (self.v, READ)},
                 is_loopy_kernel=True)
        return self.v_out


class Boundary_Method(Enum):
    """
    An Enum object storing the two types of boundary method:
    dynamics -- which corrects a field recovered into CG1.
    physics -- which corrects a field recovered into the temperature space.
    """

    dynamics = 0
    physics = 1


class Boundary_Recoverer(object):
    """
    An object that performs a `recovery` process at the domain
    boundaries that has second order accuracy. This is necessary
    because the :class:`Averager` object does not recover a field
    with sufficient accuracy at the boundaries.

    The strategy is to minimise the curvature of the function in
    the boundary cells, subject to the constraints of conserved
    mass and continuity on the interior facets. The quickest way
    to perform this is by using the analytic solution and a parloop.

    Currently this is only implemented for the (DG0, DG1, CG1)
    set of spaces, and only on a `PeriodicIntervalMesh` or
    'PeriodicUnitIntervalMesh` that has been extruded.

    :arg v_CG1: the continuous function after the first recovery
             is performed. Should be in CG1. This is correct
             on the interior of the domain.
    :arg v_DG1: the function to be output. Should be in DG1.
    :arg method: a Boundary_Method Enum object.
    :arg coords_to_adjust: a DG1 field containing 1 at locations of
                           coords that must be adjusted to give they
                           effective coords.
    """

    def __init__(self, v_CG1, v_DG1, method=Boundary_Method.physics, coords_to_adjust=None):

        self.v_DG1 = v_DG1
        self.v_CG1 = v_CG1
        self.v_DG1_old = Function(v_DG1.function_space())
        self.coords_to_adjust = coords_to_adjust

        self.method = method
        mesh = v_CG1.function_space().mesh()
        VDG0 = FunctionSpace(mesh, "DG", 0)
        VCG1 = FunctionSpace(mesh, "CG", 1)
        VDG1 = FunctionSpace(mesh, "DG", 1)

        self.num_ext = Function(VDG0)

        # check function spaces of functions
        if self.method == Boundary_Method.dynamics:
            if v_CG1.function_space() != VCG1:
                raise NotImplementedError("This boundary recovery method requires v1 to be in CG1.")
            if v_DG1.function_space() != VDG1:
                raise NotImplementedError("This boundary recovery method requires v_out to be in DG1.")
            # check whether mesh is valid
            if mesh.topological_dimension() == 2:
                # if mesh is extruded then we're fine, but if not needs to be quads
                if not VDG0.extruded and mesh.ufl_cell().cell_name() != 'quadrilateral':
                    raise NotImplementedError('For 2D meshes this recovery method requires that elements are quadrilaterals')
            elif mesh.topological_dimension() == 3:
                # assume that 3D mesh is extruded
                if mesh._base_mesh.ufl_cell().cellname() != 'quadrilateral':
                    raise NotImplementedError('For 3D extruded meshes this recovery method requires a base mesh with quadrilateral elements')
            elif mesh.topological_dimension() != 1:
                raise NotImplementedError('This boundary recovery is implemented only on certain classes of mesh.')
            if coords_to_adjust is None:
                raise ValueError('Need coords_to_adjust field for dynamics boundary methods')

        elif self.method == Boundary_Method.physics:
            # check that mesh is valid -- must be an extruded mesh
            if not VDG0.extruded:
                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)
            w_vert = FiniteElement("CG", interval, 1)
            # 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.")

        VuDG1 = VectorFunctionSpace(VDG0.mesh(), "DG", 1)
        x = SpatialCoordinate(VDG0.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
            VuDG1 = VectorFunctionSpace(mesh, "DG", 1)
            self.act_coords = Function(VuDG1).project(x)  # actual coordinates
            self.eff_coords = Function(VuDG1).project(x)  # effective coordinates

            shapes = {"nDOFs": self.v_DG1.function_space().finat_element.space_dimension(),
                      "dim": np.prod(VuDG1.shape, dtype=int)}

            num_ext_domain = ("{{[i]: 0 <= i < {nDOFs}}}").format(**shapes)
            num_ext_instructions = ("""
            <float64> SUM_EXT = 0
            for i
                SUM_EXT = SUM_EXT + EXT_V1[i]
            end

            NUM_EXT[0] = SUM_EXT
            """)

            coords_domain = ("{{[i, j, k, ii, jj, kk, ll, mm, iii, kkk]: "
                             "0 <= i < {nDOFs} and "
                             "0 <= j < {nDOFs} and 0 <= k < {dim} and "
                             "0 <= ii < {nDOFs} and 0 <= jj < {nDOFs} and "
                             "0 <= kk < {dim} and 0 <= ll < {dim} and "
                             "0 <= mm < {dim} and 0 <= iii < {nDOFs} and "
                             "0 <= kkk < {dim}}}").format(**shapes)
            coords_insts = ("""
                            <float64> sum_V1_ext = 0
                            <int> index = 100
                            <float64> dist = 0.0
                            <float64> max_dist = 0.0
                            <float64> min_dist = 0.0
                            """
                            # only do adjustment in cells with at least one DOF to adjust
                            """
                            if NUM_EXT[0] > 0
                            """
                            # find the maximum distance between DOFs in this cell, to serve as starting point for finding min distances
                            """
                                for i
                                    for j
                                        dist = 0.0
                                        for k
                                            dist = dist + pow(ACT_COORDS[i,k] - ACT_COORDS[j,k], 2.0)
                                        end
                                        dist = pow(dist, 0.5) {{id=sqrt_max_dist, dep=*}}
                                        max_dist = fmax(dist, max_dist) {{id=max_dist, dep=sqrt_max_dist}}
                                    end
                                end
                            """
                            # loop through cells and find which ones to adjust
                            """
                                for ii
                                    if EXT_V1[ii] > 0.5
                            """
                            # find closest interior node
                            """
                                        min_dist = max_dist
                                        index = 100
                                        for jj
                                            if EXT_V1[jj] < 0.5
                                                dist = 0.0
                                                for kk
                                                    dist = dist + pow(ACT_COORDS[ii,kk] - ACT_COORDS[jj,kk], 2)
                                                end
                                                dist = pow(dist, 0.5)
                                                if dist <= min_dist
                                                    index = jj
                                                end
                                                min_dist = fmin(min_dist, dist)
                                                for ll
                                                    EFF_COORDS[ii,ll] = 0.5 * (ACT_COORDS[ii,ll] + ACT_COORDS[index,ll])
                                                end
                                            end
                                        end
                                    else
                            """
                            # for DOFs that aren't exterior, use the original coordinates
                            """
                                        for mm
                                            EFF_COORDS[ii, mm] = ACT_COORDS[ii, mm]
                                        end
                                    end
                                end
                            else
                            """
                            # for interior elements, just use the original coordinates
                            """
                                for iii
                                    for kkk
                                        EFF_COORDS[iii, kkk] = ACT_COORDS[iii, kkk]
                                    end
                                end
                            end
                            """).format(**shapes)

            elimin_domain = ("{{[i, ii_loop, jj_loop, kk, ll_loop, mm, iii_loop, kkk_loop, iiii, iiiii]: "
                             "0 <= i < {nDOFs} and 0 <= ii_loop < {nDOFs} and "
                             "ii_loop + 1 <= jj_loop < {nDOFs} and ii_loop <= kk < {nDOFs} and "
                             "ii_loop + 1 <= ll_loop < {nDOFs} and ii_loop <= mm < {nDOFs} + 1 and "
                             "0 <= iii_loop < {nDOFs} and {nDOFs} - iii_loop <= kkk_loop < {nDOFs} + 1 and "
                             "0 <= iiii < {nDOFs} and 0 <= iiiii < {nDOFs}}}").format(**shapes)
            elimin_insts = ("""
                            <int> ii = 0
                            <int> jj = 0
                            <int> ll = 0
                            <int> iii = 0
                            <int> jjj = 0
                            <int> i_max = 0
                            <float64> A_max = 0.0
                            <float64> temp_f = 0.0
                            <float64> temp_A = 0.0
                            <float64> c = 0.0
                            <float64> f[{nDOFs}] = 0.0
                            <float64> a[{nDOFs}] = 0.0
                            <float64> A[{nDOFs},{nDOFs}] = 0.0
                            """
                            # We are aiming to find the vector a that solves A*a = f, for matrix A and vector f.
                            # This is done by performing row operations (swapping and scaling) to obtain A in upper diagonal form.
                            # N.B. several for loops must be executed in numerical order (loopy does not necessarily do this).
                            # For these loops we must manually iterate the index.
                            """
                            if NUM_EXT[0] > 0.0
                            """
                            # only do Gaussian elimination for elements with effective coordinates
                            """
                                for i
                            """
                            # fill f with the original field values and A with the effective coordinate values
                            """
                                    f[i] = DG1_OLD[i]
                                    A[i,0] = 1.0
                                    A[i,1] = EFF_COORDS[i,0]
                                    if {nDOFs} > 3
                                        A[i,2] = EFF_COORDS[i,1]
                                        A[i,3] = EFF_COORDS[i,0]*EFF_COORDS[i,1]
                                        if {nDOFs} > 7
                                            A[i,4] = EFF_COORDS[i,2]
                                            A[i,5] = EFF_COORDS[i,0]*EFF_COORDS[i,2]
                                            A[i,6] = EFF_COORDS[i,1]*EFF_COORDS[i,2]
                                            A[i,7] = EFF_COORDS[i,0]*EFF_COORDS[i,1]*EFF_COORDS[i,2]
                                        end
                                    end
                                end
                            """
                            # now loop through rows/columns of A
                            """
                                for ii_loop
                                    A_max = fabs(A[ii,ii])
                                    i_max = ii
                            """
                            # loop to find the largest value in the ith column
                            # set i_max as the index of the row with this largest value.
                            """
                                    jj = ii + 1
                                    for jj_loop
                                        if fabs(A[jj,ii]) > A_max
                                            i_max = jj
                                        end
                                        A_max = fmax(A_max, fabs(A[jj,ii]))
                                        jj = jj + 1
                                    end
                            """
                            # if the max value in the ith column isn't in the ith row, we must swap the rows
                            """
                                    if i_max != ii
                            """
                            # swap the elements of f
                            """
                                        temp_f = f[ii]  {{id=set_temp_f, dep=*}}
                                        f[ii] = f[i_max]  {{id=set_f_imax, dep=set_temp_f}}
                                        f[i_max] = temp_f  {{id=set_f_ii, dep=set_f_imax}}
                            """
                            # swap the elements of A
                            # N.B. kk runs from ii to (nDOFs-1) as elements below diagonal should be 0
                            """
                                        for kk
                                            temp_A = A[ii,kk]  {{id=set_temp_A, dep=*}}
                                            A[ii, kk] = A[i_max, kk]  {{id=set_A_ii, dep=set_temp_A}}
                                            A[i_max, kk] = temp_A  {{id=set_A_imax, dep=set_A_ii}}
                                        end
                                    end
                            """
                            # scale the rows below the ith row
                            """
                                    ll = ii + 1
                                    for ll_loop
                                        if ll > ii
                            """
                            # find scaling factor
                            """
                                            c = - A[ll,ii] / A[ii,ii]
                            """
                            # N.B. mm runs from ii to (nDOFs-1) as elements below diagonal should be 0
                            """
                                            for mm
                                                A[ll, mm] = A[ll, mm] + c * A[ii,mm]
                                            end
                                            f[ll] = f[ll] + c * f[ii]
                                        end
                                        ll = ll + 1
                                    end
                                    ii = ii + 1
                                end
                            """
                            # do back substitution of upper diagonal A to obtain a
                            """
                                iii = 0
                                for iii_loop
                            """
                            # jjj starts at the bottom row and works upwards
                            """
                                    jjj = {nDOFs} - iii - 1  {{id=assign_jjj, dep=*}}
                                    a[jjj] = f[jjj]   {{id=set_a, dep=assign_jjj}}
                                    for kkk_loop
                                        a[jjj] = a[jjj] - A[jjj,kkk_loop] * a[kkk_loop]
                                    end
                                    a[jjj] = a[jjj] / A[jjj,jjj]
                                    iii = iii + 1
                                end
                            """
                            # Having found a, this gives us the coefficients for the Taylor expansion with the actual coordinates.
                            """
                                for iiii
                                    if {nDOFs} == 2
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0]
                                    elif {nDOFs} == 4
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0] + a[2]*ACT_COORDS[iiii,1] + a[3]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1]
                                    elif {nDOFs} == 8
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0] + a[2]*ACT_COORDS[iiii,1] + a[3]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1] + a[4]*ACT_COORDS[iiii,2] + a[5]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,2] + a[6]*ACT_COORDS[iiii,1]*ACT_COORDS[iiii,2] + a[7]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1]*ACT_COORDS[iiii,2]
                                    end
                                end
                            """
                            # if element is not external, just use old field values.
                            """
                            else
                                for iiiii<commit_msg>correct cellname and update gauss elimination code to avoid errors with sw setups<commit_after>"""
The recovery operators used for lowest-order advection schemes.
"""
from firedrake import (expression, function, Function, FunctionSpace, Projector,
                       VectorFunctionSpace, SpatialCoordinate, as_vector,
                       dx, Interpolator, BrokenElement, interval, Constant,
                       TensorProductElement, FiniteElement, DirichletBC,
                       VectorElement)
from firedrake.utils import cached_property
from firedrake.parloops import par_loop, READ, INC, WRITE
from gusto.transport_equation import is_dg
from pyop2 import ON_TOP, ON_BOTTOM
import ufl
import numpy as np
from enum import Enum

__all__ = ["Averager", "Boundary_Method", "Boundary_Recoverer", "Recoverer"]


class Averager(object):
    """
    An object that 'recovers' a low order field (e.g. in DG0)
    into a higher order field (e.g. in CG1).
    The code is essentially that of the Firedrake Projector
    object, using the "average" method, and could possibly
    be replaced by it if it comes into the master branch.

    :arg v: the :class:`ufl.Expr` or
         :class:`.Function` to project.
    :arg v_out: :class:`.Function` to put the result in.
    """

    def __init__(self, v, v_out):

        if isinstance(v, expression.Expression) or not isinstance(v, (ufl.core.expr.Expr, function.Function)):
            raise ValueError("Can only recover UFL expression or Functions not '%s'" % type(v))

        # Check shape values
        if v.ufl_shape != v_out.ufl_shape:
            raise RuntimeError('Shape mismatch between source %s and target function spaces %s in project' % (v.ufl_shape, v_out.ufl_shape))

        self._same_fspace = (isinstance(v, function.Function) and v.function_space() == v_out.function_space())
        self.v = v
        self.v_out = v_out
        self.V = v_out.function_space()

        # Check the number of local dofs
        if self.v_out.function_space().finat_element.space_dimension() != self.v.function_space().finat_element.space_dimension():
            raise RuntimeError("Number of local dofs for each field must be equal.")

        # NOTE: Any bcs on the function self.v should just work.
        # Loop over node extent and dof extent
        self.shapes = {"nDOFs": self.V.finat_element.space_dimension(),
                       "dim": np.prod(self.V.shape, dtype=int)}
        # Averaging kernel
        average_domain = "{{[i, j]: 0 <= i < {nDOFs} and 0 <= j < {dim}}}".format(**self.shapes)
        average_instructions = ("""
                                for i
                                    for j
                                        vo[i,j] = vo[i,j] + v[i,j] / w[i,j]
                                    end
                                end
                                """)
        self._average_kernel = (average_domain, average_instructions)

    @cached_property
    def _weighting(self):
        """
        Generates a weight function for computing a projection via averaging.
        """
        w = Function(self.V)
        weight_domain = "{{[i, j]: 0 <= i < {nDOFs} and 0 <= j < {dim}}}".format(**self.shapes)
        weight_instructions = ("""
                               for i
                                   for j
                                      w[i,j] = w[i,j] + 1.0
                                   end
                               end
                               """)
        _weight_kernel = (weight_domain, weight_instructions)

        par_loop(_weight_kernel, dx, {"w": (w, INC)}, is_loopy_kernel=True)
        return w

    def project(self):
        """
        Apply the recovery.
        """

        # Ensure that the function being populated is zeroed out
        self.v_out.dat.zero()
        par_loop(self._average_kernel, dx, {"vo": (self.v_out, INC),
                                            "w": (self._weighting, READ),
                                            "v": (self.v, READ)},
                 is_loopy_kernel=True)
        return self.v_out


class Boundary_Method(Enum):
    """
    An Enum object storing the two types of boundary method:
    dynamics -- which corrects a field recovered into CG1.
    physics -- which corrects a field recovered into the temperature space.
    """

    dynamics = 0
    physics = 1


class Boundary_Recoverer(object):
    """
    An object that performs a `recovery` process at the domain
    boundaries that has second order accuracy. This is necessary
    because the :class:`Averager` object does not recover a field
    with sufficient accuracy at the boundaries.

    The strategy is to minimise the curvature of the function in
    the boundary cells, subject to the constraints of conserved
    mass and continuity on the interior facets. The quickest way
    to perform this is by using the analytic solution and a parloop.

    Currently this is only implemented for the (DG0, DG1, CG1)
    set of spaces, and only on a `PeriodicIntervalMesh` or
    'PeriodicUnitIntervalMesh` that has been extruded.

    :arg v_CG1: the continuous function after the first recovery
             is performed. Should be in CG1. This is correct
             on the interior of the domain.
    :arg v_DG1: the function to be output. Should be in DG1.
    :arg method: a Boundary_Method Enum object.
    :arg coords_to_adjust: a DG1 field containing 1 at locations of
                           coords that must be adjusted to give they
                           effective coords.
    """

    def __init__(self, v_CG1, v_DG1, method=Boundary_Method.physics, coords_to_adjust=None):

        self.v_DG1 = v_DG1
        self.v_CG1 = v_CG1
        self.v_DG1_old = Function(v_DG1.function_space())
        self.coords_to_adjust = coords_to_adjust

        self.method = method
        mesh = v_CG1.function_space().mesh()
        VDG0 = FunctionSpace(mesh, "DG", 0)
        VCG1 = FunctionSpace(mesh, "CG", 1)
        VDG1 = FunctionSpace(mesh, "DG", 1)

        self.num_ext = Function(VDG0)

        # check function spaces of functions
        if self.method == Boundary_Method.dynamics:
            if v_CG1.function_space() != VCG1:
                raise NotImplementedError("This boundary recovery method requires v1 to be in CG1.")
            if v_DG1.function_space() != VDG1:
                raise NotImplementedError("This boundary recovery method requires v_out to be in DG1.")
            # check whether mesh is valid
            if mesh.topological_dimension() == 2:
                # if mesh is extruded then we're fine, but if not needs to be quads
                if not VDG0.extruded and mesh.ufl_cell().cellname() != 'quadrilateral':
                    raise NotImplementedError('For 2D meshes this recovery method requires that elements are quadrilaterals')
            elif mesh.topological_dimension() == 3:
                # assume that 3D mesh is extruded
                if mesh._base_mesh.ufl_cell().cellname() != 'quadrilateral':
                    raise NotImplementedError('For 3D extruded meshes this recovery method requires a base mesh with quadrilateral elements')
            elif mesh.topological_dimension() != 1:
                raise NotImplementedError('This boundary recovery is implemented only on certain classes of mesh.')
            if coords_to_adjust is None:
                raise ValueError('Need coords_to_adjust field for dynamics boundary methods')

        elif self.method == Boundary_Method.physics:
            # check that mesh is valid -- must be an extruded mesh
            if not VDG0.extruded:
                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)
            w_vert = FiniteElement("CG", interval, 1)
            # 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.")

        VuDG1 = VectorFunctionSpace(VDG0.mesh(), "DG", 1)
        x = SpatialCoordinate(VDG0.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
            VuDG1 = VectorFunctionSpace(mesh, "DG", 1)
            self.act_coords = Function(VuDG1).project(x)  # actual coordinates
            self.eff_coords = Function(VuDG1).project(x)  # effective coordinates

            shapes = {"nDOFs": self.v_DG1.function_space().finat_element.space_dimension(),
                      "dim": np.prod(VuDG1.shape, dtype=int)}

            num_ext_domain = ("{{[i]: 0 <= i < {nDOFs}}}").format(**shapes)
            num_ext_instructions = ("""
            <float64> SUM_EXT = 0
            for i
                SUM_EXT = SUM_EXT + EXT_V1[i]
            end

            NUM_EXT[0] = SUM_EXT
            """)

            coords_domain = ("{{[i, j, k, ii, jj, kk, ll, mm, iii, kkk]: "
                             "0 <= i < {nDOFs} and "
                             "0 <= j < {nDOFs} and 0 <= k < {dim} and "
                             "0 <= ii < {nDOFs} and 0 <= jj < {nDOFs} and "
                             "0 <= kk < {dim} and 0 <= ll < {dim} and "
                             "0 <= mm < {dim} and 0 <= iii < {nDOFs} and "
                             "0 <= kkk < {dim}}}").format(**shapes)
            coords_insts = ("""
                            <float64> sum_V1_ext = 0
                            <int> index = 100
                            <float64> dist = 0.0
                            <float64> max_dist = 0.0
                            <float64> min_dist = 0.0
                            """
                            # only do adjustment in cells with at least one DOF to adjust
                            """
                            if NUM_EXT[0] > 0
                            """
                            # find the maximum distance between DOFs in this cell, to serve as starting point for finding min distances
                            """
                                for i
                                    for j
                                        dist = 0.0
                                        for k
                                            dist = dist + pow(ACT_COORDS[i,k] - ACT_COORDS[j,k], 2.0)
                                        end
                                        dist = pow(dist, 0.5) {{id=sqrt_max_dist, dep=*}}
                                        max_dist = fmax(dist, max_dist) {{id=max_dist, dep=sqrt_max_dist}}
                                    end
                                end
                            """
                            # loop through cells and find which ones to adjust
                            """
                                for ii
                                    if EXT_V1[ii] > 0.5
                            """
                            # find closest interior node
                            """
                                        min_dist = max_dist
                                        index = 100
                                        for jj
                                            if EXT_V1[jj] < 0.5
                                                dist = 0.0
                                                for kk
                                                    dist = dist + pow(ACT_COORDS[ii,kk] - ACT_COORDS[jj,kk], 2)
                                                end
                                                dist = pow(dist, 0.5)
                                                if dist <= min_dist
                                                    index = jj
                                                end
                                                min_dist = fmin(min_dist, dist)
                                                for ll
                                                    EFF_COORDS[ii,ll] = 0.5 * (ACT_COORDS[ii,ll] + ACT_COORDS[index,ll])
                                                end
                                            end
                                        end
                                    else
                            """
                            # for DOFs that aren't exterior, use the original coordinates
                            """
                                        for mm
                                            EFF_COORDS[ii, mm] = ACT_COORDS[ii, mm]
                                        end
                                    end
                                end
                            else
                            """
                            # for interior elements, just use the original coordinates
                            """
                                for iii
                                    for kkk
                                        EFF_COORDS[iii, kkk] = ACT_COORDS[iii, kkk]
                                    end
                                end
                            end
                            """).format(**shapes)

            elimin_domain = ("{{[i, ii_loop, jj_loop, kk, ll_loop, mm, iii_loop, kkk_loop, iiii, iiiii]: "
                             "0 <= i < {nDOFs} and 0 <= ii_loop < {nDOFs} and "
                             "ii_loop + 1 <= jj_loop < {nDOFs} and ii_loop <= kk < {nDOFs} and "
                             "ii_loop + 1 <= ll_loop < {nDOFs} and ii_loop <= mm < {nDOFs} + 1 and "
                             "0 <= iii_loop < {nDOFs} and {nDOFs} - iii_loop <= kkk_loop < {nDOFs} + 1 and "
                             "0 <= iiii < {nDOFs} and 0 <= iiiii < {nDOFs}}}").format(**shapes)
            elimin_insts = ("""
                            <int> ii = 0
                            <int> jj = 0
                            <int> ll = 0
                            <int> iii = 0
                            <int> jjj = 0
                            <int> i_max = 0
                            <float64> A_max = 0.0
                            <float64> temp_f = 0.0
                            <float64> temp_A = 0.0
                            <float64> c = 0.0
                            <float64> f[{nDOFs}] = 0.0
                            <float64> a[{nDOFs}] = 0.0
                            <float64> A[{nDOFs},{nDOFs}] = 0.0
                            """
                            # We are aiming to find the vector a that solves A*a = f, for matrix A and vector f.
                            # This is done by performing row operations (swapping and scaling) to obtain A in upper diagonal form.
                            # N.B. several for loops must be executed in numerical order (loopy does not necessarily do this).
                            # For these loops we must manually iterate the index.
                            """
                            if NUM_EXT[0] > 0.0
                            """
                            # only do Gaussian elimination for elements with effective coordinates
                            """
                                for i
                            """
                            # fill f with the original field values and A with the effective coordinate values
                            """
                                    f[i] = DG1_OLD[i]
                                    A[i,0] = 1.0
                                    A[i,1] = EFF_COORDS[i,0]
                                    if {nDOFs} > 3
                                        A[i,2] = EFF_COORDS[i,1]
                                        A[i,3] = EFF_COORDS[i,0]*EFF_COORDS[i,1]
                                        if {nDOFs} > 7
                                            A[i,4] = EFF_COORDS[i,{dim}-1]
                                            A[i,5] = EFF_COORDS[i,0]*EFF_COORDS[i,{dim}-1]
                                            A[i,6] = EFF_COORDS[i,1]*EFF_COORDS[i,{dim}-1]
                                            A[i,7] = EFF_COORDS[i,0]*EFF_COORDS[i,1]*EFF_COORDS[i,{dim}-1]
                                        end
                                    end
                                end
                            """
                            # now loop through rows/columns of A
                            """
                                for ii_loop
                                    A_max = fabs(A[ii,ii])
                                    i_max = ii
                            """
                            # loop to find the largest value in the ith column
                            # set i_max as the index of the row with this largest value.
                            """
                                    jj = ii + 1
                                    for jj_loop
                                        if fabs(A[jj,ii]) > A_max
                                            i_max = jj
                                        end
                                        A_max = fmax(A_max, fabs(A[jj,ii]))
                                        jj = jj + 1
                                    end
                            """
                            # if the max value in the ith column isn't in the ith row, we must swap the rows
                            """
                                    if i_max != ii
                            """
                            # swap the elements of f
                            """
                                        temp_f = f[ii]  {{id=set_temp_f, dep=*}}
                                        f[ii] = f[i_max]  {{id=set_f_imax, dep=set_temp_f}}
                                        f[i_max] = temp_f  {{id=set_f_ii, dep=set_f_imax}}
                            """
                            # swap the elements of A
                            # N.B. kk runs from ii to (nDOFs-1) as elements below diagonal should be 0
                            """
                                        for kk
                                            temp_A = A[ii,kk]  {{id=set_temp_A, dep=*}}
                                            A[ii, kk] = A[i_max, kk]  {{id=set_A_ii, dep=set_temp_A}}
                                            A[i_max, kk] = temp_A  {{id=set_A_imax, dep=set_A_ii}}
                                        end
                                    end
                            """
                            # scale the rows below the ith row
                            """
                                    ll = ii + 1
                                    for ll_loop
                                        if ll > ii
                            """
                            # find scaling factor
                            """
                                            c = - A[ll,ii] / A[ii,ii]
                            """
                            # N.B. mm runs from ii to (nDOFs-1) as elements below diagonal should be 0
                            """
                                            for mm
                                                A[ll, mm] = A[ll, mm] + c * A[ii,mm]
                                            end
                                            f[ll] = f[ll] + c * f[ii]
                                        end
                                        ll = ll + 1
                                    end
                                    ii = ii + 1
                                end
                            """
                            # do back substitution of upper diagonal A to obtain a
                            """
                                iii = 0
                                for iii_loop
                            """
                            # jjj starts at the bottom row and works upwards
                            """
                                    jjj = {nDOFs} - iii - 1  {{id=assign_jjj, dep=*}}
                                    a[jjj] = f[jjj]   {{id=set_a, dep=assign_jjj}}
                                    for kkk_loop
                                        a[jjj] = a[jjj] - A[jjj,kkk_loop] * a[kkk_loop]
                                    end
                                    a[jjj] = a[jjj] / A[jjj,jjj]
                                    iii = iii + 1
                                end
                            """
                            # Having found a, this gives us the coefficients for the Taylor expansion with the actual coordinates.
                            """
                                for iiii
                                    if {nDOFs} == 2
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0]
                                    elif {nDOFs} == 4
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0] + a[2]*ACT_COORDS[iiii,1] + a[3]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1]
                                    elif {nDOFs} == 8
                                        DG1[iiii] = a[0] + a[1]*ACT_COORDS[iiii,0] + a[2]*ACT_COORDS[iiii,1] + a[3]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1] + a[4]*ACT_COORDS[iiii,{dim}-1] + a[5]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,{dim}-1] + a[6]*ACT_COORDS[iiii,1]*ACT_COORDS[iiii,{dim}-1] + a[7]*ACT_COORDS[iiii,0]*ACT_COORDS[iiii,1]*ACT_COORDS[iiii,{dim}-1]
                                    end
                                end
                            """
                            # if element is not external, just use old field values.
                            """
                            else
                                for iiiii
{'fraction_non_alphanumeric': 0.08090723682736053, 'fraction_numerical': 0.011340460342006769, 'mean_word_length': 0.9569317135213368, 'pattern_counts': {'":': 16, '<': 145, '<?xml version=': 0, '>': 55, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 318, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3552738', 'n_tokens_mistral': 11488, 'n_tokens_neox': 10848, 'n_words': 3133}