import numpy as np
from scipy import sparse as sp
from SimPEG.Utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap


def checkBC(bc):
    """

        Checks if boundary condition 'bc' is valid.

        Each bc must be either 'dirichlet' or 'neumann'

    """
    if(type(bc) is str):
        bc = [bc, bc]
    assert type(bc) is list, 'bc must be a list'
    assert len(bc) == 2, 'bc must have two elements'

    for bc_i in bc:
        assert type(bc_i) is str, "each bc must be a string"
        assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'"
    return bc


def ddxCellGrad(n, bc):
    """
        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1

        For Cell-Centered **Dirichlet**, use a ghost point::

            (u_1 - u_g)/hf = grad

                u_g       u_1      u_2
                 *    |    *   |    *     ...
                      ^
                      0

            u_g = - u_1
            grad = 2*u1/dx
            negitive on the other side.

        For Cell-Centered **Neumann**, use a ghost point::

            (u_1 - u_g)/hf = 0

                u_g       u_1      u_2
                 *    |    *   |    *     ...

            u_g = u_1
            grad = 0;  put a zero in.

    """
    bc = checkBC(bc)

    D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
    # Set the first side
    if(bc[0] == 'dirichlet'):
        D[0, 0] = 2
    elif(bc[0] == 'neumann'):
        D[0, 0] = 0
    # Set the second side
    if(bc[1] == 'dirichlet'):
        D[-1, -1] = -2
    elif(bc[1] == 'neumann'):
        D[-1, -1] = 0
    return D

def ddxCellGradBC(n, bc):
    """

        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1

        For Cell-Centered **Dirichlet**, use a ghost point::

            (u_1 - u_g)/hf = grad

             u_g       u_1      u_2
              *    |    *   |    *     ...
                   ^
                  u_b

        We know the value at the boundary (u_b)::

            (u_g+u_1)/2 = u_b               (the average)
            u_g = 2*u_b - u_1

            So plug in to gradient:

            (u_1 - (2*u_b - u_1))/hf = grad
            2*(u_1-u_b)/hf = grad

        Separate, because BC are known (and can move to RHS later)::

            ( 2/hf )*u_1 + ( -2/hf )*u_b = grad

                           (   ^   ) JUST RETURN THIS


    """
    bc = checkBC(bc)

    ij   = (np.array([0, n]),np.array([0, 1]))
    vals = np.zeros(2)

    # Set the first side
    if(bc[0] == 'dirichlet'):
        vals[0] = -2
    elif(bc[0] == 'neumann'):
        vals[0] = 0
    # Set the second side
    if(bc[1] == 'dirichlet'):
        vals[1] = 2
    elif(bc[1] == 'neumann'):
        vals[1] = 0
    D = sp.csr_matrix((vals, ij), shape=(n+1,2))
    return D


