diff --git a/fdelmodc3D/Makefile b/fdelmodc3D/Makefile index e3f21d2bd2a8f7636a8bde0f7bf7ce8784177db4..40cd0621e9f9f25afb3f1c51033fa323fc8c9770 100644 --- a/fdelmodc3D/Makefile +++ b/fdelmodc3D/Makefile @@ -8,6 +8,7 @@ ALLINC = -I. LIBS += -L$L -lgenfft -lm $(LIBSM) #LIBS += -L$L -lgenfft -lm -lc #OPTC = -g -Wall -fsignaling-nans -O0 -fopenmp +#OPTC += -g #OPTC += -fopenmp -Waddress #OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC)) #PGI options for compiler feedback diff --git a/fdelmodc3D/acoustic2_3D.c b/fdelmodc3D/acoustic2_3D.c index 3793a59c72b62d68547ee636df9700f1da082bcd..6767c9b76f9c4a52a3316eca6b73445f506da8cf 100644 --- a/fdelmodc3D/acoustic2_3D.c +++ b/fdelmodc3D/acoustic2_3D.c @@ -7,7 +7,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, @@ -19,17 +19,17 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose) + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose) { /********************************************************************* @@ -123,7 +123,7 @@ long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, #pragma ivdep for (iz=ioXz; iz<ieXz; iz++) { vx[iy*n2*n1+ix*n1+iz] -= - rox[iy*n2*n1+ix*n1+iz]*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]); + rox[iy][ix][iz]*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]); } } } @@ -135,7 +135,7 @@ long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, #pragma ivdep for (iz=ioYz; iz<ieYz; iz++) { vy[iy*n2*n1+ix*n1+iz] -= - roy[iy*n2*n1+ix*n1+iz]*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]); + roy[iy][ix][iz]*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]); } } } @@ -147,7 +147,7 @@ long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, #pragma ivdep for (iz=ioZz; iz<ieZz; iz++) { vz[iy*n2*n1+ix*n1+iz] -= - roz[iy*n2*n1+ix*n1+iz]*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]); + roz[iy][ix][iz]*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]); } } } @@ -175,7 +175,7 @@ long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioPx; ix<iePx; ix++) { #pragma ivdep for (iz=ioPz; iz<iePz; iz++) { - p[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*( + p[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*( (vx[iy*n2*n1+(ix+1)*n1+iz] - vx[iy*n2*n1+ix*n1+iz]) + (vy[(iy+1)*n2*n1+ix*n1+iz] - vy[iy*n2*n1+ix*n1+iz]) + (vz[iy*n2*n1+ix*n1+iz+1] - vz[iy*n2*n1+ix*n1+iz])); diff --git a/fdelmodc3D/acoustic4_3D.c b/fdelmodc3D/acoustic4_3D.c index 36316d652c7b0bbb1c4006368d18d0ad13003a1b..6029d539497639edf40968dcc59910352a31f541 100644 --- a/fdelmodc3D/acoustic4_3D.c +++ b/fdelmodc3D/acoustic4_3D.c @@ -8,7 +8,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose); @@ -17,15 +17,15 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose); long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, - float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *mul, long itime, long verbose); + float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, - float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *mul, long itime, long verbose); + float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***mul, long itime, long verbose); long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, - float *vx, float *vy, float *vz, float *p, float *rox, float *roy, float *roz, float *l2m, long verbose) + float *vx, float *vy, float *vz, float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose) { /********************************************************************* COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: @@ -89,7 +89,7 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo for (ix=mod.ioXx; ix<mod.ieXx; ix++) { #pragma ivdep for (iz=mod.ioXz; iz<mod.ieXz; iz++) { - vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*( + vx[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]) + c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz])); } @@ -102,7 +102,7 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo for (ix=mod.ioYx; ix<mod.ieYx; ix++) { #pragma ivdep for (iz=mod.ioYz; iz<mod.ieYz; iz++) { - vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*( + vy[iy*n2*n1+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]) + c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*n1+iz])); } @@ -115,7 +115,7 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo for (ix=mod.ioZx; ix<mod.ieZx; ix++) { #pragma ivdep for (iz=mod.ioZz; iz<mod.ieZz; iz++) { - vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*( + vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]) + c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2])); } @@ -148,7 +148,7 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo for (iy=mod.ioPy; iy<mod.iePy; iy++) { #pragma ivdep for (iz=mod.ioPz; iz<mod.iePz; iz++) { - p[iy*n1*n2+ix*n1+iz] -= l2m[iy*n1*n2+ix*n1+iz]*( + p[iy*n1*n2+ix*n1+iz] -= l2m[iy][ix][iz]*( c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) + c2*(vx[iy*n1*n2+(ix+2)*n1+iz] - vx[iy*n1*n2+(ix-1)*n1+iz]) + c1*(vy[(iy+1)*n1*n2+ix*n1+iz] - vy[iy*n1*n2+ix*n1+iz]) + diff --git a/fdelmodc3D/acoustic4_qr_3D.c b/fdelmodc3D/acoustic4_qr_3D.c index 08d8d1b4910d9040fdf756c84355ee9d326823eb..5a87306257273549b8c2ee61b5b306b07c9d2660 100644 --- a/fdelmodc3D/acoustic4_qr_3D.c +++ b/fdelmodc3D/acoustic4_qr_3D.c @@ -9,7 +9,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, @@ -21,17 +21,17 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose) + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose) { /********************************************************************* COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: @@ -211,7 +211,7 @@ long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, } #pragma ivdep for (iz=ioXz; iz<ieXz; iz++) { - vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*( + vx[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]) + c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz])); } @@ -232,7 +232,7 @@ long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, } #pragma ivdep for (iz=ioYz; iz<ieYz; iz++) { - vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*( + vy[iy*n2*n1+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]) + c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*n1+iz])); } @@ -253,7 +253,7 @@ long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, } #pragma ivdep for (iz=ioZz; iz<ieZz; iz++) { - vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*( + vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]) + c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2])); } @@ -283,7 +283,7 @@ long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, } #pragma ivdep for (iz=ioPz; iz<iePz; iz++) { - p[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*( + p[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*( c1*(vx[iy*n2*n1+(ix+1)*n1+iz] - vx[iy*n2*n1+ix*n1+iz]) + c2*(vx[iy*n2*n1+(ix+2)*n1+iz] - vx[iy*n2*n1+(ix-1)*n1+iz]) + c1*(vy[(iy+1)*n2*n1+ix*n1+iz] - vy[iy*n2*n1+ix*n1+iz]) + diff --git a/fdelmodc3D/acoustic6_3D.c b/fdelmodc3D/acoustic6_3D.c index 1ae3f359516d271ef583f621db30c5db67a97e30..a28ed894429edbba0db7865bb428c876473b0f70 100644 --- a/fdelmodc3D/acoustic6_3D.c +++ b/fdelmodc3D/acoustic6_3D.c @@ -7,7 +7,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, @@ -19,17 +19,17 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose) + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose) { /********************************************************************* COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: @@ -125,7 +125,7 @@ long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioXx; ix<ieXx; ix++) { #pragma ivdep for (iz=ioXz; iz<ieXz; iz++) { - vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*( + vx[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]) + c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz]) + c3*(p[iy*n2*n1+(ix+2)*n1+iz] - p[iy*n2*n1+(ix-3)*n1+iz])); @@ -139,7 +139,7 @@ long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioYx; ix<ieYx; ix++) { #pragma ivdep for (iz=ioYz; iz<ieYz; iz++) { - vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*( + vy[iy*n2*n1+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]) + c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*n1+iz]) + c3*(p[(iy+2)*n2*n1+ix*n1+iz] - p[(iy-3)*n2*n1+ix*n1+iz])); @@ -153,7 +153,7 @@ long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioZx; ix<ieZx; ix++) { #pragma ivdep for (iz=ioZz; iz<ieZz; iz++) { - vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*( + vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]) + c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2]) + c3*(p[iy*n2*n1+ix*n1+iz+2] - p[iy*n2*n1+ix*n1+iz-3])); @@ -184,7 +184,7 @@ long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioPx; ix<iePx; ix++) { #pragma ivdep for (iz=ioPz; iz<iePz; iz++) { - p[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*( + p[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*( c1*(vx[iy*n2*n1+(ix+1)*n1+iz] - vx[iy*n2*n1+ix*n1+iz]) + c2*(vx[iy*n2*n1+(ix+2)*n1+iz] - vx[iy*n2*n1+(ix-1)*n1+iz]) + c3*(vx[iy*n2*n1+(ix+3)*n1+iz] - vx[iy*n2*n1+(ix-2)*n1+iz]) + diff --git a/fdelmodc3D/acousticSH4_3D.c b/fdelmodc3D/acousticSH4_3D.c index dc7c5648e0ce178afdb04a14bed3c4f985398071..3353892584e2a6fd9ac7636ec3e0d7442b2261f7 100644 --- a/fdelmodc3D/acousticSH4_3D.c +++ b/fdelmodc3D/acousticSH4_3D.c @@ -7,7 +7,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, @@ -19,17 +19,17 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose); long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *tx, float *ty, float *tz, - float *vz, float *rox, float *roy, float *roz, float *mul, long verbose) + float *vz, float ***rox, float ***roy, float ***roz, float ***mul, long verbose) { /********************************************************************* COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: @@ -123,7 +123,7 @@ long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioXx; ix<ieXx; ix++) { #pragma ivdep for (iz=ioXz; iz<ieXz; iz++) { - tx[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*( + tx[iy*n2*n1+ix*n1+iz] -= mul[iy][ix][iz]*( c1*(vz[iy*n2*n1+ix*n1+iz] - vz[iy*n2*n1+(ix-1)*n1+iz]) + c2*(vz[iy*n2*n1+(ix+1)*n1+iz] - vz[iy*n2*n1+(ix-2)*n1+iz])); } @@ -136,7 +136,7 @@ long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioYx; ix<ieYx; ix++) { #pragma ivdep for (iz=ioYz; iz<ieYz; iz++) { - ty[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*( + ty[iy*n2*n1+ix*n1+iz] -= mul[iy][ix][iz]*( c1*(vz[iy*n2*n1+ix*n1+iz] - vz[(iy-1)*n2*n1+ix*n1+iz]) + c2*(vz[(iy+1)*n2*n1+ix*n1+iz] - vz[(iy-2)*n2*n1+ix*n1+iz])); } @@ -149,7 +149,7 @@ long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioZx; ix<ieZx; ix++) { #pragma ivdep for (iz=ioZz; iz<ieZz; iz++) { - tz[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*( + tz[iy*n2*n1+ix*n1+iz] -= mul[iy][ix][iz]*( c1*(vz[iy*n2*n1+ix*n1+iz] - vz[iy*n2*n1+ix*n1+iz-1]) + c2*(vz[iy*n2*n1+ix*n1+iz+1] - vz[iy*n2*n1+ix*n1+iz-2])); } @@ -171,7 +171,7 @@ long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, for (ix=ioPx; ix<iePx; ix++) { #pragma ivdep for (iz=ioPz; iz<iePz; iz++) { - vz[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*( + vz[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tx[iy*n2*n1+(ix+1)*n1+iz] - tx[iy*n2*n1+ix*n1+iz]) + c2*(tx[iy*n2*n1+(ix+2)*n1+iz] - tx[iy*n2*n1+(ix-1)*n1+iz]) + c1*(ty[iy*n2*n1+(ix+1)*n1+iz] - ty[iy*n2*n1+ix*n1+iz]) + diff --git a/fdelmodc3D/alloc3D.c b/fdelmodc3D/alloc3D.c index d82748fabacda9dd5605dae2aa35815a77d15997..139d1b126485f0e137a2c9f756dadcc96d364dd4 100644 --- a/fdelmodc3D/alloc3D.c +++ b/fdelmodc3D/alloc3D.c @@ -5,10 +5,10 @@ /* All rights reserved. */ /* allocate a 3-d array of floats */ -float ***alloc3float(modPar mod) +void ***alloc3float(modPar mod) { size_t i3,i2,n1,n2,n3,size; - float ***p; + void ***p; n1 = mod.naz; n2 = mod.nax; @@ -16,9 +16,9 @@ float ***alloc3float(modPar mod) size = sizeof(float); /* allocate pointer arrays of n3 and n2 */ - if ((p=(float***)malloc(n3*sizeof(float**)))==NULL) + if ((p=(void***)malloc(n3*sizeof(void**)))==NULL) return NULL; - if ((p[0]=(float**)malloc(n3*n2*sizeof(float*)))==NULL) { + if ((p[0]=(void**)malloc(n3*n2*sizeof(void*)))==NULL) { free(p); return NULL; } @@ -26,31 +26,31 @@ float ***alloc3float(modPar mod) /* allocate data size and assign pointers to n1, n2 and n3 */ if (mod.nfx==1 && mod.nfy==1 ) { // 1D model fprintf(stderr,"Allocating 1D model \n"); - if ((p[0][0]=(float*)malloc(n1*size))==NULL) { + if ((p[0][0]=(void*)malloc(n1*size))==NULL) { free(p); return NULL; } for (i3=0; i3<n3; i3++) { - p[i3] = p[0]; + p[i3] = p[0]+n2*i3; for (i2=0; i2<n2; i2++) - p[i3][i2] = (float*)p[0][0]; + p[i3][i2] = (void*)p[0][0]; } } else if (mod.nfy==1) { // 2D model - fprintf(stderr,"Allocating 2D model \n"); - if ((p[0][0]=(float*)malloc(n2*n1*size))==NULL) { + fprintf(stderr,"Allocating 2D model ny=%d with nx=%d and nz=%d\n",n3,n2,n1); + if ((p[0][0]=(void*)malloc(n2*n1*size))==NULL) { free(p); return NULL; } for (i3=0; i3<n3; i3++) { - p[i3] = p[0]; + p[i3] = p[0]+n2*i3; for (i2=0; i2<n2; i2++) - p[i3][i2] = (float*)p[0]+n1*i2; + p[i3][i2] = (char*)p[0][0]+size*n1*i2; } } else { // 3D model fprintf(stderr,"Allocating 3D model \n"); - if ((p[0][0]=(float*)malloc(n3*n2*n1*size))==NULL) { + if ((p[0][0]=(void*)malloc(n3*n2*n1*size))==NULL) { free(p[0]); free(p); return NULL; @@ -58,13 +58,13 @@ float ***alloc3float(modPar mod) for (i3=0; i3<n3; i3++) { p[i3] = p[0]+n2*i3; for (i2=0; i2<n2; i2++) - p[i3][i2] = (float*)p[0][0]+n1*(i2+n2*i3); + p[i3][i2] = (char*)p[0][0]+size*n1*(i2+n2*i3); } } return p; } -void free3float(float ***p) +void free3float(void ***p) { free(p[0][0]); free(p[0]); diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c index 0a7d1f1165a4eb2d7c7f4c3bdcaa374fd205c1c0..51c21f35b62421ced0b6d98d15da04dec6dee9de 100644 --- a/fdelmodc3D/applySource3D.c +++ b/fdelmodc3D/applySource3D.c @@ -22,7 +22,7 @@ void vmess(char *fmt, ...); **********************************************************************/ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose) + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose) { long is0, ibndz, ibndy, ibndx; long isrc, ix, iy, iz, n1, n2; @@ -147,7 +147,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l /* in older version added factor 2.0 to be compliant with defined Green's functions in Marchenko algorithm */ /* this is now set to 1.0 */ - src_ampl *= (1.0/(mod.dx*mod.dx))*l2m[iy*n1*n2+ix*n1+iz]; + src_ampl *= (1.0/(mod.dx*mod.dx))*l2m[iy][ix][iz]; if (verbose>5) { vmess("Source %li at grid [ix=%li,iy=%li,iz=%li] at itime %li has value %e",isrc, ix, iy, iz, itime, src_ampl); @@ -156,12 +156,12 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l /* Force source */ if (src.type == 6) { - vx[iy*n1*n2+ix*n1+iz] += src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]); + vx[iy*n1*n2+ix*n1+iz] += src_ampl*rox[iy][ix][iz]/(l2m[iy][ix][iz]); /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */ //vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]); } else if (src.type == 7) { - vz[iy*n1*n2+ix*n1+iz] += src_ampl*roz[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]); + vz[iy*n1*n2+ix*n1+iz] += src_ampl*roz[iy][ix][iz]/(l2m[iy][ix][iz]); /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */ /* stable implementation changes amplitude and more work is needed */ //vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]); @@ -274,7 +274,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l * Curl S-pot = CURL(F) = dF_x/dz - dF_z/dx ***********************************************************************/ else if(src.type == 5) { - src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]); + src_ampl = src_ampl*rox[iy][ix][iz]/(l2m[iy][ix][iz]); if (src.orient == 3) src_ampl = -src_ampl; /* first order derivatives */ vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx; @@ -322,7 +322,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l * Divergence P-pot = DIV(F) = dF_x/dx + dF_z/dz ***********************************************************************/ else if(src.type == 8) { - src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]); + src_ampl = src_ampl*rox[iy][ix][iz]/(l2m[iy][ix][iz]); if (src.orient == 3) src_ampl = -src_ampl; vx[iy*n1*n2+(ix+1)*n1+iz] += src_ampl*sdx; vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx; @@ -369,4 +369,4 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l } /* loop over isrc */ return 0; -} \ No newline at end of file +} diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c index 8281c66e6db3dc7b1e34c946381aeec9c3711839..12905a61a4de54fd597d9f7bf6da00b89f4dc58a 100644 --- a/fdelmodc3D/boundaries3D.c +++ b/fdelmodc3D/boundaries3D.c @@ -8,7 +8,7 @@ void vmess(char *fmt, ...); long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose) + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose) { /********************************************************************* Joeri's original filling order @@ -223,7 +223,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -244,7 +244,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -264,7 +264,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -285,7 +285,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -314,7 +314,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -333,7 +333,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -354,7 +354,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -379,7 +379,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -401,7 +401,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -427,7 +427,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -448,7 +448,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -466,7 +466,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -488,7 +488,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -517,7 +517,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -536,7 +536,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -556,7 +556,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -579,7 +579,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -599,7 +599,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -625,7 +625,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -644,7 +644,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -662,7 +662,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -681,7 +681,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -710,7 +710,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -728,7 +728,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -747,7 +747,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -770,7 +770,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -790,7 +790,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -818,7 +818,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -843,7 +843,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -868,7 +868,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -894,7 +894,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -928,7 +928,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -952,7 +952,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -978,7 +978,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1008,7 +1008,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1035,7 +1035,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1066,7 +1066,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1092,7 +1092,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1115,7 +1115,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1142,7 +1142,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1175,7 +1175,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1198,7 +1198,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1222,7 +1222,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1249,7 +1249,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1273,7 +1273,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -1303,7 +1303,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1326,7 +1326,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1348,7 +1348,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1371,7 +1371,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1404,7 +1404,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1426,7 +1426,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1449,7 +1449,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1476,7 +1476,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1500,7 +1500,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -1539,7 +1539,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1558,7 +1558,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1576,7 +1576,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1595,7 +1595,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1624,7 +1624,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1642,7 +1642,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1661,7 +1661,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1684,7 +1684,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1704,7 +1704,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -1730,7 +1730,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1749,7 +1749,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1767,7 +1767,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1786,7 +1786,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1814,7 +1814,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1832,7 +1832,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1851,7 +1851,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1874,7 +1874,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1894,7 +1894,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -1919,7 +1919,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib]; @@ -1939,7 +1939,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -1958,7 +1958,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -1980,7 +1980,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2010,7 +2010,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2028,7 +2028,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2047,7 +2047,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2070,7 +2070,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2090,7 +2090,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2118,7 +2118,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2141,7 +2141,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2163,7 +2163,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2186,7 +2186,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2219,7 +2219,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2241,7 +2241,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2264,7 +2264,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2291,7 +2291,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2315,7 +2315,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2345,7 +2345,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2368,7 +2368,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2390,7 +2390,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2413,7 +2413,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2445,7 +2445,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2467,7 +2467,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2490,7 +2490,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2517,7 +2517,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2541,7 +2541,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -2570,7 +2570,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2594,7 +2594,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2617,7 +2617,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2643,7 +2643,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2677,7 +2677,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2699,7 +2699,7 @@ MID left mid mid for (iy=iyo; iy<iye; iy++) { #pragma ivdep for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2722,7 +2722,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2749,7 +2749,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2773,7 +2773,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -2815,7 +2815,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -2834,7 +2834,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -2854,7 +2854,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -2879,7 +2879,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -2898,7 +2898,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -2918,7 +2918,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -2944,7 +2944,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2963,7 +2963,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -2983,7 +2983,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3013,7 +3013,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3036,7 +3036,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3060,7 +3060,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3089,7 +3089,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3112,7 +3112,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3136,7 +3136,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3166,7 +3166,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3189,7 +3189,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3213,7 +3213,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3254,7 +3254,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -3274,7 +3274,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -3295,7 +3295,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -3321,7 +3321,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -3341,7 +3341,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -3362,7 +3362,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -3389,7 +3389,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3409,7 +3409,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3430,7 +3430,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3459,7 +3459,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3483,7 +3483,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3508,7 +3508,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3538,7 +3538,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3562,7 +3562,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3587,7 +3587,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3618,7 +3618,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3642,7 +3642,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3667,7 +3667,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3707,7 +3707,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -3731,7 +3731,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -3755,7 +3755,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3782,7 +3782,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3810,7 +3810,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -3838,7 +3838,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3876,7 +3876,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); @@ -3900,7 +3900,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[(iy-1)*n1*n2+ix*n1+iz]) + c2*(tzz[(iy+1)*n1*n2+ix*n1+iz] - tzz[(iy-2)*n1*n2+ix*n1+iz])); @@ -3924,7 +3924,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); @@ -3951,7 +3951,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vx[iy*n1*n2+ix*n1+iz] -= rox[iy*n1*n2+ix*n1+iz]*( + vx[iy*n1*n2+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -3979,7 +3979,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+iz]*( + vy[iy*n1*n2+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -4007,7 +4007,7 @@ MID left mid mid #pragma ivdep for (iy=iyo; iy<iye; iy++) { for (iz=izo; iz<ize; iz++) { - vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+iz]*( + vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -4040,7 +4040,7 @@ MID left mid mid return 0; } -long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose) +long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long itime, long verbose) { /********************************************************************* @@ -4146,7 +4146,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[iy*n2*n1+ix*n1+iz+2] - vz[iy*n2*n1+ix*n1+iz-1]); Jz = RA[ipml2]*RA[ipml]*dvz - RA[ipml2]*RA[ipml]*dt*Pzpml[ix*npml+ipml]; Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz+dvx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz+dvx); ipml--; } } @@ -4166,7 +4166,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml]; Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx+dvz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx+dvz); ipml--; } } @@ -4183,7 +4183,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml]; Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx); ipml--; } } @@ -4196,7 +4196,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[ix*npml+ipml]; Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz); ipml--; } } @@ -4215,7 +4215,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml]; Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx+dvz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx+dvz); ipml++; } } @@ -4232,7 +4232,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml]; Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx); ipml++; } } @@ -4245,7 +4245,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[ix*npml+ipml]; Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz); ipml--; } } @@ -4264,7 +4264,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml]; Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz+dvx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz+dvx); ipml++; } } @@ -4281,7 +4281,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml]; Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz); ipml++; } } @@ -4294,7 +4294,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml]; Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx); //p[ix*n1+iz] -= l2m[ix*n1+iz]*(dvx); ipml++; } @@ -4312,7 +4312,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]); Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml]; Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jz); ipml++; } } @@ -4325,7 +4325,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]); Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml]; Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx; - p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx); + p[ix*n1+iz] -= l2m[iy][ix][iz]*(Jx); ipml--; } } diff --git a/fdelmodc3D/demo/model.scr b/fdelmodc3D/demo/model.scr index 41b1f5c21d5f8be0c68682981569198bb58dbddd..473fa61b865ee06e278de115f18a2a2c6aae86d4 100755 --- a/fdelmodc3D/demo/model.scr +++ b/fdelmodc3D/demo/model.scr @@ -1,9 +1,11 @@ #!/bin/bash #SBATCH -J OpenMP-test #SBATCH --nodes=1 -#SBATCH --ntasks-per-node=40 +#SBATCH --ntasks-per-node=8 #SBATCH --time=0:15:00 +cd /vardim/home/thorbcke/src/OpenSource/fdelmodc3D/demo + cp=2000 rho=1000 dx=2.5 @@ -17,6 +19,7 @@ makemod sizex=1000 sizez=500 dx=$dx dz=$dx cp0=$cp cs0=$cs ro0=$rho \ #extend to 3D ny=201 +nx=401 dy=2500 rm ro3d.su rm cp3d.su @@ -32,17 +35,29 @@ for (( i=0; i<$ny; i++ ));do sushw < syncl_ro.su key=gy,scalel a=$yloc,-1000 >> ro3d.su done + +suwind key=gx min=0 max=0 < syncl_cp.su > cp1d.su +suwind key=gx min=0 max=0 < syncl_ro.su > ro1d.su makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 # -export KMP_AFFINITY=verbose,granularity=fine,compact,1,0 +export KMP_AFFINITY=granularity=fine,scatter,0,0 export OMP_NUM_THREADS=8 +# file_cp=cp3d.su ischeme=1 \ +# file_den=ro3d.su \ + +# file_cp=syncl_cp.su ischeme=1 \ +# file_den=syncl_ro.su \ +# ny=$ny \ + time ../fdelmodc3D \ - file_cp=cp3d.su ischeme=1 \ - file_den=ro3d.su \ + file_cp=cp1d.su ischeme=1 \ + file_den=ro1d.su \ + ny=$ny \ + nx=$nx \ file_src=wave.su \ file_rcv=shot3d_fd.su \ src_type=7 \ @@ -60,5 +75,5 @@ time ../fdelmodc3D \ zrcv1=0 zrcv2=0 \ xsrc=0 zsrc=0 ysrc=0 \ ntaper=101 \ - left=2 right=2 bottom=2 top=1 + left=4 right=4 bottom=4 top=1 diff --git a/fdelmodc3D/elastic4dc_3D.c b/fdelmodc3D/elastic4dc_3D.c index 4a86c78ea9b55b97fd81c476ad4c1a2923f6b4f2..f2639ad0c7a9dd868adc4799467d62c1eb3f4a8a 100644 --- a/fdelmodc3D/elastic4dc_3D.c +++ b/fdelmodc3D/elastic4dc_3D.c @@ -8,7 +8,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose); long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose); @@ -17,16 +17,16 @@ long reStoreSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd, long ixsrc, lo float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, long verbose); long boundariesP3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, - float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *mul, long itime, long verbose); + float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***mul, long itime, long verbose); long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, - float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *mul, long itime, long verbose); + float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***mul, long itime, long verbose); long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long verbose) + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long verbose) { /********************************************************************* COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID: @@ -90,7 +90,7 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l for (ix=mod.ioXx; ix<mod.ieXx; ix++) { #pragma ivdep for (iz=mod.ioXz; iz<mod.ieXz; iz++) { - vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*( + vx[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*( c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] + txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) + @@ -107,7 +107,7 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l for (ix=mod.ioYx; ix<mod.ieYx; ix++) { #pragma ivdep for (iz=mod.ioYz; iz<mod.ieYz; iz++) { - vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*( + vy[iy*n2*n1+ix*n1+iz] -= roy[iy][ix][iz]*( c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] + tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] + txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) + @@ -124,7 +124,7 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l for (ix=mod.ioZx; ix<mod.ieZx; ix++) { #pragma ivdep for (iz=mod.ioZz; iz<mod.ieZz; iz++) { - vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*( + vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*( c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] + tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] + txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) + @@ -155,9 +155,9 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l c2*(vy[(iy+2)*n2*n1+ix*n1+iz] - vy[(iy-1)*n2*n1+ix*n1+iz]); dvz = c1*(vz[iy*n2*n1+ix*n1+iz+1] - vz[iy*n2*n1+ix*n1+iz]) + c2*(vz[iy*n2*n1+ix*n1+iz+2] - vz[iy*n2*n1+ix*n1+iz-1]); - txx[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*dvx + lam[iy*n2*n1+ix*n1+iz]*(dvy + dvz); - tyy[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*dvy + lam[iy*n2*n1+ix*n1+iz]*(dvz + dvx); - tzz[iy*n2*n1+ix*n1+iz] -= l2m[iy*n2*n1+ix*n1+iz]*dvz + lam[iy*n2*n1+ix*n1+iz]*(dvx + dvy); + txx[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*dvx + lam[iy][ix][iz]*(dvy + dvz); + tyy[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*dvy + lam[iy][ix][iz]*(dvz + dvx); + tzz[iy*n2*n1+ix*n1+iz] -= l2m[iy][ix][iz]*dvz + lam[iy][ix][iz]*(dvx + dvy); } } } @@ -177,4 +177,4 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l reStoreSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, verbose); return 0; -} \ No newline at end of file +} diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c index 21b3f5858add01f7eca9c2d5d8597fbbd27926c0..c9b31306fde496f3f4a108de882f6016caeaa386 100644 --- a/fdelmodc3D/fdelmodc3D.c +++ b/fdelmodc3D/fdelmodc3D.c @@ -17,11 +17,14 @@ double wallclock_time(void); void threadAffinity(void); +float ***alloc3float(modPar mod); +void free3float(float ***p); + long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, long verbose); -long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *muu, float *tss, float *tes, float *tep); +long readModel3D(modPar mod, bndPar bnd, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***muu, float *tss, float *tes, float *tep); long defineSource3D(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_nwav, long reverse, long verbose); @@ -30,32 +33,32 @@ long writeSrcRecPos3D(modPar *mod, recPar *rec, srcPar *src, shotPar *shot); long acoustic6_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose); + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose); long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose); + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose); long acoustic2_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose); + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose); long acoustic4_qr_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, - float *p, float *rox, float *roy, float *roz, float *l2m, long verbose); + float *p, float ***rox, float ***roy, float ***roz, float ***l2m, long verbose); long acousticSH4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *tx, float *ty, float *tz, - float *vz, float *rox, float *roy, float *roz, float *mul, long verbose); + float *vz, float ***rox, float ***roy, float ***roz, float ***mul, long verbose); long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, - float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long verbose); + float ***rox, float ***roy, float ***roz, float ***l2m, float ***lam, float ***mul, long verbose); long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, - float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz, + float *txz, float *txy, float *tyz, float ***l2m, float ***rox, float ***roy, float ***roz, float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_txy, float *rec_tyz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, long verbose); @@ -285,8 +288,7 @@ int main(int argc, char **argv) bndPar bnd; shotPar shot; float **src_nwav; - float *rox, *roy, *roz, *l2m, *lam, *mul; - float ***rox3, ***roy3, ***roz3, ***l2m3, ***lam3, ***mul3; + float ***rox, ***roy, ***roz, ***l2m, ***lam, ***mul; float *tss, *tes, *tep, *p, *q, *r; float *vx, *vy, *vz, *tzz, *tyy, *txz, *txy, *tyz, *txx; float *rec_vx, *rec_vy, *rec_vz, *rec_p; @@ -330,28 +332,18 @@ int main(int argc, char **argv) n2 = mod.nax; sizem=mod.nax*mod.naz*mod.nay; - rox3 = (float ***)alloc3float(mod); - roy3 = (float ***)alloc3float(mod); - roz3 = (float ***)alloc3float(mod); - l2m3 = (float ***)alloc3float(mod); - rox = (float *)rox3[0][0]; - roy = (float *)roy3[0][0]; - roz = (float *)roz3[0][0]; - l2m = (float *)l2m3[0][0]; -/* - rox = (float *)calloc(sizem,sizeof(float)); - roy = (float *)calloc(sizem,sizeof(float)); - roz = (float *)calloc(sizem,sizeof(float)); - l2m = (float *)calloc(sizem,sizeof(float)); -*/ + rox = (float ***)alloc3float(mod); + roy = (float ***)alloc3float(mod); + roz = (float ***)alloc3float(mod); + l2m = (float ***)alloc3float(mod); if (mod.ischeme==2) { tss = (float *)calloc(sizem,sizeof(float)); tep = (float *)calloc(sizem,sizeof(float)); q = (float *)calloc(sizem,sizeof(float)); } if (mod.ischeme>2) { - lam = (float *)calloc(sizem,sizeof(float)); - mul = (float *)calloc(sizem,sizeof(float)); + lam = (float ***)alloc3float(mod); + mul = (float ***)alloc3float(mod); } if (mod.ischeme==4) { tss = (float *)calloc(sizem,sizeof(float)); @@ -395,7 +387,6 @@ int main(int argc, char **argv) defineSource3D(wav, src, mod, rec, src_nwav, mod.grid_dir, verbose); - /* allocate arrays for wavefield and receiver arrays */ vx = (float *)calloc(sizem,sizeof(float)); @@ -428,9 +419,9 @@ int main(int argc, char **argv) rec_udp = (float *)calloc(mod.nax*rec.nt,sizeof(float)); } /* get velocity and density at first receiver location */ - ir = mod.ioZz + rec.z[0]+(rec.x[0]+mod.ioZx)*n1+(rec.y[0]+mod.ioZy)*n1*n2; - rec.rho = mod.dt/(mod.dx*roz[ir]); - rec.cp = sqrt(l2m[ir]*(roz[ir]))*mod.dx/mod.dt; + //ir = mod.ioZz + rec.z[0]+(rec.x[0]+mod.ioZx)*n1+(rec.y[0]+mod.ioZy)*n1*n2; + rec.rho = mod.dt/(mod.dx*roz[rec.y[0]+mod.ioZy][rec.x[0]+mod.ioZx][mod.ioZz+rec.z[0]]); + rec.cp = sqrt(l2m[rec.y[0]+mod.ioZy][rec.x[0]+mod.ioZx][mod.ioZz+rec.z[0]]*(roz[rec.y[0]+mod.ioZy][rec.x[0]+mod.ioZx][mod.ioZz+rec.z[0]]))*mod.dx/mod.dt; if(sna.beam) { @@ -473,7 +464,7 @@ int main(int argc, char **argv) if (bnd.lef==4 || bnd.lef==2) ioPx += bnd.ntap; if (bnd.fro==4 || bnd.fro==2) ioPy += bnd.ntap; if (bnd.top==4 || bnd.top==2) ioPz += bnd.ntap; - if (rec.sinkvel) sinkvel=l2m[(rec.y[0]+ioPy)*n1*n2+(rec.x[0]+ioPx)*n1+rec.z[0]+ioPz]; + if (rec.sinkvel) sinkvel=l2m[rec.y[0]+ioPy][rec.x[0]+ioPx][rec.z[0]+ioPz]; else sinkvel = 0.0; /* sink receivers to value different than sinkvel */ @@ -481,7 +472,7 @@ int main(int argc, char **argv) iz = rec.z[ir]; iy = rec.y[ir]; ix = rec.x[ir]; - while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz+ioPz] == sinkvel) iz++; + while(l2m[iy+ioPy][ix+ioPx][iz+ioPz] == sinkvel) iz++; rec.z[ir]=iz+rec.sinkdepth; rec.zr[ir]=rec.zr[ir]+(rec.z[ir]-iz)*mod.dz; if (verbose>3) vmess("receiver position %li at grid[ix=%li, iy=%li iz=%li] = (x=%f y=%f z=%f)", ir, ix+ioPx, iy+ioPy, rec.z[ir]+ioPz, rec.xr[ir]+mod.x0, rec.yr[ir]+mod.y0, rec.zr[ir]+mod.z0); @@ -493,7 +484,7 @@ int main(int argc, char **argv) iz = shot.z[ishot]; iy = shot.y[ishot]; ix = shot.x[ishot]; - while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz+ioPz] == 0.0) iz++; + while(l2m[iy+ioPy][ix+ioPx][iz+ioPz] == 0.0) iz++; shot.z[ishot]=iz+src.sinkdepth; } @@ -501,7 +492,7 @@ int main(int argc, char **argv) for (iy=0; iy<mod.ny; iy++) { for (ix=0; ix<mod.nx; ix++) { iz = ioPz; - while(l2m[(iy+ioPy)*n1*n2+(ix+ioPx)*n1+iz] == 0.0) iz++; + while(l2m[iy+ioPy][ix+ioPx][iz] == 0.0) iz++; bnd.surface[(iy+ioPy)*n2+ix+ioPx] = iz; if ((verbose>3) && (iz != ioPz)) vmess("Topography surface x=%.2f y=%.2f z=%.2f", mod.x0+mod.dx*ix, mod.y0+mod.dy*iy, mod.z0+mod.dz*(iz-ioPz)); } @@ -594,6 +585,8 @@ int main(int argc, char **argv) /* Main loop over the number of time steps */ for (it=it0; it<it1; it++) { + vmess("running time step time = %d.",it); + #pragma omp parallel default (shared) \ shared (rox, roy, roz, l2m, lam, mul, txx, tyy, tzz, txz, tyz, txy, vx, vy, vz) \ shared (tss, tep, tes, r, q, p) \ @@ -684,7 +677,6 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos if (sna.nsnap) { writeSnapTimes3D(mod, sna, bnd, wav, ixsrc, iysrc, izsrc, it, vx, vy, vz, tzz, tyy, txx, txz, tyz, txy, verbose); - } /* calculate beams */ @@ -751,10 +743,10 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos /* free arrays */ initargs(argc,argv); /* this will free the arg arrays declared */ - free(rox); - free(roy); - free(roz); - free(l2m); + free3float(rox); + free3float(roy); + free3float(roz); + free3float(l2m); free(src_nwav[0]); free(src_nwav); free(vx); @@ -799,8 +791,8 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos free(q); } if (mod.ischeme>2) { - free(lam); - free(mul); + free3float(lam); + free3float(mul); free(txz); free(tyz); free(txy); diff --git a/fdelmodc3D/getModelInfo3D.c b/fdelmodc3D/getModelInfo3D.c index 9705935cb990e1e1b6a00d68efa3883ed06d621d..211000d5213c980c7fed43eb40dd7f2a434739bf 100644 --- a/fdelmodc3D/getModelInfo3D.c +++ b/fdelmodc3D/getModelInfo3D.c @@ -53,9 +53,6 @@ long getModelInfo3D(char *file_name, long *n1, long *n2, long *n3, ny = 1; gy = hdr.gy; - if ( NINT(100.0*((*d1)/(*d2)))!=100 ) { - verr("dx and dz are different in the model !"); - } if ( NINT(1000.0*(*d1))==0 ) { if(!getparfloat("dx",d1)) { verr("dx is equal to zero use parameter dx= to set value"); @@ -66,6 +63,12 @@ long getModelInfo3D(char *file_name, long *n1, long *n2, long *n3, ntraces = (long) (bytes/trace_sz); *n2 = ntraces; + if (*n2==1) *d2 = *d1; + + if ( NINT(100.0*((*d1)/(*d2)))!=100 ) { + verr("dx and dz are different in the model !"); + } + /* check to find out min and max values gather */ one_shot = 1; @@ -112,8 +115,12 @@ long getModelInfo3D(char *file_name, long *n1, long *n2, long *n3, *n3 = ny; *n2 = ntraces/ny; - //*d3 = ((float)(gy-gy0))/((float)ny); //correcting - *d3 = ((float)(gy-gy0))/(((float)(ny-1))*1000); //del + if (ny==1) { + *d3 = *d2; + } + else { + *d3 = ((float)(gy-gy0))/(((float)(ny-1))*1000); //del + } if ( NINT(100.0*((*d1)/(*d3)))!=100 ) { verr("dx and dy are different in the model !"); diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c index 836ce4cf0ed7e9806a9fde628031286f2caf64d6..25b7f221f2f82b9752414f3e33f3d15b5c7325c4 100644 --- a/fdelmodc3D/getParameters3D.c +++ b/fdelmodc3D/getParameters3D.c @@ -148,10 +148,14 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar if (!getparlong("ny",&ny)) ny=nx; dx=dz; dy=dx; + sub_x0=-0.5*(nx-1)*(dx); + sub_y0=-0.5*(ny-1)*(dy); } else if (ny==1) { // 2D model if (!getparlong("ny",&ny)) ny=nx; dy=dx; + sub_y0=-0.5*(ny-1)*(dy); + fprintf(stderr,"get param ny=%d dy=%f y0=%f\n", ny, dy, sub_y0); } mod->dz = dz; mod->dx = dx; @@ -271,6 +275,15 @@ criteria we have imposed.*/ vmess("*******************************************"); vmess("*************** model info ****************"); vmess("*******************************************"); + if (mod.nfx == 1) { + vmess(" 1-dimensional model is read from file"); + } + else if (mod.nfy == 1) { + vmess(" 2-dimensional model is read from file"); + } + else { + vmess(" 3-dimensional model is read from file"); + } vmess("nz = %li ny = %li nx = %li", nz, ny, nx); vmess("dz = %8.4f dy = %8.4f dx = %8.4f", dz, dy, dx); vmess("zmin = %8.4f zmax = %8.4f", sub_z0, zmax); diff --git a/fdelmodc3D/getRecTimes3D.c b/fdelmodc3D/getRecTimes3D.c index 4c19786a8028eb10c0bc3183fe3c31c72e3603e8..16ffef78040ad2bf0226820517b795758d902410 100644 --- a/fdelmodc3D/getRecTimes3D.c +++ b/fdelmodc3D/getRecTimes3D.c @@ -19,7 +19,7 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, - float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz, + float *txz, float *txy, float *tyz, float ***l2m, float ***rox, float ***roy, float ***roz, float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_txy, float *rec_tyz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, long verbose) @@ -361,14 +361,14 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]); dvz = c1*(vz[iy*n1*n2+ix*n1+iz+1] - vz[iy*n1*n2+ix*n1+iz]) + c2*(vz[iy*n1*n2+ix*n1+iz+2] - vz[iy*n1*n2+ix*n1+iz-1]); - field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz]*(dvx+dvy+dvz); + field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy][ix][iz]*(dvx+dvy+dvz); dvx = c1*(vx[iy*n1*n2+(ix+1)*n1+iz1] - vx[iy*n1*n2+ix*n1+iz1]) + c2*(vx[iy*n1*n2+(ix+2)*n1+iz1] - vx[iy*n1*n2+(ix-1)*n1+iz1]); dvy = c1*(vy[(iy+1)*n1*n2+ix*n1+iz1] - vy[iy*n1*n2+ix*n1+iz1]) + c2*(vy[(iy+2)*n1*n2+ix*n1+iz1] - vy[(iy-1)*n1*n2+ix*n1+iz1]); dvz = c1*(vz[iy*n1*n2+ix*n1+iz1+1] - vz[iy*n1*n2+ix*n1+iz1]) + c2*(vz[iy*n1*n2+ix*n1+iz1+2] - vz[iy*n1*n2+ix*n1+iz1-1]); - field += tzz[iy*n1*n2+ix*n1+iz1] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz1]*(dvx+dvy+dvz); + field += tzz[iy*n1*n2+ix*n1+iz1] + (1.0/2.0)*l2m[iy][ix][iz1]*(dvx+dvy+dvz); rec_p[irec*rec.nt+isam] = 0.5*field; } else { @@ -383,14 +383,14 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]); dvz = c1*(vz[iy*n1*n2+ix*n1+iz+1] - vz[iy*n1*n2+ix*n1+iz]) + c2*(vz[iy*n1*n2+ix*n1+iz+2] - vz[iy*n1*n2+ix*n1+iz-1]); - field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix*n1+iz]*(dvx+dvy+dvz); + field = tzz[iy*n1*n2+ix*n1+iz] + (1.0/2.0)*l2m[iy][ix][iz]*(dvx+dvy+dvz); dvx = c1*(vx[iy*n1*n2+(ix1+1)*n1+iz] - vx[iy*n1*n2+ix1*n1+iz]) + c2*(vx[iy*n1*n2+(ix1+2)*n1+iz] - vx[iy*n1*n2+(ix1-1)*n1+iz]); dvy = c1*(vy[(iy+1)*n1*n2+ix1*n1+iz] - vy[iy*n1*n2+ix1*n1+iz]) + c2*(vy[(iy+2)*n1*n2+ix1*n1+iz] - vy[(iy-1)*n1*n2+ix1*n1+iz]); dvz = c1*(vz[iy*n1*n2+ix1*n1+iz+1] - vz[iy*n1*n2+ix1*n1+iz]) + c2*(vz[iy*n1*n2+ix1*n1+iz+2] - vz[iy*n1*n2+ix1*n1+iz-1]); - field += tzz[iy*n1*n2+ix1*n1+iz] + (1.0/2.0)*l2m[iy*n1*n2+ix1*n1+iz]*(dvx+dvy+dvz); + field += tzz[iy*n1*n2+ix1*n1+iz] + (1.0/2.0)*l2m[iy][ix][iz]*(dvx+dvy+dvz); rec_p[irec*rec.nt+isam] = 0.5*field; } else { @@ -438,10 +438,10 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, /* interpolate vz to Txx/Tzz position by taking the mean of 2 values */ else if (rec.int_vz == 2) { if (mod.ischeme == 1) { /* interpolate Vz times +1/2 Dt forward to P times */ - field = vz[iy*n1*n2+ix*n1+iz] - 0.5*roz[iy*n1*n2+ix*n1+iz]*( + field = vz[iy*n1*n2+ix*n1+iz] - 0.5*roz[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+ix*n1+iz-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2])); - field += vz[iy*n1*n2+ix*n1+iz2] - 0.5*roz[iy*n1*n2+ix*n1+iz2]*( + field += vz[iy*n1*n2+ix*n1+iz2] - 0.5*roz[iy][ix][iz2]*( c1*(tzz[iy*n1*n2+ix*n1+iz2] - tzz[iy*n1*n2+ix*n1+iz2-1]) + c2*(tzz[iy*n1*n2+ix*n1+iz2+1] - tzz[iy*n1*n2+ix*n1+iz2-2])); rec_vz[irec*rec.nt+isam] = 0.5*field; @@ -468,10 +468,10 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam, /* interpolate vx to Txx/Tzz position by taking the mean of 2 values */ else if (rec.int_vx == 2) { if (mod.ischeme == 1) { /* interpolate Vx times +1/2 Dt forward to P times */ - field = vx[iy*n1*n2+ix*n1+iz] - 0.5*rox[iy*n1*n2+ix*n1+iz]*( + field = vx[iy*n1*n2+ix*n1+iz] - 0.5*rox[iy][ix][iz]*( c1*(tzz[iy*n1*n2+ix*n1+iz] - tzz[iy*n1*n2+(ix-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix+1)*n1+iz] - tzz[iy*n1*n2+(ix-2)*n1+iz])); - field += vx[ix2*n1+iz] - 0.5*rox[ix2*n1+iz]*( + field += vx[ix2*n1+iz] - 0.5*rox[iy][ix2][iz]*( c1*(tzz[iy*n1*n2+ix2*n1+iz] - tzz[iy*n1*n2+(ix2-1)*n1+iz]) + c2*(tzz[iy*n1*n2+(ix2+1)*n1+iz] - tzz[iy*n1*n2+(ix2-2)*n1+iz])); rec_vx[irec*rec.nt+isam] = 0.5*field; diff --git a/fdelmodc3D/readModel3D.c b/fdelmodc3D/readModel3D.c index 157f7a1790d4bcb1c8539aebf37385b170679e63..61dcea102609f5b9d568f8b86d96cea5daba37c2 100644 --- a/fdelmodc3D/readModel3D.c +++ b/fdelmodc3D/readModel3D.c @@ -15,6 +15,9 @@ #define MIN(x,y) ((x) < (y) ? (x) : (y)) #define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5)) +float ***alloc3float(modPar mod); +void free3float(float ***p); + /** * Reads gridded model files and compute from them medium parameters used in the FD kernels. * The files read in contain the P (and S) wave velocity and density. @@ -26,33 +29,35 @@ **/ -long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, - float *l2m, float *lam, float *muxz, float *tss, float *tes, float *tep) +long readModel3D(modPar mod, bndPar bnd, float ***rox, float ***roy, float ***roz, + float ***l2m, float ***lam, float ***muxz, float *tss, float *tes, float *tep) { FILE *fpcp, *fpcs, *fpro; FILE *fpqp=NULL, *fpqs=NULL; size_t nread; long i, j, l, itmp, tracesToDo; - long n1, n2, n3, ix, iy, iz, nz, ny, nx; + long n1, n2, n3, ix, iy, iz, nz, ny, nx, nfz, nfx, nfy; long ixo, iyo, izo, ixe, iye, ize; long ioXx, ioXy, ioXz, ioYx, ioYy, ioYz, ioZz, ioZy, ioZx, ioPx, ioPy, ioPz, ioTx, ioTy, ioTz; float cp2, cs2, cs111, cs112, cs121, cs211, cs122, cs212, cs221, mul, mu, lamda2mu, lamda; float cs2c, cs2b, cs2a, cpx, cpy, cpz, bx, by, bz, fac; - float *cp, *cs, *ro, *qp, *qs; + float ***cp, ***cs, ***ro, *qp, *qs; float a, b; segy hdr; /* grid size and start positions for the components */ nz = mod.nz; - ny = mod.ny; nx = mod.nx; + ny = mod.ny; n1 = mod.naz; n2 = mod.nax; n3 = mod.nay; + nfz = mod.nfz; + nfx = mod.nfx; + nfy = mod.nfy; fac = mod.dt/mod.dx; - /* Vx: rox */ ioXx=mod.ioXx; ioXy=mod.ioXy; @@ -88,27 +93,29 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* open files and read first header */ - cp = (float *)malloc(nz*ny*nx*sizeof(float)); + //cp = (float *)malloc(nz*ny*nx*sizeof(float)); + cp = (float ***)alloc3float(mod); fpcp = fopen( mod.file_cp, "r" ); assert( fpcp != NULL); nread = fread(&hdr, 1, TRCBYTES, fpcp); assert(nread == TRCBYTES); - ro = (float *)malloc(nz*ny*nx*sizeof(float)); + //ro = (float *)malloc(nz*ny*nx*sizeof(float)); + ro = (float ***)alloc3float(mod); fpro = fopen( mod.file_ro, "r" ); assert( fpro != NULL); nread = fread(&hdr, 1, TRCBYTES, fpro); assert(nread == TRCBYTES); - cs = (float *)calloc(nz*ny*nx,sizeof(float)); + //cs = (float *)calloc(nz*ny*nx,sizeof(float)); if (mod.ischeme>2 && mod.ischeme!=5) { + cs = (float ***)alloc3float(mod); fpcs = fopen( mod.file_cs, "r" ); assert( fpcs != NULL); nread = fread(&hdr, 1, TRCBYTES, fpcs); assert(nread == TRCBYTES); } - /* for visco acoustic/elastic media open Q file(s) if given as parameter */ if (mod.file_qp != NULL && (mod.ischeme==2 || mod.ischeme==4)) { @@ -130,15 +137,16 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* read all traces */ tracesToDo = mod.nx*mod.ny; - for (iy=0; iy<ny; iy++) { - for (ix=0; ix<nx; ix++ ) { - i = iy*nx+ix; - nread = fread(&cp[i*nz], sizeof(float), hdr.ns, fpcp); + for (iy=0; iy<nfy; iy++) { + for (ix=0; ix<nfx; ix++ ) { + //i = iy*nx+ix; +// fprintf(stderr,"iy=%d ix=%d\n", iy, ix); + nread = fread(&cp[iy][ix][0], sizeof(float), hdr.ns, fpcp); assert (nread == hdr.ns); - nread = fread(&ro[i*nz], sizeof(float), hdr.ns, fpro); + nread = fread(&ro[iy][ix][0], sizeof(float), hdr.ns, fpro); assert (nread == hdr.ns); if (mod.ischeme>2 && mod.ischeme!=5) { - nread = fread(&cs[i*nz], sizeof(float), hdr.ns, fpcs); + nread = fread(&cs[iy][ix][0], sizeof(float), hdr.ns, fpcs); assert (nread == hdr.ns); } @@ -218,10 +226,14 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* check for zero densities */ - for (i=0;i<nz*nx*ny;i++) { - if (ro[i]==0.0) { - vwarn("Zero density for trace=%li sample=%li", i/nz, i%nz); - verr("ERROR zero density is not a valid value, program exit"); + for (iy=0;iy<ny;iy++) { + for (ix=0;ix<nx;ix++) { + for (iz=0;iz<nz;iz++) { + if (ro[iy][ix][iz]==0.0) { + vwarn("Zero density for trace=[%li][%li][%li]", iy, ix, iz); + verr("ERROR zero density is not a valid value, program exit"); + } + } } } @@ -236,13 +248,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (ix=0;ix<nx-1;ix++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - cs212 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2*ro[iy][ix][iz]; + cs211 = cs2a*ro[iy][ix+1][iz]; + cs212 = cs2a*ro[iy][ix+1][iz]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -252,13 +264,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs112 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs122 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs112 = cs2a*ro[iy+1][ix][iz]; + // cs122 = cs2a*ro[iy+1][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -268,15 +280,15 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs2b = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - // cs2c = cs[(iy+1)*nx*nz+(ix+1)*nz+iz]*cs[(iy+1)*nx*nz+(ix+1)*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs211 = cs2b*ro[iy*nx*nz+(ix+1)*nz+iz]; - // cs221 = cs2c*ro[(iy+1)*nx*nz+(ix+1)*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs2b = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + // cs2c = cs[iy+1][ix+1][iz]*cs[iy+1][ix+1][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2a*ro[iy+1][ix][iz]; + // cs211 = cs2b*ro[iy][ix+1][iz]; + // cs221 = cs2c*ro[iy+1][ix+1][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -284,19 +296,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } } @@ -305,15 +317,15 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iz=0;iz<nz-1;iz++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - cs2c = cs[iy*nx*nz+(ix+1)*nz+iz+1]*cs[iy*nx*nz+(ix+1)*nz+iz+1]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2b*ro[iy*nx*nz+ix*nz+iz+1]; - cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - cs212 = cs2c*ro[iy*nx*nz+(ix+1)*nz+iz+1]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + cs2b = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + cs2c = cs[iy][ix+1][iz+1]*cs[iy][ix+1][iz+1]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2b*ro[iy][ix][iz+1]; + cs211 = cs2a*ro[iy][ix+1][iz]; + cs212 = cs2c*ro[iy][ix+1][iz+1]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -323,13 +335,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - // cs122 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs112 = cs2a*ro[iy][ix][iz+1]; + // cs122 = cs2a*ro[iy][ix][iz+1]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -339,13 +351,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - // cs221 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs211 = cs2a*ro[iy][ix+1][iz]; + // cs221 = cs2a*ro[iy][ix+1][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -353,19 +365,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = ro[iy*nx*nz+ix*nz+iz]; - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = ro[iy][ix][iz]; + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } } @@ -374,13 +386,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iz=0;iz<nz-1;iz++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs212 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2a*ro[iy][ix][iz+1]; + cs211 = cs2*ro[iy][ix][iz]; + cs212 = cs2a*ro[iy][ix][iz+1]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -390,15 +402,15 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs2c = cs[(iy+1)*nx*nz+ix*nz+iz+1]*cs[(iy+1)*nx*nz+ix*nz+iz+1]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - // cs122 = cs2c*ro[(iy+1)*nx*nz+ix*nz+iz+1]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + // cs2b = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs2c = cs[iy+1][ix][iz+1]*cs[iy+1][ix][iz+1]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2b*ro[iy+1][ix][iz]; + // cs112 = cs2a*ro[iy][ix][iz+1]; + // cs122 = cs2c*ro[iy+1][ix][iz+1]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -408,13 +420,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs221 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2a*ro[iy+1][ix][iz]; + // cs211 = cs2*ro[iy][ix][iz]; + // cs221 = cs2a*ro[iy+1][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -422,19 +434,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = ro[iy*nx*nz+ix*nz+iz]; - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = ro[iy][ix][iz]; + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } } @@ -443,13 +455,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (ix=0;ix<nx-1;ix++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - cs212 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2*ro[iy][ix][iz]; + cs211 = cs2a*ro[iy][ix+1][iz]; + cs212 = cs2a*ro[iy][ix+1][iz]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -459,12 +471,12 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs112 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs122 = cs2*ro[iy*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iiy][ix][iz*nz+iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs112 = cs2*ro[iy][ix][iz]; + // cs122 = cs2*ro[iy][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -474,13 +486,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - // cs221 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs211 = cs2a*ro[iy][ix+1][iz]; + // cs221 = cs2a*ro[iy][ix+1][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -488,19 +500,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = ro[iy*nx*nz+ix*nz+iz]; - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = ro[iy][ix][iz]; + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } ix = nx-1; @@ -508,12 +520,12 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=0;iy<ny-1;iy++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs212 = cs2*ro[iy*nx*nz+ix*nz+iz]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2*ro[iy][ix][iz]; + cs211 = cs2*ro[iy][ix][iz]; + cs212 = cs2*ro[iy][ix][iz]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -523,13 +535,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs112 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs122 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2a*ro[iy+1][ix][iz]; + // cs112 = cs2*ro[iy][ix][iz]; + // cs122 = cs2a*ro[iy+1][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -539,13 +551,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs221 = cs2a*ro[(iy+1)*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2a*ro[iy+1][ix][iz]; + // cs211 = cs2*ro[iy][ix][iz]; + // cs221 = cs2a*ro[iy+1][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -553,19 +565,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = ro[iy*nx*nz+ix*nz+iz]; - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/bx; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = ro[iy][ix][iz]; + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/bx; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } ix = nx-1; @@ -573,13 +585,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iz=0;iz<nz-1;iz++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs212 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2a*ro[iy][ix][iz+1]; + cs211 = cs2*ro[iy][ix][iz]; + cs212 = cs2a*ro[iy][ix][iz+1]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -589,13 +601,13 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - // cs122 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs112 = cs2a*ro[iy][ix][iz+1]; + // cs122 = cs2a*ro[iy][ix][iz+1]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -605,12 +617,12 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs211 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs221 = cs2*ro[iy*nx*nz+ix*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2*ro[iy][ix][iz]; + // cs211 = cs2*ro[iy][ix][iz]; + // cs221 = cs2*ro[iy][ix][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -618,53 +630,53 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = ro[iy*nx*nz+ix*nz+iz]; - by = ro[iy*nx*nz+ix*nz+iz]; - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/bx; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = ro[iy][ix][iz]; + by = ro[iy][ix][iz]; + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/bx; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } ix=nx-1; iy=ny-1; iz=nz-1; - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = ro[iy*nx*nz+ix*nz+iz]; - by = ro[iy*nx*nz+ix*nz+iz]; - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mu; + bx = ro[iy][ix][iz]; + by = ro[iy][ix][iz]; + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mu; for (iy=0;iy<ny-1;iy++) { for (ix=0;ix<nx-1;ix++) { for (iz=0;iz<nz-1;iz++) { /* for muxz field */ /* csxyz */ - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - cs2c = cs[iy*nx*nz+(ix+1)*nz+iz+1]*cs[iy*nx*nz+(ix+1)*nz+iz+1]; - cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - cs112 = cs2b*ro[iy*nx*nz+ix*nz+iz+1]; - cs211 = cs2a*ro[iy*nx*nz+ix*nz+iz]; - cs212 = cs2c*ro[iy*nx*nz+ix*nz+iz+1]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + cs2b = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + cs2c = cs[iy][ix+1][iz+1]*cs[iy][ix+1][iz+1]; + cs111 = cs2*ro[iy][ix][iz]; + cs112 = cs2b*ro[iy][ix][iz+1]; + cs211 = cs2a*ro[iy][ix][iz]; + cs212 = cs2c*ro[iy][ix][iz+1]; if (cs111 > 0.0) { mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212); } @@ -674,15 +686,15 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muyz field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1]; - // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs2c = cs[(iy+1)*nx*nz+ix*nz+iz+1]*cs[(iy+1)*nx*nz+ix*nz+iz+1]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs112 = cs2a*ro[iy*nx*nz+ix*nz+iz+1]; - // cs122 = cs2c*ro[(iy+1)*nx*nz+ix*nz+iz+1]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix][iz+1]*cs[iy][ix][iz+1]; + // cs2b = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs2c = cs[iy+1][ix][iz+1]*cs[iy+1][ix][iz+1]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2b*ro[iy+1][ix][iz]; + // cs112 = cs2a*ro[iy][ix][iz+1]; + // cs122 = cs2c*ro[iy+1][ix][iz+1]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122); // } @@ -692,15 +704,15 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, /* for muxy field IN PROGRESS!!!!!!!!!!!!!!!!! */ /* csxyz */ - // cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - // cs2 = cs[iy*nx*nz+ix*nz+iz]*cs[iy*nx*nz+ix*nz+iz]; - // cs2a = cs[iy*nx*nz+(ix+1)*nz+iz]*cs[iy*nx*nz+(ix+1)*nz+iz]; - // cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz]; - // cs2c = cs[(iy+1)*nx*nz+(ix+1)*nz+iz]*cs[(iy+1)*nx*nz+(ix+1)*nz+iz]; - // cs111 = cs2*ro[iy*nx*nz+ix*nz+iz]; - // cs121 = cs2b*ro[(iy+1)*nx*nz+ix*nz+iz]; - // cs211 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz]; - // cs221 = cs2c*ro[(iy+1)*nx*nz+(ix+1)*nz+iz]; + // cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + // cs2 = cs[iy][ix][iz]*cs[iy][ix][iz]; + // cs2a = cs[iy][ix+1][iz]*cs[iy][ix+1][iz]; + // cs2b = cs[iy+1][ix][iz]*cs[iy+1][ix][iz]; + // cs2c = cs[iy+1][ix+1][iz]*cs[iy+1][ix+1][iz]; + // cs111 = cs2*ro[iy][ix][iz]; + // cs121 = cs2b*ro[iy+1][ix][iz]; + // cs211 = cs2a*ro[iy][ix+1][iz]; + // cs221 = cs2c*ro[iy+1][ix+1][iz]; // if (cs111 > 0.0) { // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221); // } @@ -708,19 +720,19 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, // mul = 0.0; // } - mu = cs2*ro[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + mu = cs2*ro[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; lamda = lamda2mu - 2*mu; - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; - lam[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda; - muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul; + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; + lam[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda; + muxz[iy+ioTy][ix+ioTx][iz+ioTz]=fac*mul; } } } @@ -730,123 +742,123 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, iz = nz-1; for (iy=0;iy<ny-1;iy++) { for (ix=0;ix<nx-1;ix++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } } iy = ny-1; for (iz=0;iz<nz-1;iz++) { for (ix=0;ix<nx-1;ix++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = ro[iy*nx*nz+ix*nz+iz]; - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = ro[iy][ix][iz]; + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } } ix = nx-1; for (iz=0;iz<nz-1;iz++) { for (iy=0;iy<ny-1;iy++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = ro[iy*nx*nz+ix*nz+iz]; - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = ro[iy][ix][iz]; + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } } iz = nz-1; iy = ny-1; for (ix=0;ix<nx-1;ix++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = ro[iy*nx*nz+ix*nz+iz]; - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = ro[iy][ix][iz]; + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } iz = nz-1; ix = nx-1; for (iy=0;iy<ny-1;iy++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = ro[iy*nx*nz+ix*nz+iz]; - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = ro[iy][ix][iz]; + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } ix = nx-1; iy = ny-1; for (iz=0;iz<nz-1;iz++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - - bx = ro[iy*nx*nz+ix*nz+iz]; - by = ro[iy*nx*nz+ix*nz+iz]; - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + + bx = ro[iy][ix][iz]; + by = ro[iy][ix][iz]; + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } ix=nx-1; iy=ny-1; iz=nz-1; - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; - bx = ro[iy*nx*nz+ix*nz+iz]; - by = ro[iy*nx*nz+ix*nz+iz]; - bz = ro[iy*nx*nz+ix*nz+iz]; - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; + bx = ro[iy][ix][iz]; + by = ro[iy][ix][iz]; + bz = ro[iy][ix][iz]; + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; for (iy=0; iy<ny-1; iy++) { for (ix=0; ix<nx-1; ix++) { for (iz=0; iz<nz-1; iz++) { - cp2 = cp[iy*nx*nz+ix*nz+iz]*cp[iy*nx*nz+ix*nz+iz]; - lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz]; + cp2 = cp[iy][ix][iz]*cp[iy][ix][iz]; + lamda2mu = cp2*ro[iy][ix][iz]; - bx = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+(ix+1)*nz+iz]); - by = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[(iy+1)*nx*nz+ix*nz+iz]); - bz = 0.5*(ro[iy*nx*nz+ix*nz+iz]+ro[iy*nx*nz+ix*nz+iz+1]); - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=fac/bx; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=fac/by; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=fac/bz; - l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu; + bx = 0.5*(ro[iy][ix][iz]+ro[iy][ix+1][iz]); + by = 0.5*(ro[iy][ix][iz]+ro[iy+1][ix][iz]); + bz = 0.5*(ro[iy][ix][iz]+ro[iy][ix][iz+1]); + rox[iy+ioXy][ix+ioXx][iz+ioXz]=fac/bx; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=fac/by; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=fac/bz; + l2m[iy+ioPy][ix+ioPx][iz+ioPz]=fac*lamda2mu; } } } @@ -857,10 +869,10 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=0; iy<ny; iy++) { for (ix=0; ix<nx; ix++) { for (iz=0; iz<nz; iz++) { - if (l2m[(iy+ioPy)*n2*n1+(ix+ioPx)*n1+iz+ioPz]==0.0) { - rox[(iy+ioXy)*n2*n1+(ix+ioXx)*n1+iz+ioXz]=0.0; - roy[(iy+ioYy)*n2*n1+(ix+ioYx)*n1+iz+ioYz]=0.0; - roz[(iy+ioZy)*n2*n1+(ix+ioZx)*n1+iz+ioZz]=0.0; + if (l2m[iy+ioPy][ix+ioPx][iz+ioPz]==0.0) { + rox[iy+ioXy][ix+ioXx][iz+ioXz]=0.0; + roy[iy+ioYy][ix+ioYx][iz+ioYz]=0.0; + roz[iy+ioZy][ix+ioZx][iz+ioZz]=0.0; } } } @@ -884,7 +896,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[iy*n2*n1+ixe*n1+iz]; + rox[iy][ix][iz] = rox[iy][ixe][iz]; } } } @@ -899,7 +911,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[iy*n2*n1+ixe*n1+iz]; + roy[iy][ix][iz] = roy[iy][ixe][iz]; } } } @@ -914,7 +926,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[iy*n2*n1+ixe*n1+iz]; + roz[iy][ix][iz] = roz[iy][ixe][iz]; } } } @@ -929,7 +941,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[iy*n2*n1+ixe*n1+iz]; + l2m[iy][ix][iz] = l2m[iy][ixe][iz]; } } } @@ -945,7 +957,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[iy*n2*n1+ixe*n1+iz]; + lam[iy][ix][iz] = lam[iy][ixe][iz]; } } } @@ -959,7 +971,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ixe*n1+iz]; + muxz[iy][ix][iz] = muxz[iy][ixe][iz]; } } } @@ -1014,7 +1026,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[iy*n2*n1+(ixo-1)*n1+iz]; + rox[iy][ix][iz] = rox[iy][ixo-1][iz]; } } } @@ -1029,7 +1041,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[iy*n2*n1+(ixo-1)*n1+iz]; + roy[iy][ix][iz] = roy[iy][ixo-1][iz]; } } } @@ -1044,7 +1056,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[iy*n2*n1+(ixo-1)*n1+iz]; + roz[iy][ix][iz] = roz[iy][ixo-1][iz]; } } } @@ -1059,7 +1071,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[iy*n2*n1+(ixo-1)*n1+iz]; + l2m[iy][ix][iz] = l2m[iy][ixo-1][iz]; } } } @@ -1075,7 +1087,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[iy*n2*n1+(ixo-1)*n1+iz]; + lam[iy][ix][iz] = lam[iy][ixo-1][iz]; } } } @@ -1090,7 +1102,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+(ixo-1)*n1+iz]; + muxz[iy][ix][iz] = muxz[iy][ixo-1][iz]; } } } @@ -1147,7 +1159,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[iye*n2*n1+ix*n1+iz]; + rox[iy][ix][iz] = rox[iye][ix][iz]; } } } @@ -1164,7 +1176,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[iye*n2*n1+ix*n1+iz]; + roy[iy][ix][iz] = roy[iye][ix][iz]; } } } @@ -1181,7 +1193,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[iye*n2*n1+ix*n1+iz]; + roz[iy][ix][iz] = roz[iye][ix][iz]; } } } @@ -1200,7 +1212,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[iye*n2*n1+ix*n1+iz]; + l2m[iy][ix][iz] = l2m[iye][ix][iz]; } } } @@ -1218,7 +1230,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[iye*n2*n1+ix*n1+iz]; + lam[iy][ix][iz] = lam[iye][ix][iz]; } } } @@ -1234,7 +1246,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[iye*n2*n1+ix*n1+iz]; + muxz[iy][ix][iz] = muxz[iye][ix][iz]; } } } @@ -1295,7 +1307,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[(iyo-1)*n2*n1+ix*n1+iz]; + rox[iy][ix][iz] = rox[iyo-1][ix][iz]; } } } @@ -1312,7 +1324,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[(iyo-1)*n2*n1+ix*n1+iz]; + roy[iy][ix][iz] = roy[iyo-1][ix][iz]; } } } @@ -1329,7 +1341,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[(iyo-1)*n2*n1+ix*n1+iz]; + roz[iy][ix][iz] = roz[iyo-1][ix][iz]; } } } @@ -1348,7 +1360,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[(iyo-1)*n2*n1+ix*n1+iz]; + l2m[iy][ix][iz] = l2m[iyo-1][ix][iz]; } } } @@ -1366,7 +1378,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[(iyo-1)*n2*n1+ix*n1+iz]; + lam[iy][ix][iz] = lam[iyo-1][ix][iz]; } } } @@ -1383,7 +1395,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[(iyo-1)*n2*n1+ix*n1+iz]; + muxz[iy][ix][iz] = muxz[iyo-1][ix][iz]; } } } @@ -1446,7 +1458,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[iy*n2*n1+ix*n1+ize]; + rox[iy][ix][iz] = rox[iy][ix][ize]; } } } @@ -1465,7 +1477,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[iy*n2*n1+ix*n1+ize]; + roy[iy][ix][iz] = roy[iy][ix][ize]; } } } @@ -1484,7 +1496,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[iy*n2*n1+ix*n1+ize]; + roz[iy][ix][iz] = roz[iy][ix][ize]; } } } @@ -1499,7 +1511,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[iy*n2*n1+ix*n1+ize]; + l2m[iy][ix][iz] = l2m[iy][ix][ize]; } } } @@ -1515,7 +1527,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[iy*n2*n1+ix*n1+ize]; + lam[iy][ix][iz] = lam[iy][ix][ize]; } } } @@ -1530,7 +1542,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+ize]; + muxz[iy][ix][iz] = muxz[iy][ix][ize]; } } } @@ -1589,7 +1601,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - rox[iy*n2*n1+ix*n1+iz] = rox[iy*n2*n1+ix*n1+izo-1]; + rox[iy][ix][iz] = rox[iy][ix][izo-1]; } } } @@ -1608,7 +1620,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roy[iy*n2*n1+ix*n1+iz] = roy[iy*n2*n1+ix*n1+izo-1]; + roy[iy][ix][iz] = roy[iy][ix][izo-1]; } } } @@ -1627,7 +1639,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - roz[iy*n2*n1+ix*n1+iz] = roz[iy*n2*n1+ix*n1+izo-1]; + roz[iy][ix][iz] = roz[iy][ix][izo-1]; } } } @@ -1641,7 +1653,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - l2m[iy*n2*n1+ix*n1+iz] = l2m[iy*n2*n1+ix*n1+izo-1]; + l2m[iy][ix][iz] = l2m[iy][ix][izo-1]; } } } @@ -1657,7 +1669,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - lam[iy*n2*n1+ix*n1+iz] = lam[iy*n2*n1+ix*n1+izo-1]; + lam[iy][ix][iz] = lam[iy][ix][izo-1]; } } } @@ -1672,7 +1684,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, for (iy=iyo; iy<iye; iy++) { for (ix=ixo; ix<ixe; ix++) { for (iz=izo; iz<ize; iz++) { - muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+izo-1]; + muxz[iy][ix][iz] = muxz[iy][ix][izo-1]; } } } @@ -1713,9 +1725,11 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz, } - free(cp); - free(ro); - free(cs); + free3float(cp); + free3float(ro); + if (mod.ischeme>2 && mod.ischeme!=5) { + free3float(cs); + } return 0; } diff --git a/marchenko/demo/invisible/primaries.scr b/marchenko/demo/invisible/primaries.scr index 26d4483628eadaaa521a6691576c7bcdb0f29344..d50529cd373cc4a6a270762068a1295fdf132792 100755 --- a/marchenko/demo/invisible/primaries.scr +++ b/marchenko/demo/invisible/primaries.scr @@ -16,7 +16,7 @@ makewave fp=25 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1 suwind key=fldr min=$select max=$select < shotsdx4_rp.su > shot${select}.su fconv file_in1=shot${select}.su file_in2=wave.su file_out=shotw.su -marchenko_primaries file_shot=shotsdx4_rp.su ishot=$select \ +../../marchenko_primaries file_shot=shotsdx4_rp.su ishot=$select \ file_src=wave.su file_rr=pred_inv.su \ verbose=2 istart=21 iend=751 fmax=90 \ niterskip=1024 niter=21 shift=20 diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c index 8c76552db29f7ee3066f8fb41cdb55a985522c81..677d6738bde64a1306d79e90c5fbe818dbaed60b 100644 --- a/marchenko/marchenko_primaries.c +++ b/marchenko/marchenko_primaries.c @@ -104,7 +104,7 @@ int main (int argc, char **argv) int size, n1, n2, ntap, tap, di, ntraces; int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn; int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft; - int iter, niter, niterec, niterskip, niterrun, tracf, *muteW; + int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW; int hw, ii, ishot, istart, iend; int smooth, *ixpos, npos, ix, m, pad, T; int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift; @@ -487,6 +487,7 @@ int main (int argc, char **argv) /* once every 'niterskip' time-steps start from fresh G_d and do niter (~20) iterations */ if ( ((ii-istart)%niterskip==0) || (ii==istart) ) { niterrun=niter; + recur=0; if (verbose) vmess("Doing %d iterations for time-sample %d\n",niterrun,ii); for (l = 0; l < Nfoc; l++) { for (i = 0; i < nxs; i++) { @@ -507,24 +508,23 @@ int main (int argc, char **argv) } } } - else { /* use f1min as stating point and do only 2 iterations */ + else { /* use f1min from previous iteration as starting point and do niterec iterations */ niterrun=niterec; + recur=1; if (verbose>1) vmess("Doing %d iterations for time-sample %d",niterrun,ii); for (l = 0; l < Nfoc; l++) { for (i = 0; i < npos; i++) { j=0; ix = ixpos[i]; - G_d[l*nxs*nts+i*nts+j] = +f1min[l*nxs*nts+i*nts+j]; + G_d[l*nxs*nts+i*nts+j] = -f1min[l*nxs*nts+i*nts+j]; for (j = 1; j < nts; j++) { - G_d[l*nxs*nts+i*nts+j] = +f1min[l*nxs*nts+i*nts+nts+j]; + G_d[l*nxs*nts+i*nts+j] = -f1min[l*nxs*nts+i*nts+nts-j]; } /* org G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j]; for (j = 1; j < nts; j++) { G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j]; } -*/ -/* for (j = 0; j < nts-ii+T*shift; j++) { G_d[l*nxs*nts+i*nts+j] = 0.0; } @@ -578,7 +578,7 @@ int main (int argc, char **argv) for (l = 0; l < Nfoc; l++) { for (i = 0; i < npos; i++) { ix = ixpos[i]; - if (niterrun==niterec) { + if (recur==1) { /* use f1min from previous iteration */ for (j = 0; j < nts; j++) { Ni[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j]; } @@ -595,6 +595,7 @@ int main (int argc, char **argv) f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j]; } } + /* apply mute window */ for (j = nts-shift; j < nts; j++) { Ni[l*nxs*nts+i*nts+j] = 0.0; }