#include "macrodef.h"
#include "share_func.h"

/*
 * kodis_sh — Kreiss-Oliger dissipation on shell patches.
 * Same stencil coefficients as Cartesian kodis. Uses symmetry_stbd.
 */
extern "C" void kodis_sh_(const int ex[3],
              const double *X, const double *Y, const double *Z,
              const double *f, double *f_rhs,
              const double SoAi[2],
              int Symmetry, double eps, int sst)
{
    (void)sst;
    const double ZEO=0.0;
    const int ex1=ex[0], ex2=ex[1], ex3=ex[2];
    const double dX=X[1]-X[0], dY=Y[1]-Y[0], dZ=Z[1]-Z[0];
    const int imaxF=ex1, jmaxF=ex2, kmaxF=ex3;
    const double SoA[2]={SoAi[0],SoAi[1]};

#if (ghost_width == 2)
    {
        const int ord=2, r=2;
        const double cof=16.0, F4=4.0, F6=6.0;
        const int NO_SYMM=0, OCTANT=2;
        int iminF=1,jminF=1,kminF=1;
        if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-1;
        if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-1;
        if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-1;

        const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
        double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
        symmetry_stbd(ord,ex,f,fh,SoA);

        const int i0_lo=(iminF+1>0)?iminF+1:0,j0_lo=(jminF+1>0)?jminF+1:0,k0_lo=0;
        const int i0_hi=imaxF-3,j0_hi=jmaxF-3,k0_hi=kmaxF-3;
        if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
            for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
            for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
            for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
                const double Dx=((fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])-F4*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
                const double Dy=((fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])-F4*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
                const double Dz=((fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])-F4*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
                f_rhs[p]-=(eps/cof)*(Dx+Dy+Dz);
            }}}
        }
        free(fh);return;
    }
#elif (ghost_width == 3)
    {
        const int ord=3, r=3;
        const double cof=64.0,SIX=6.0,FIT=15.0,TWT=20.0;
        const int NO_SYMM=0,OCTANT=2;
        int iminF=1,jminF=1,kminF=1;
        if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-2;
        if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-2;
        if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-2;

        const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
        double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
        symmetry_stbd(ord,ex,f,fh,SoA);

        const int i0_lo=(iminF+2>0)?iminF+2:0,j0_lo=(jminF+2>0)?jminF+2:0,k0_lo=0;
        const int i0_hi=imaxF-4,j0_hi=jmaxF-4,k0_hi=kmaxF-4;
        if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
            for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
            for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
            for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
                const double Dx=((fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])-SIX*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])+FIT*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
                const double Dy=((fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])-SIX*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])+FIT*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
                const double Dz=((fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])-SIX*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])+FIT*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
                f_rhs[p]+=(eps/cof)*(Dx+Dy+Dz);
            }}}
        }
        free(fh);return;
    }
#elif (ghost_width == 4)
    {
        const int ord=4, r=4;
        const double cof=256.0,F8=8.0,F28=28.0,F56=56.0,F70=70.0;
        const int NO_SYMM=0,OCTANT=2;
        int iminF=1,jminF=1,kminF=1;
        if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-3;
        if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-3;
        if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-3;

        const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
        double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
        symmetry_stbd(ord,ex,f,fh,SoA);

        const int i0_lo=(iminF+3>0)?iminF+3:0,j0_lo=(jminF+3>0)?jminF+3:0,k0_lo=0;
        const int i0_hi=imaxF-5,j0_hi=jmaxF-5,k0_hi=kmaxF-5;
        if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
            for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
            for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
            for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
                const double Dx=((fh[idx_fh_stbd(iF-4,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+4,jF,kF,ord,ex)])-F8*(fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])+F28*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])-F56*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
                const double Dy=((fh[idx_fh_stbd(iF,jF-4,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+4,kF,ord,ex)])-F8*(fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])+F28*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])-F56*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
                const double Dz=((fh[idx_fh_stbd(iF,jF,kF-4,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+4,ord,ex)])-F8*(fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])+F28*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])-F56*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
                f_rhs[p]-=(eps/cof)*(Dx+Dy+Dz);
            }}}
        }
        free(fh);return;
    }
#elif (ghost_width == 5)
    {
        const int ord=5, r=5;
        const double cof=1024.0,F10=10.0,F45k=45.0,F120=120.0,F210=210.0,F252=252.0;
        const int NO_SYMM=0,OCTANT=2;
        int iminF=1,jminF=1,kminF=1;
        if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-4;
        if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-4;
        if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-4;

        const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
        double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
        symmetry_stbd(ord,ex,f,fh,SoA);

        const int i0_lo=(iminF+4>0)?iminF+4:0,j0_lo=(jminF+4>0)?jminF+4:0,k0_lo=0;
        const int i0_hi=imaxF-6,j0_hi=jmaxF-6,k0_hi=kmaxF-6;
        if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
            for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
            for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
            for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
                const double Dx=((fh[idx_fh_stbd(iF-5,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+5,jF,kF,ord,ex)])-F10*(fh[idx_fh_stbd(iF-4,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+4,jF,kF,ord,ex)])+F45k*(fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])-F120*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])+F210*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
                const double Dy=((fh[idx_fh_stbd(iF,jF-5,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+5,kF,ord,ex)])-F10*(fh[idx_fh_stbd(iF,jF-4,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+4,kF,ord,ex)])+F45k*(fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])-F120*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])+F210*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
                const double Dz=((fh[idx_fh_stbd(iF,jF,kF-5,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+5,ord,ex)])-F10*(fh[idx_fh_stbd(iF,jF,kF-4,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+4,ord,ex)])+F45k*(fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])-F120*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])+F210*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
                f_rhs[p]+=(eps/cof)*(Dx+Dy+Dz);
            }}}
        }
        free(fh);return;
    }
#else
#error "kodiss_sh_c.C: unsupported ghost_width"
#endif
}