class DiffOperators(object):
    """
        Class creates the differential operators that you need!
    """
    def __init__(self):
        raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.')

    def faceDiv():
        doc = "Construct divergence operator (face-stg to cell-centres)."

        def fget(self):
            if(self._faceDiv is None):
                # The number of cell centers in each direction
                n = self.vnC
                # Compute faceDivergence operator on faces
                if(self.dim == 1):
                    D = ddx(n[0])
                elif(self.dim == 2):
                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
                    D = sp.hstack((D1, D2), format="csr")
                elif(self.dim == 3):
                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
                    D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
                    D = sp.hstack((D1, D2, D3), format="csr")
                # Compute areas of cell faces & volumes
                S = self.area
                V = self.vol
                self._faceDiv = sdiag(1/V)*D*sdiag(S)

            return self._faceDiv
        return locals()
    _faceDiv = None
    faceDiv = property(**faceDiv())

    def faceDivx():
        doc = "Construct divergence operator in the x component (face-stg to cell-centres)."

        def fget(self):
            if(self._faceDivx is None):
                # The number of cell centers in each direction
                n = self.vnC
                # Compute faceDivergence operator on faces
                if(self.dim == 1):
                    D1 = ddx(n[0])
                elif(self.dim == 2):
                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
                elif(self.dim == 3):
                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
                # Compute areas of cell faces & volumes
                S = self.r(self.area, 'F', 'Fx', 'V')
                V = self.vol
                self._faceDivx = sdiag(1/V)*D1*sdiag(S)

            return self._faceDivx
        return locals()
    _faceDivx = None
    faceDivx = property(**faceDivx())

    def faceDivy():
        doc = "Construct divergence operator in the y component (face-stg to cell-centres)."

        def fget(self):
            if(self.dim < 2): return None
            if(self._faceDivy is None):
                # The number of cell centers in each direction
                n = self.vnC
                # Compute faceDivergence operator on faces
                if(self.dim == 2):
                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
                elif(self.dim == 3):
                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
                # Compute areas of cell faces & volumes
                S = self.r(self.area, 'F', 'Fy', 'V')
                V = self.vol
                self._faceDivy = sdiag(1/V)*D2*sdiag(S)

            return self._faceDivy
        return locals()
    _faceDivy = None
    faceDivy = property(**faceDivy())

    def faceDivz():
        doc = "Construct divergence operator in the z component (face-stg to cell-centres)."

        def fget(self):
            if(self.dim < 3): return None
            if(self._faceDivz is None):
                # The number of cell centers in each direction
                n = self.vnC
                # Compute faceDivergence operator on faces
                D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
                # Compute areas of cell faces & volumes
                S = self.r(self.area, 'F', 'Fz', 'V')
                V = self.vol
                self._faceDivz = sdiag(1/V)*D3*sdiag(S)

            return self._faceDivz
        return locals()
    _faceDivz = None
    faceDivz = property(**faceDivz())

    def nodalGrad():
        doc = "Construct gradient operator (nodes to edges)."

        def fget(self):
            if(self._nodalGrad is None):
                # The number of cell centers in each direction
                n = self.vnC
                # Compute divergence operator on faces
                if(self.dim == 1):
                    G = ddx(n[0])
                elif(self.dim == 2):
                    D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
                    D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
                    G = sp.vstack((D1, D2), format="csr")
                elif(self.dim == 3):
                    D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
                    D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
                    D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
                    G = sp.vstack((D1, D2, D3), format="csr")
                # Compute lengths of cell edges
                L = self.edge
                self._nodalGrad = sdiag(1/L)*G
            return self._nodalGrad
        return locals()
    _nodalGrad = None
    nodalGrad = property(**nodalGrad())

    def nodalLaplacian():
        doc = "Construct laplacian operator (nodes to edges)."

        def fget(self):
            if(self._nodalLaplacian is None):
                print 'Warning: Laplacian has not been tested rigorously.'
                # The number of cell centers in each direction
                n = self.vnC
                # Compute divergence operator on faces
                if(self.dim == 1):
                    D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
                    L  = - D1.T*D1
                elif(self.dim == 2):
                    D1 = sdiag(1./self.hx) * ddx(n[0])
                    D2 = sdiag(1./self.hy) * ddx(n[1])
                    L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
                    L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
                    L  = L1 + L2
                elif(self.dim == 3):
                    D1 = sdiag(1./self.hx) * ddx(n[0])
                    D2 = sdiag(1./self.hy) * ddx(n[1])
                    D3 = sdiag(1./self.hz) * ddx(n[2])
                    L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
                    L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
                    L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
                    L  = L1 + L2 + L3
                self._nodalLaplacian = L
            return self._nodalLaplacian
        return locals()
    _nodalLaplacian = None
    nodalLaplacian = property(**nodalLaplacian())

    def setCellGradBC(self, BC):
        """
        Function that sets the boundary conditions for cell-centred derivative operators.

        Examples::

            BC = 'neumann'                                            # Neumann in all directions
            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
                                                                      #     Dirichlet else

        """
        if(type(BC) is str):
            BC = [BC]*self.dim
        if(type(BC) is list):
            assert len(BC) == self.dim, 'BC list must be the size of your mesh'
        else:
            raise Exception("BC must be a str or a list.")

        for i, bc_i in enumerate(BC):
            BC[i] = checkBC(bc_i)

        # ensure we create a new gradient next time we call it
        self._cellGrad        = None
        self._cellGradBC      = None
        self._cellGradBC_list = BC
        return BC
    _cellGradBC_list = 'neumann'

    def cellGrad():
        doc = "The cell centered Gradient, takes you to cell faces."

        def fget(self):
            if(self._cellGrad is None):
                BC = self.setCellGradBC(self._cellGradBC_list)
                n = self.vnC
                if(self.dim == 1):
                    G = ddxCellGrad(n[0], BC[0])
                elif(self.dim == 2):
                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
                    G = sp.vstack((G1, G2), format="csr")
                elif(self.dim == 3):
                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
                    G = sp.vstack((G1, G2, G3), format="csr")
                # Compute areas of cell faces & volumes
                S = self.area
                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
                self._cellGrad = sdiag(S/V)*G
            return self._cellGrad
        return locals()
    _cellGrad = None
    cellGrad = property(**cellGrad())

    def cellGradBC():
        doc = "The cell centered Gradient boundary condition matrix"

        def fget(self):
            if(self._cellGradBC is None):
                BC = self.setCellGradBC(self._cellGradBC_list)
                n = self.vnC
                if(self.dim == 1):
                    G = ddxCellGradBC(n[0], BC[0])
                elif(self.dim == 2):
                    G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                    G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                    G = sp.block_diag((G1, G2), format="csr")
                elif(self.dim == 3):
                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                    G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                    G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
                    G = sp.block_diag((G1, G2, G3), format="csr")
                # Compute areas of cell faces & volumes
                S = self.area
                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
                self._cellGradBC = sdiag(S/V)*G
            return self._cellGradBC
        return locals()
    _cellGradBC = None
    cellGradBC = property(**cellGradBC())

    def cellGradx():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."

        def fget(self):
            if getattr(self, '_cellGradx', None) is None:
                BC = ['neumann', 'neumann']
                n = self.vnC
                if(self.dim == 1):
                    G1 = ddxCellGrad(n[0], BC)
                elif(self.dim == 2):
                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
                elif(self.dim == 3):
                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fx', 'V')
                self._cellGradx = sdiag(L)*G1
            return self._cellGradx
        return locals()
    cellGradx = property(**cellGradx())

    def cellGrady():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
        def fget(self):
            if self.dim < 2: return None
            if getattr(self, '_cellGrady', None) is None:
                BC = ['neumann', 'neumann']
                n = self.vnC
                if(self.dim == 2):
                    G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
                elif(self.dim == 3):
                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fy', 'V')
                self._cellGrady = sdiag(L)*G2
            return self._cellGrady
        return locals()
    cellGrady = property(**cellGrady())

    def cellGradz():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
        def fget(self):
            if self.dim < 3: return None
            if getattr(self, '_cellGradz', None) is None:
                BC = ['neumann', 'neumann']
                n = self.vnC
                G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fz', 'V')
                self._cellGradz = sdiag(L)*G3
            return self._cellGradz
        return locals()
    cellGradz = property(**cellGradz())

    def edgeCurl():
        doc = "Construct the 3D curl operator."

        def fget(self):
            if(self._edgeCurl is None):
                assert self.dim > 2, "Edge Curl only programed for 3D."
                # The number of cell centers in each direction
                n1 = self.nCx
                n2 = self.nCy
                n3 = self.nCz

                # Compute lengths of cell edges
                L = self.edge

                # Compute areas of cell faces
                S = self.area

                # Compute divergence operator on faces
                d1 = ddx(n1)
                d2 = ddx(n2)
                d3 = ddx(n3)

                D32 = kron3(d3, speye(n2), speye(n1+1))
                D23 = kron3(speye(n3), d2, speye(n1+1))
                D31 = kron3(d3, speye(n2+1), speye(n1))
                D13 = kron3(speye(n3), speye(n2+1), d1)
                D21 = kron3(speye(n3+1), d2, speye(n1))
                D12 = kron3(speye(n3+1), speye(n2), d1)

                O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])

                C = sp.vstack((sp.hstack((O1, -D32, D23)),
                               sp.hstack((D31, O2, -D13)),
                               sp.hstack((-D21, D12, O3))), format="csr")

                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
            return self._edgeCurl
        return locals()
    _edgeCurl = None
    edgeCurl = property(**edgeCurl())

    def getBCProjWF(self, BC, discretization='CC'):
        """

        The weak form boundary condition projection matrices.

        Examples::

            BC = 'neumann'                                            # Neumann in all directions
            BC = ['neumann', 'dirichlet', 'neumann']                  # 3D, Dirichlet in y Neumann else
            BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
                                                                      #     Dirichlet else

        """

        if discretization is not 'CC':
            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')

        if(type(BC) is str):
            BC = [BC for _ in self.vnC]  # Repeat the str self.dim times
        elif(type(BC) is list):
            assert len(BC) == self.dim, 'BC list must be the size of your mesh'
        else:
            raise Exception("BC must be a str or a list.")

        for i, bc_i in enumerate(BC):
            BC[i] = checkBC(bc_i)


        def projDirichlet(n, bc):
            bc = checkBC(bc)
            ij = ([0,n], [0,1])
            vals = [0,0]
            if(bc[0] == 'dirichlet'):
                vals[0] = -1
            if(bc[1] == 'dirichlet'):
                vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(n+1,2))

        def projNeumannIn(n, bc):
            bc = checkBC(bc)
            P = sp.identity(n+1).tocsr()
            if(bc[0] == 'neumann'):
                P = P[1:,:]
            if(bc[1] == 'neumann'):
                P = P[:-1,:]
            return P

        def projNeumannOut(n, bc):
            bc = checkBC(bc)
            ij   = ([0, 1],[0, n])
            vals = [0,0]
            if(bc[0] == 'neumann'):
                vals[0] = 1
            if(bc[1] == 'neumann'):
                vals[1] = 1
            return sp.csr_matrix((vals, ij), shape=(2,n+1))

        n = self.vnC
        indF = self.faceBoundaryInd
        if(self.dim == 1):
            Pbc = projDirichlet(n[0], BC[0])
            indF = indF[0] | indF[1]
            Pbc = Pbc*sdiag(self.area[indF])

            Pin = projNeumannIn(n[0], BC[0])

            Pout = projNeumannOut(n[0], BC[0])
        elif(self.dim == 2):
            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
            Pbc = Pbc*sdiag(self.area[indF])

            P1 = sp.kron(speye(n[1]), projNeumannIn(n[0], BC[0]))
            P2 = sp.kron(projNeumannIn(n[1], BC[1]), speye(n[0]))
            Pin = sp.block_diag((P1, P2), format="csr")

            P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
            P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
            Pout = sp.block_diag((P1, P2), format="csr")
        elif(self.dim == 3):
            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
            Pbc = Pbc*sdiag(self.area[indF])

            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
            P2 = kron3(speye(n[2]), projNeumannIn(n[1], BC[1]), speye(n[0]))
            P3 = kron3(projNeumannIn(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pin = sp.block_diag((P1, P2, P3), format="csr")

            P1 = kron3(speye(n[2]), speye(n[1]), projNeumannOut(n[0], BC[0]))
            P2 = kron3(speye(n[2]), projNeumannOut(n[1], BC[1]), speye(n[0]))
            P3 = kron3(projNeumannOut(n[2], BC[2]), speye(n[1]), speye(n[0]))
            Pout = sp.block_diag((P1, P2, P3), format="csr")

        return Pbc, Pin, Pout


    # --------------- Averaging ---------------------

    @property
    def aveF2CC(self):
        "Construct the averaging operator on cell faces to cell centers."
        if getattr(self, '_aveF2CC', None) is None:
            n = self.vnC
            if(self.dim == 1):
                self._aveF2CC = av(n[0])
            elif(self.dim == 2):
                self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
                                                 sp.kron(av(n[1]), speye(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
                                                   kron3(speye(n[2]), av(n[1]), speye(n[0])),
                                                   kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
        return self._aveF2CC


    @property
    def aveF2CCV(self):
        "Construct the averaging operator on cell faces to cell centers."
        if getattr(self, '_aveF2CCV', None) is None:
            n = self.vnC
            if(self.dim == 1):
                self._aveF2CCV = av(n[0])
            elif(self.dim == 2):
                self._aveF2CCV = sp.block_diag((sp.kron(speye(n[1]), av(n[0])),
                                               sp.kron(av(n[1]), speye(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveF2CCV = sp.block_diag((kron3(speye(n[2]), speye(n[1]), av(n[0])),
                                               kron3(speye(n[2]), av(n[1]), speye(n[0])),
                                               kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
        return self._aveF2CCV

    @property
    def aveCC2F(self):
        "Construct the averaging operator on cell cell centers to faces."
        if getattr(self, '_aveCC2F', None) is None:
            n = self.vnC
            if(self.dim == 1):
                self._aveCC2F = avExtrap(n[0])
            elif(self.dim == 2):
                self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
                                           sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
                                           kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
                                           kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
        return self._aveCC2F

    @property
    def aveE2CC(self):
        "Construct the averaging operator on cell edges to cell centers."
        if getattr(self, '_aveE2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
            elif(self.dim == 2):
                self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
                                               sp.kron(speye(n[1]), av(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
                                                  kron3(av(n[2]), speye(n[1]), av(n[0])),
                                                  kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
        return self._aveE2CC

    @property
    def aveE2CCV(self):
        "Construct the averaging operator on cell edges to cell centers."
        if getattr(self, '_aveE2CCV', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
            elif(self.dim == 2):
                self._aveE2CCV = sp.block_diag((sp.kron(av(n[1]), speye(n[0])),
                                                sp.kron(speye(n[1]), av(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveE2CCV = sp.block_diag((kron3(av(n[2]), av(n[1]), speye(n[0])),
                                                kron3(av(n[2]), speye(n[1]), av(n[0])),
                                                kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
        return self._aveE2CCV

    @property
    def aveN2CC(self):
        "Construct the averaging operator on cell nodes to cell centers."
        if getattr(self, '_aveN2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveN2CC = av(n[0])
            elif(self.dim == 2):
                self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
            elif(self.dim == 3):
                self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
        return self._aveN2CC

    @property
    def aveN2E(self):
        "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."

        if getattr(self, '_aveN2E', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveN2E = av(n[0])
            elif(self.dim == 2):
                self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
                                          sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
            elif(self.dim == 3):
                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
                                          kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
                                          kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
        return self._aveN2E

    @property
    def aveN2F(self):
        "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
        if getattr(self, '_aveN2F', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveN2F = av(n[0])
            elif(self.dim == 2):
                self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
                                          sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
            elif(self.dim == 3):
                self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
                                          kron3(av(n[2]), speye(n[1]+1), av(n[0])),
                                          kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
        return self._aveN2F