diff --git a/fdelmodc3D/Makefile b/fdelmodc3D/Makefile
index cb9e49d998d477e2350c36ae8ee5c12aead90fa3..e3f21d2bd2a8f7636a8bde0f7bf7ce8784177db4 100644
--- a/fdelmodc3D/Makefile
+++ b/fdelmodc3D/Makefile
@@ -25,6 +25,7 @@ SRCC = $(PRG).c \
acoustic6_3D.c \
acousticSH4_3D.c \
acoustic4_qr_3D.c \
+ elastic4dc_3D.c \
defineSource3D.c \
getParameters3D.c \
getWaveletInfo3D.c \
diff --git a/fdelmodc3D/acoustic4_3D.c b/fdelmodc3D/acoustic4_3D.c
index f3c2c346f37b0b935d0f32ab3d8c466bf3b16ade..36316d652c7b0bbb1c4006368d18d0ad13003a1b 100644
--- a/fdelmodc3D/acoustic4_3D.c
+++ b/fdelmodc3D/acoustic4_3D.c
@@ -6,22 +6,26 @@
#define MIN(x,y) ((x) < (y) ? (x) : (y))
-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);
-//long applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, long verbose);
+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);
-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);
-//long storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, 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);
-long reStoreSourceOnSurface3D(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);
-//long reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
+long reStoreSourceOnSurface3D(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);
-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);
-//long boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, 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);
-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 boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, 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);
-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)
+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)
{
/*********************************************************************
COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
@@ -82,9 +86,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
/* calculate vx for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
-//#pragma ivdep
for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
- #pragma ivdep //deletar
+ #pragma ivdep
for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]) +
@@ -96,9 +99,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
/* calculate vy for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
-//#pragma ivdep
for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
- #pragma ivdep //deletar
+ #pragma ivdep
for (iz=mod.ioYz; iz<mod.ieYz; iz++) {
vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]) +
@@ -110,9 +112,8 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
/* calculate vz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iy, iz) schedule(guided,1)
for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
-//#pragma ivdep
for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
- #pragma ivdep //deletar
+ #pragma ivdep
for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]) +
@@ -142,12 +143,10 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
/* calculate p/tzz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iy, iz) schedule(guided,1)
-//#pragma omp for private (ix, iz) schedule(dynamic)
#pragma ivdep
for (ix=mod.ioPx; ix<mod.iePx; ix++) {
-//#pragma ivdep
for (iy=mod.ioPy; iy<mod.iePy; iy++) {
- #pragma ivdep //deletar
+ #pragma ivdep
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
p[iy*n1*n2+ix*n1+iz] -= l2m[iy*n1*n2+ix*n1+iz]*(
c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) +
@@ -156,11 +155,6 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]) +
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]));
-
- //checking //deletar
- if(iz==(izsrc+mod.ioPz+bnd.ntap) && ix==(ixsrc+mod.ioPx+bnd.ntap) && iy==(iysrc+mod.ioPy+bnd.ntap) ){
- //vmess("inloop--p[ix=%li,iy=%li,iz=%li] at itime %li has value %e", ix, iy, iz, itime, p[iy*n1*n2+ix*n1+iz]); //deletar
- }
}
}
}
diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
index fa6437dab6aea1d48056c02a62228145728b9816..0a7d1f1165a4eb2d7c7f4c3bdcaa374fd205c1c0 100644
--- a/fdelmodc3D/applySource3D.c
+++ b/fdelmodc3D/applySource3D.c
@@ -20,14 +20,15 @@ void vmess(char *fmt, ...);
* The Netherlands
*
**********************************************************************/
-// applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
-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)
+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)
{
long is0, ibndz, ibndy, ibndx;
long isrc, ix, iy, iz, n1, n2;
long id1, id2, id3;
float src_ampl, time, scl, dt, sdx;
- float Mxx, Myy, Mzz, Mxz, Myz, Mxy, rake;
+ float Mxx, Myy, Mzz, Mxz, Myz, Mxy;
static long first=1;
if (src.type==6) {
@@ -350,19 +351,22 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
}
}
else if(src.type == 9) {
- rake = 0.5*M_PI;
- Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(rake)*sin(src.strike)*sin(src.strike));
- Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(src.strike)+cos(src.dip*2.0)*sin(rake)*sin(src.strike));
- Mzz = sin(src.dip*2.0)*sin(rake);
+ Mxx = -1.0*(sin(src.dip)*cos(src.rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(src.rake)*sin(src.strike)*sin(src.strike));
+ Myy = sin(src.dip)*cos(src.rake)*sin(2.0*src.strike)-sin(src.dip*2.0)*sin(src.rake)*cos(src.strike)*cos(src.strike);
+ Mzz = sin(src.dip*2.0)*sin(src.rake);
+ Mxz = -1.0*(cos(src.dip)*cos(src.rake)*cos(src.strike)+cos(src.dip*2.0)*sin(src.rake)*sin(src.strike));
+ Mxy = sin(src.dip)*cos(src.rake)*cos(src.strike*2.0)+0.5*(sin(src.dip*2.0)*sin(src.rake)*sin(src.strike*2.0));
+ Myz = -1.0*(cos(src.dip)*cos(src.rake)*sin(src.strike)-cos(src.dip*2.0)*sin(src.rake)*cos(src.strike));
txx[iy*n1*n2+ix*n1+iz] -= Mxx*src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz] -= Myy*src_ampl;
tzz[iy*n1*n2+ix*n1+iz] -= Mzz*src_ampl;
txz[iy*n1*n2+ix*n1+iz] -= Mxz*src_ampl;
+ txy[iy*n1*n2+ix*n1+iz] -= Mxy*src_ampl;
+ tyz[iy*n1*n2+ix*n1+iz] -= Myz*src_ampl;
} /* src.type */
} /* ischeme */
} /* loop over isrc */
-
-
return 0;
-}
+}
\ No newline at end of file
diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 180e774d6ae6be26f7d07f14bfd3aef92eaeb617..8281c66e6db3dc7b1e34c946381aeec9c3711839 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -6,7 +6,9 @@
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)
+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)
{
/*********************************************************************
Joeri's original filling order
@@ -77,7 +79,7 @@ MID left mid mid
dy = mod.dy;
dt = mod.dt;
fac = dt/dx;
- if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+ if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) || (bnd.fro==2) || (bnd.bac==2) ) pml=1;
else pml=0;
ibnd = mod.iorder/2-1;
@@ -102,7 +104,6 @@ MID left mid mid
if (bnd.top==3) { /* rigid surface at top */
#pragma omp for private (ix, iy, iz) nowait
-//#pragma ivdep
for (iy=1; iy<=ny; iy++) {
#pragma ivdep
for (ix=1; ix<=nx; ix++) {
@@ -116,7 +117,6 @@ MID left mid mid
}
if (bnd.rig==3) { /* rigid surface at right */
#pragma omp for private (ix, iy, iz) nowait
-//#pragma ivdep
for (iy=1; iy<=ny; iy++) {
#pragma ivdep
for (iz=1; iz<=nz; iz++) {
@@ -220,7 +220,6 @@ MID left mid mid
ibz = (bnd.ntap+izo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -242,7 +241,6 @@ MID left mid mid
ibx = (ixo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -263,7 +261,6 @@ MID left mid mid
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -285,7 +282,6 @@ MID left mid mid
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -315,7 +311,6 @@ MID left mid mid
ibx = (bnd.ntap+ixo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -335,7 +330,6 @@ MID left mid mid
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -357,7 +351,6 @@ MID left mid mid
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -383,7 +376,6 @@ MID left mid mid
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -406,7 +398,6 @@ MID left mid mid
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -433,7 +424,6 @@ MID left mid mid
ib = (bnd.ntap+izo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -455,7 +445,6 @@ MID left mid mid
ibx = (ixo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -474,7 +463,6 @@ MID left mid mid
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -497,7 +485,6 @@ MID left mid mid
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
#pragma ivdep
for (iz=izo; iz<ize; iz++) {
@@ -813,92 +800,109 @@ MID left mid mid
}
}
}
-
- }
- /*********/
- /* Bottom */
- /*********/
- if (bnd.bot==4) {
-
- if (mod.ischeme <= 2) { /* Acoustic scheme */
+ else { /* Elastic scheme */
/* Vx field */
- /* mid bottom mid vx*/
+ /* mid top mid vx */
ixo = mod.ioXx;
ixe = mod.ieXx;
iyo = mod.ioXy;
iye = mod.ieXy;
- izo = mod.ieXz;
- ize = mod.ieXz+bnd.ntap;
-
- ib = (ize-bnd.ntap);
+ izo = mod.ioXz-bnd.ntap;
+ ize = mod.ioXz;
+
+ ibz = (bnd.ntap+izo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[iz-ib];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ibz-iz];
}
}
}
- /* right bottom mid corner vx */
+
+
+ /* right top mid corner vx */
if (bnd.rig==4) {
ixo = mod.ieXx;
ixe = mod.ieXx+bnd.ntap;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
ibx = (ixo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
- /* right bottom front corner vx */
+
+
+ /* right top front corner vx */
if (bnd.fro==4) {
iyo = mod.ioXy-bnd.ntap;
iye = mod.ioXy;
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* right bottom back corner vx */
+
+
+ /* right top back corner vx */
if (bnd.bac==4) {
iyo = mod.ieXy;
iye = mod.ieXy+bnd.ntap;
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
@@ -909,187 +913,244 @@ MID left mid mid
ixe = mod.ieXx;
iyo = mod.ioXy;
iye = mod.ieXy;
- izo = mod.ieXz;
- ize = mod.ieXz+bnd.ntap;
+ izo = mod.ioXz-bnd.ntap;
+ ize = mod.ioXz;
- /* left bottom mid corner vx*/
+ /* left top mid corner vx */
if (bnd.lef==4) {
ixo = mod.ioXx-bnd.ntap;
ixe = mod.ioXx;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
ibx = (bnd.ntap+ixo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
- /* left bottom front corner vx */
+
+ /* left top front corner vx */
if (bnd.fro==4) {
iyo = mod.ioXy-bnd.ntap;
iye = mod.ioXy;
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* left bottom back corner vx*/
+
+
+ /* left top back corner vx */
if (bnd.bac==4) {
iyo = mod.ieXy;
iye = mod.ieXy+bnd.ntap;
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
}
- /* mid bottom front corner vx */
+
+
+ /* mid top front corner vx */
if (bnd.fro==4) {
ixo = mod.ioXx;
ixe = mod.ieXx;
iyo = mod.ioXy-bnd.ntap;
iye = mod.ioXy;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* mid bottom back corner vx */
+
+
+ /* mid top back corner vx */
if (bnd.bac==4) {
iyo = mod.ieXy;
iye = mod.ieXy+bnd.ntap;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
}
}
}
- }
-
+ }
+
/* Vy field */
- /* mid bottom mid vy */
+ /* mid top mid vy */
ixo = mod.ioYx;
ixe = mod.ieYx;
iyo = mod.ioYy;
iye = mod.ieYy;
- izo = mod.ieYz;
- ize = mod.ieYz+bnd.ntap;
+ izo = mod.ioYz-bnd.ntap;
+ ize = mod.ioYz;
- ib = (ize-bnd.ntap);
+ ib = (bnd.ntap+izo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iz-ib];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[ib-iz];
}
}
}
- /* right bottom mid corner vy */
+
+
+ /* right top mid corner vy */
if (bnd.rig==4) {
ixo = mod.ieYx;
- ixe = ixo+bnd.ntap;
- ibz = (izo);
+ ixe = mod.ieYx+bnd.ntap;
+ ibz = (bnd.ntap+izo-1);
ibx = (ixo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
- /* right bottom front corner vy */
+ /* right top front corner vy */
if (bnd.fro==4) {
iyo = mod.ioYy-bnd.ntap;
iye = mod.ioYy;
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* right bottom back corner vy */
+
+
+
+ /* right top back corner vy */
if (bnd.bac==4) {
iyo = mod.ieYy;
iye = mod.ieYy+bnd.ntap;
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-#pragma ivdep
+//#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
@@ -1099,14 +1160,15 @@ MID left mid mid
ixe = mod.ieYx;
iyo = mod.ioYy;
iye = mod.ieYy;
- izo = mod.ieYz;
- ize = mod.ieYz+bnd.ntap;
+ izo = mod.ioYz-bnd.ntap;
+ ize = mod.ioYz;
- /* left bottom mid corner vy */
+
+ /* left top mid corner vy */
if (bnd.lef==4) {
ixo = mod.ioYx-bnd.ntap;
ixe = mod.ioYx;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
ibx = (bnd.ntap+ixo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1114,14 +1176,19 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
- /* left bottom front corner vy */
+
+ /* left top front corner vy */
if (bnd.fro==4) {
iyo = mod.ioYy-bnd.ntap;
iye = mod.ioYy;
@@ -1132,15 +1199,20 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* left bottom back corner vy */
+
+ /* left top back corner vy */
if (bnd.bac==4) {
iyo = mod.ieYy;
iye = mod.ieYy+bnd.ntap;
@@ -1151,22 +1223,26 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
}
- /* mid bottom front corner vy */
+ /* front top mid corner vy */
if (bnd.fro==4) {
ixo = mod.ioYx;
ixe = mod.ieYx;
iyo = mod.ioYy-bnd.ntap;
iye = mod.ioYy;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1174,19 +1250,23 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* mid bottom back corner vy */
+ /* back top mid corner vy */
if (bnd.bac==4) {
iyo = mod.ieYy;
iye = mod.ieYy+bnd.ntap;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1194,10 +1274,14 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
}
}
}
@@ -1205,49 +1289,56 @@ MID left mid mid
/* Vz field */
+ /* mid top mid vz*/
ixo = mod.ioZx;
ixe = mod.ieZx;
iyo = mod.ioZy;
iye = mod.ieZy;
- izo = mod.ieZz;
- ize = mod.ieZz+bnd.ntap;
-
- ib = (ize-bnd.ntap);
+ izo = mod.ioZz-bnd.ntap;
+ ize = mod.ioZz;
+
+ ib = (bnd.ntap+izo-1);
#pragma omp for private (ix, iy, iz)
for (ix=ixo; ix<ixe; ix++) {
#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]*(
- 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];
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-iz];
}
}
}
-
-
- /* right bottom corner */
+ /* right top mid corner vz */
if (bnd.rig==4) {
ixo = mod.ieZx;
ixe = ixo+bnd.ntap;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
ibx = (ixo);
-#pragma omp for private(ix,iy,iz)
+#pragma omp for private(ix,iz)
for (ix=ixo; ix<ixe; ix++) {
#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]*(
- 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.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
-
- /* right bottom front corner */
+ /* right top front corner vz */
if (bnd.fro==4) {
iyo = mod.ioZy-bnd.ntap;
iye = mod.ioZy;
@@ -1258,18 +1349,19 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
}
-
-
-
- /* right bottom back corner */
+ /* right top back corner vz */
if (bnd.bac==4) {
iyo = mod.ieZy;
iye = mod.ieZy+bnd.ntap;
@@ -1280,63 +1372,74 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(ibz-iz)];
}
}
}
}
}
-
+
ixo = mod.ioZx;
ixe = mod.ieZx;
iyo = mod.ioZy;
iye = mod.ieZy;
- izo = mod.ieZz;
- ize = mod.ieZz+bnd.ntap;
+ izo = mod.ioZz-bnd.ntap;
+ ize = mod.ioZz;
- /* left bottom corner */
+ /* left top mid corner vz */
if (bnd.lef==4) {
ixo = mod.ioZx-bnd.ntap;
ixe = mod.ioZx;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
ibx = (bnd.ntap+ixo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
-//#pragma ivdep
+#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
- /* left bottom front corner */
+ /* left top front corner vz */
if (bnd.fro==4) {
iyo = mod.ioZy-bnd.ntap;
iye = mod.ioZy;
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
+#pragma ivdep
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* left bottom back corner */
+ /* left top back corner vz */
if (bnd.bac==4) {
iyo = mod.ieZy;
iye = mod.ieZy+bnd.ntap;
@@ -1347,22 +1450,26 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(ibz-iz)];
}
}
}
}
}
- /* front bottom corner */
+ /* mid top front corner vz */
if (bnd.fro==4) {
ixo = mod.ioZx;
ixe = mod.ieZx;
iyo = mod.ioZy-bnd.ntap;
iye = mod.ioZy;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (bnd.ntap+iyo-1);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1370,19 +1477,23 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibz-iz)];
}
}
}
}
- /* Back bottom corner */
+ /* mid top back corner vz */
if (bnd.bac==4) {
iyo = mod.ieZy;
iye = mod.ieZy+bnd.ntap;
- ibz = (izo);
+ ibz = (bnd.ntap+izo-1);
iby = (iyo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1390,37 +1501,1753 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
- vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibz-iz)];
}
}
}
}
-
}
+
}
-
-
/*********/
- /* Left */
+ /* Bottom */
+ /*********/
+ if (bnd.bot==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ /* mid bottom mid vx*/
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ieXz;
+ ize = mod.ieXz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[iz-ib];
+ }
+ }
+ }
+ /* right bottom mid corner vx */
+ if (bnd.rig==4) {
+ ixo = mod.ieXx;
+ ixe = mod.ieXx+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* right bottom front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* right bottom back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ieXz;
+ ize = mod.ieXz+bnd.ntap;
+
+ /* left bottom mid corner vx*/
+ if (bnd.lef==4) {
+ ixo = mod.ioXx-bnd.ntap;
+ ixe = mod.ioXx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner vx*/
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* mid bottom front corner vx */
+ if (bnd.fro==4) {
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* mid bottom back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+ /* Vy field */
+ /* mid bottom mid vy */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ieYz;
+ ize = mod.ieYz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iz-ib];
+ }
+ }
+ }
+ /* right bottom mid corner vy */
+ if (bnd.rig==4) {
+ ixo = mod.ieYx;
+ ixe = ixo+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* right bottom front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* right bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ieYz;
+ ize = mod.ieYz+bnd.ntap;
+
+ /* left bottom mid corner vy */
+ if (bnd.lef==4) {
+ ixo = mod.ioYx-bnd.ntap;
+ ixe = mod.ioYx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* mid bottom front corner vy */
+ if (bnd.fro==4) {
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* mid bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+ /* Vz field */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ieZz;
+ ize = mod.ieZz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private (ix, iy, iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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];
+ }
+ }
+ }
+
+
+ /* right bottom corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieZx;
+ ixe = ixo+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+
+ /* right bottom front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+
+ /* right bottom back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ieZz;
+ ize = mod.ieZz+bnd.ntap;
+
+ /* left bottom corner */
+ if (bnd.lef==4) {
+ ixo = mod.ioZx-bnd.ntap;
+ ixe = mod.ioZx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+//#pragma ivdep
+ 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]*(
+ 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.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+ 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]*(
+ 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.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* front bottom corner */
+ if (bnd.fro==4) {
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* Back bottom corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+ }
+ else { /* Elastic scheme */
+
+ /* Vx field */
+ /* mid bottom mid vx*/
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ieXz;
+ ize = mod.ieXz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[iz-ib];
+ }
+ }
+ }
+ /* right bottom mid corner vx */
+ if (bnd.rig==4) {
+ ixo = mod.ieXx;
+ ixe = mod.ieXx+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* right bottom front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* right bottom back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ieXz;
+ ize = mod.ieXz+bnd.ntap;
+
+ /* left bottom mid corner vx*/
+ if (bnd.lef==4) {
+ ixo = mod.ioXx-bnd.ntap;
+ ixe = mod.ioXx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner vx*/
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* mid bottom front corner vx */
+ if (bnd.fro==4) {
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* mid bottom back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+ /* Vy field */
+ /* mid bottom mid vy */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ieYz;
+ ize = mod.ieYz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iz-ib];
+ }
+ }
+ }
+ /* right bottom mid corner vy */
+ if (bnd.rig==4) {
+ ixo = mod.ieYx;
+ ixe = ixo+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* right bottom front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* right bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ieYz;
+ ize = mod.ieYz+bnd.ntap;
+
+ /* left bottom mid corner vy */
+ if (bnd.lef==4) {
+ ixo = mod.ioYx-bnd.ntap;
+ ixe = mod.ioYx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* mid bottom front corner vy */
+ if (bnd.fro==4) {
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* mid bottom back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+ /* Vz field */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ieZz;
+ ize = mod.ieZz+bnd.ntap;
+
+ ib = (ize-bnd.ntap);
+#pragma omp for private (ix, iy, iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[iz-ib];
+ }
+ }
+ }
+
+
+ /* right bottom corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieZx;
+ ixe = ixo+bnd.ntap;
+ ibz = (izo);
+ ibx = (ixo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+
+ /* right bottom front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+
+
+ /* right bottom back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ieZz;
+ ize = mod.ieZz+bnd.ntap;
+
+ /* left bottom corner */
+ if (bnd.lef==4) {
+ ixo = mod.ioZx-bnd.ntap;
+ ixe = mod.ioZx;
+ ibz = (izo);
+ ibx = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+//#pragma ivdep
+ 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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ /* left bottom front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+ 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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iby-iy)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* left bottom back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxyz[(iy-iby)*bnd.ntap*bnd.ntap+(ibx-ix)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ }
+ /* front bottom corner */
+ if (bnd.fro==4) {
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ ibz = (izo);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+ /* Back bottom corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ ibz = (izo);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ }
+ }
+
+ }
+ }
+
+
+
+ /*********/
+ /* Left */
+ /*********/
+ if (bnd.lef==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ /* left mid mid vx */
+ ixo = mod.ioXx-bnd.ntap;
+ ixe = mod.ioXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ioXz;
+ ize = mod.ieXz;
+
+ ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ /* Vy field */
+ /* left mid mid vy */
+ ixo = mod.ioYx-bnd.ntap;
+ ixe = mod.ioYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz;
+ ize = mod.ieYz;
+
+ ib = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]));
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ /* Vz field */
+ /* left mid mid vz */
+ ixo = mod.ioZx-bnd.ntap;
+ ixe = mod.ioZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ioZz;
+ ize = mod.ieZz;
+
+ ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private (ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vz*/
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vz */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ }
+
+ else { /* Elastic scheme */
+
+ /* Vx field */
+ /* left mid mid vx */
+ ixo = mod.ioXx-bnd.ntap;
+ ixe = mod.ioXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ioXz;
+ ize = mod.ieXz;
+
+ ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vx */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vx */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ /* Vy field */
+ /* left mid mid vy */
+ ixo = mod.ioYx-bnd.ntap;
+ ixe = mod.ioYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz;
+ ize = mod.ieYz;
+
+ ib = (bnd.ntap+ixo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vy */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vy */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ /* Vz field */
+ /* left mid mid vz */
+ ixo = mod.ioZx-bnd.ntap;
+ ixe = mod.ioZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ioZz;
+ ize = mod.ieZz;
+
+ ib = (bnd.ntap+ixo-1);
+
+#pragma omp for private (ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ib-ix];
+ }
+ }
+ }
+ /* left mid front corner vz*/
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+ /* left mid back corner vz */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ ibx = (bnd.ntap+ixo-1);
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ }
+ }
+ }
+ }
+
+ }
+
+ }
+
+ /*********/
+ /* Right */
/*********/
- if (bnd.lef==4) {
+ if (bnd.rig==4) {
if (mod.ischeme <= 2) { /* Acoustic scheme */
/* Vx field */
- /* left mid mid vx */
- ixo = mod.ioXx-bnd.ntap;
- ixe = mod.ioXx;
+ /* right mid mid vx */
+ ixo = mod.ieXx;
+ ixe = mod.ieXx+bnd.ntap;
iyo = mod.ioXy;
iye = mod.ieXy;
izo = mod.ioXz;
ize = mod.ieXz;
-
- ib = (bnd.ntap+ixo-1);
+
+ ib = (ixo);
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
@@ -1431,16 +3258,17 @@ MID left mid mid
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]));
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib];
}
}
}
- /* left mid front corner vx */
+ /* right mid front corner vx */
if (bnd.fro==4) {
iyo = mod.ioXy-bnd.ntap;
iye = mod.ioXy;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (bnd.ntap+iyo-1);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1450,17 +3278,18 @@ MID left mid mid
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]));
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
}
}
}
- /* left mid back corner vx */
+ /* right mid back corner vx*/
if (bnd.bac==4) {
iyo = mod.ieXy;
iye = mod.ieXy+bnd.ntap;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (iyo);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1470,22 +3299,23 @@ MID left mid mid
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]));
- vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
}
}
}
/* Vy field */
- /* left mid mid vy */
- ixo = mod.ioYx-bnd.ntap;
- ixe = mod.ioYx;
+ /* right mid mid vy */
+ ixo = mod.ieYx;
+ ixe = mod.ieYx+bnd.ntap;
iyo = mod.ioYy;
iye = mod.ieYy;
izo = mod.ioYz;
ize = mod.ieYz;
-
- ib = (bnd.ntap+ixo-1);
+
+ ib = (ixo);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1494,17 +3324,18 @@ MID left mid mid
vy[iy*n1*n2+ix*n1+iz] -= roy[iy*n1*n2+ix*n1+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]));
-
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ib-ix];
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib];
}
}
}
- /* left mid front corner vy */
+ /* right mid front corner vy */
if (bnd.fro==4) {
iyo = mod.ioYy-bnd.ntap;
iye = mod.ioYy;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (bnd.ntap+iyo-1);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1514,17 +3345,18 @@ MID left mid mid
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]));
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
}
}
}
- /* left mid back corner vy */
+ /* right mid back corner vy */
if (bnd.bac==4) {
iyo = mod.ieYy;
iye = mod.ieYy+bnd.ntap;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (iyo);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1534,24 +3366,25 @@ MID left mid mid
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]));
- vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
}
}
}
-
+
/* Vz field */
- /* left mid mid vz */
- ixo = mod.ioZx-bnd.ntap;
- ixe = mod.ioZx;
+ /* right mid mid vz */
+
+ ixo = mod.ieZx;
+ ixe = mod.ieZx+bnd.ntap;
iyo = mod.ioZy;
iye = mod.ieZy;
izo = mod.ioZz;
ize = mod.ieZz;
+
+ ib = (ixo);
- ib = (bnd.ntap+ixo-1);
-
-#pragma omp for private (ix,iy,iz)
+#pragma omp for private (ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
@@ -1559,17 +3392,18 @@ MID left mid mid
vz[iy*n1*n2+ix*n1+iz] -= roz[iy*n1*n2+ix*n1+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[ib-ix];
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ix-ib];
}
}
}
- /* left mid front corner vz*/
+ /* right mid front corner vz */
if (bnd.fro==4) {
iyo = mod.ioZy-bnd.ntap;
iye = mod.ioZy;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (bnd.ntap+iyo-1);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1579,17 +3413,18 @@ MID left mid mid
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.tapxz[(iby-iy)*bnd.ntap+(ibx-ix)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
}
}
}
- /* left mid back corner vz */
+ /* right mid back corner vz */
if (bnd.bac==4) {
iyo = mod.ieZy;
iye = mod.ieZy+bnd.ntap;
- ibx = (bnd.ntap+ixo-1);
+ ibx = (ixo);
iby = (iyo);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -1599,22 +3434,14 @@ MID left mid mid
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.tapxz[(iy-iby)*bnd.ntap+(ibx-ix)];
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
}
}
}
-
- }
-
- }
-
- /*********/
- /* Right */
- /*********/
- if (bnd.rig==4) {
- if (mod.ischeme <= 2) { /* Acoustic scheme */
+ }
+ else { /* Elastic scheme */
/* Vx field */
/* right mid mid vx */
@@ -1633,8 +3460,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
vx[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib];
}
@@ -1653,8 +3484,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
@@ -1674,8 +3509,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
vx[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
@@ -1700,11 +3539,14 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
vy[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib];
- //vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[ix-ib];
}
}
}
@@ -1721,8 +3563,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
@@ -1742,8 +3588,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
vy[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
@@ -1769,11 +3619,14 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
vz[iy*n1*n2+ix*n1+iz] *= bnd.tapz[ix-ib];
- //vz[iy*n1*n2+ix*n1+iz] *= bnd.tapx[ix-ib];
}
}
}
@@ -1790,8 +3643,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iby-iy)*bnd.ntap+(ix-ibx)];
}
@@ -1811,8 +3668,12 @@ MID left mid mid
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]*(
- 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]));
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
vz[iy*n1*n2+ix*n1+iz] *= bnd.tapxz[(iy-iby)*bnd.ntap+(ix-ibx)];
}
@@ -1903,6 +3764,93 @@ MID left mid mid
}
}
}
+
+ else { /* Elastic scheme */
+
+ /* Vx field */
+ /* mid mid front vx */
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ izo = mod.ioXz;
+ ize = mod.ieXz;
+
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iby-iy];
+ }
+ }
+ }
+
+ /* Vy field */
+ /* mid mid front vy */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ izo = mod.ioYz;
+ ize = mod.ieYz;
+
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iby-iy];
+ }
+ }
+ }
+
+ /* Vz field */
+ /* mid mid front vz */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ izo = mod.ioZz;
+ ize = mod.ieZz;
+
+ iby = (bnd.ntap+iyo-1);
+#pragma omp for private (ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iby-iy];
+ }
+ }
+ }
+ }
}
@@ -1985,6 +3933,93 @@ MID left mid mid
}
}
}
+
+ else { /* Elastic scheme */
+
+ /* Vx field */
+ /* mid mid back vx */
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ izo = mod.ioXz;
+ ize = mod.ieXz;
+
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+
+ vx[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iy-iby];
+ }
+ }
+ }
+
+ /* Vy field */
+ /* mid mid back vy */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ izo = mod.ioYz;
+ ize = mod.ieYz;
+
+ iby = (iyo);
+#pragma omp for private(ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vy[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iy-iby];
+ }
+ }
+ }
+
+ /* Vz field */
+ /* mid mid back vz */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ izo = mod.ioZz;
+ ize = mod.ieZz;
+
+ iby = (iyo);
+#pragma omp for private (ix,iy,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+
+ vz[iy*n1*n2+ix*n1+iz] *= bnd.tapy[iy-iby];
+ }
+ }
+ }
+ }
}
@@ -1993,6 +4028,7 @@ MID left mid mid
{
if (allocated) {
free(Vxpml);
+ free(Vypml);
free(Vzpml);
free(sigmu);
free(RA);
@@ -2039,7 +4075,7 @@ long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, floa
dy = mod.dy;
dt = mod.dt;
fac = dt/dx;
- if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+ if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) || (bnd.fro==2) || (bnd.bac==2) ) pml=1;
else pml=0;
/************************************************************/
diff --git a/fdelmodc3D/elastic4_3D.c b/fdelmodc3D/elastic4_3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..630ee02ac8c8d2499225963e2bbe6898a2c6a7be
--- /dev/null
+++ b/fdelmodc3D/elastic4_3D.c
@@ -0,0 +1,227 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+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);
+
+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);
+
+long reStoreSourceOnSurface3D(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);
+
+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);
+
+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 elastic4_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)
+{
+/*********************************************************************
+ COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
+
+ The captial symbols T (=Txx,Tzz) Txz,Vx,Vz represent the actual grid
+ The indices ix,iz are related to the T grid, so the capital T
+ symbols represent the actual modelled grid.
+
+ one cel (iz,ix)
+ |
+ V extra column of vx,txz
+ |
+ ------- V
+ | txz vz| txz vz txz vz txz vz txz vz txz vz txz
+ | |
+ | vx t | vx t vx t vx t vx t vx t vx
+ -------
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz--Txz-Vz--Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx t vx t vx t vx t vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz <--|
+ |
+ extra row of txz/vz |
+
+ AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+
+***********************************************************************/
+
+ float c1, c2;
+ float dvx, dvy, dvz;
+ long ix, iy, iz;
+ long n1, n2;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ n1 = mod.naz;
+ n2 = mod.nax;
+
+ /* calculate vx for all grid points except on the virtual boundary*/
+#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
+ for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+ }
+ }
+ }
+
+ /* calculate vy for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+ }
+ }
+ }
+
+ /* calculate vz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+ }
+ }
+ }
+
+ /* Add force source */
+ if (src.type > 5) {
+ applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, src_nwav, verbose);
+ }
+
+ /* boundary condition clears velocities on boundaries */
+ boundariesP3D(mod, bnd, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, itime, verbose);
+
+ /* calculate Txx/tzz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz, dvx, dvy, dvz) nowait schedule(guided,1)
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ dvx = 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]);
+ dvy = c1*(vy[(iy+1)*n2*n1+ix*n1+iz] - vy[iy*n2*n1+ix*n1+iz]) +
+ 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);
+ }
+ }
+ }
+
+
+
+ /* calculate Txz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioTy; iy<mod.ieTy; iy++) {
+ for (ix=mod.ioTx; ix<mod.ieTx; ix++) {
+ #pragma ivdep
+ for (iz=mod.ioTz; iz<mod.ieTz; iz++) {
+ txz[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*(
+ c1*(vx[iy*n2*n1+ix*n1+iz] - vx[iy*n2*n1+ix*n1+iz-1] +
+ vz[iy*n2*n1+ix*n1+iz] - vz[iy*n2*n1+(ix-1)*n1+iz]) +
+ c2*(vx[iy*n2*n1+ix*n1+iz+1] - vx[iy*n2*n1+ix*n1+iz-2] +
+ vz[iy*n2*n1+(ix+1)*n1+iz] - vz[iy*n2*n1+(ix-2)*n1+iz]) );
+ }
+ }
+ }
+
+ /* calculate Txy for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioTy; iy<mod.ieTy; iy++) {
+ for (ix=mod.ioTx; ix<mod.ieTx; ix++) {
+ #pragma ivdep
+ for (iz=mod.ioTz; iz<mod.ieTz; iz++) {
+ txy[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*(
+ c1*(vx[iy*n2*n1+ix*n1+iz] - vx[(iy-1)*n2*n1+ix*n1+iz] +
+ vy[iy*n2*n1+ix*n1+iz] - vy[iy*n2*n1+(ix-1)*n1+iz]) +
+ c2*(vx[(iy+1)*n2*n1+ix*n1+iz] - vx[(iy-2)*n2*n1+ix*n1+iz] +
+ vy[iy*n2*n1+(ix+1)*n1+iz] - vy[iy*n2*n1+(ix-2)*n1+iz]) );
+ }
+ }
+ }
+
+ /* calculate Tyz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioTy; iy<mod.ieTy; iy++) {
+ for (ix=mod.ioTx; ix<mod.ieTx; ix++) {
+ #pragma ivdep
+ for (iz=mod.ioTz; iz<mod.ieTz; iz++) {
+ tyz[iy*n2*n1+ix*n1+iz] -= mul[iy*n2*n1+ix*n1+iz]*(
+ c1*(vz[iy*n2*n1+ix*n1+iz] - vz[(iy-1)*n2*n1+ix*n1+iz] +
+ vy[iy*n2*n1+ix*n1+iz] - vy[iy*n2*n1+ix*n1+iz-1]) +
+ c2*(vz[(iy+1)*n2*n1+ix*n1+iz] - vz[(iy-2)*n2*n1+ix*n1+iz] +
+ vy[iy*n2*n1+ix*n1+iz+1] - vy[iy*n2*n1+ix*n1+iz-2]) );
+ }
+ }
+ }
+
+ /* Add stress source */
+ if (src.type < 6) {
+ applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, src_nwav, verbose);
+ }
+
+ /* check if there are sources placed on the boundaries */
+ storeSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, verbose);
+
+ /* Free surface: calculate free surface conditions for stresses */
+ boundariesV3D(mod, bnd, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, itime, verbose);
+
+ /* restore source positions on the edge */
+ 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/elastic4dc_3D.c b/fdelmodc3D/elastic4dc_3D.c
new file mode 100644
index 0000000000000000000000000000000000000000..4a86c78ea9b55b97fd81c476ad4c1a2923f6b4f2
--- /dev/null
+++ b/fdelmodc3D/elastic4dc_3D.c
@@ -0,0 +1,180 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+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);
+
+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);
+
+long reStoreSourceOnSurface3D(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);
+
+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);
+
+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 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)
+{
+/*********************************************************************
+ COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
+
+ The captial symbols T (=Txx,Tzz) Txz,Vx,Vz represent the actual grid
+ The indices ix,iz are related to the T grid, so the capital T
+ symbols represent the actual modelled grid.
+
+ one cel (iz,ix)
+ |
+ V extra column of vx,txz
+ |
+ ------- V
+ | txz vz| txz vz txz vz txz vz txz vz txz vz txz
+ | |
+ | vx t | vx t vx t vx t vx t vx t vx
+ -------
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz--Txz-Vz--Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+ | | | | | | |
+ txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
+ | | | | | | |
+ vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz
+
+ vx t vx t vx t vx t vx t vx t vx
+
+ txz vz txz vz txz vz txz vz txz vz txz vz txz <--|
+ |
+ extra row of txz/vz |
+
+ AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+
+***********************************************************************/
+
+ float c1, c2;
+ float dvx, dvy, dvz;
+ long ix, iy, iz;
+ long n1, n2;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ n1 = mod.naz;
+ n2 = mod.nax;
+
+ /* calculate vx for all grid points except on the virtual boundary*/
+#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
+ for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+ }
+ }
+ }
+
+ /* calculate vy for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+ }
+ }
+ }
+
+ /* calculate vz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
+ for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
+ 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]*(
+ 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]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+ }
+ }
+ }
+
+ /* Add force source */
+ if (src.type > 5) {
+ applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, src_nwav, verbose);
+ }
+
+ /* boundary condition clears velocities on boundaries */
+ boundariesP3D(mod, bnd, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, itime, verbose);
+
+ /* calculate Txx/tzz for all grid points except on the virtual boundary */
+#pragma omp for private (ix, iy, iz, dvx, dvy, dvz) nowait schedule(guided,1)
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+#pragma ivdep
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ dvx = 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]);
+ dvy = c1*(vy[(iy+1)*n2*n1+ix*n1+iz] - vy[iy*n2*n1+ix*n1+iz]) +
+ 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);
+ }
+ }
+ }
+
+ /* Add stress source */
+ if (src.type < 6) {
+ applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, src_nwav, verbose);
+ }
+
+ /* check if there are sources placed on the boundaries */
+ storeSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, verbose);
+
+ /* Free surface: calculate free surface conditions for stresses */
+ boundariesV3D(mod, bnd, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, itime, verbose);
+
+ /* restore source positions on the edge */
+ 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 b/fdelmodc3D/fdelmodc3D
deleted file mode 100755
index d167a03b98e669688a662a2490e725c3ec7b92ed..0000000000000000000000000000000000000000
Binary files a/fdelmodc3D/fdelmodc3D and /dev/null differ
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index 5db4c61fd9bc183e21a60b84c116a9464e675006..21b3f5858add01f7eca9c2d5d8597fbbd27926c0 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -48,6 +48,11 @@ 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);
+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);
+
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,
@@ -492,7 +497,6 @@ int main(int argc, char **argv)
shot.z[ishot]=iz+src.sinkdepth;
}
-
/* scan for free surface boundary in case it has a topography */
for (iy=0; iy<mod.ny; iy++) {
for (ix=0; ix<mod.nx; ix++) {
@@ -642,7 +646,8 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
vmess("Visco-Elastic order 4 not yet available");
break;
case 5 : /* Elastic FD kernel with S-velocity set to zero*/
- vmess("DC-Elastic order 4 not yet available");
+ elastic4dc_3D(mod, src, wav, bnd, it, ixsrc, iysrc, izsrc, src_nwav,
+ vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, verbose);
break;
}
diff --git a/fdelmodc3D/fdelmodc3D.h b/fdelmodc3D/fdelmodc3D.h
index fff9d9592c7804123079602d86c0e0c795f7e403..ed4ce7b49a22fd040aeb2530f1f925546814677e 100644
--- a/fdelmodc3D/fdelmodc3D.h
+++ b/fdelmodc3D/fdelmodc3D.h
@@ -179,6 +179,7 @@ typedef struct _sourcePar { /* Source Array Parameters */
float amplitude;
float dip;
float strike;
+ float rake;
long distribution;
long window;
long injectionrate;
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index fcbfe4a8f0a70be2937b26d9bb1b3a97dde59604..836ce4cf0ed7e9806a9fde628031286f2caf64d6 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -551,10 +551,11 @@ criteria we have imposed.*/
if (!getparfloat("dyshot",&dyshot)) dyshot=dy;
if (!getparfloat("dzshot",&dzshot)) dzshot=0.0;
if (!getparfloat("dip",&src->dip)) src->dip=0.0;
- if (!getparfloat("strike",&src->strike)) src->strike=1.0;
- if (src->strike>=0) src->strike=0.5*M_PI;
- else src->strike = -0.5*M_PI;
+ if (!getparfloat("strike",&src->strike)) src->strike=0.0;
+ if (!getparfloat("rake",&src->rake)) src->rake=0.0;
src->dip = M_PI*(src->dip/180.0);
+ src->strike = M_PI*(src->strike/180.0);
+ src->rake = M_PI*(src->rake/180.0);
if (shot->n>1) {
idxshot=MAX(0,NINT(dxshot/dx));
@@ -969,6 +970,7 @@ criteria we have imposed.*/
vmess("*******************************************");
vmess("wav_nt = %6li wav_nx = %li", wav->ns, wav->nx);
vmess("src_type = %6li src_orient = %li", src->type, src->orient);
+ vmess("number of sources = %li", shot->n);
vmess("fmax = %8.2f", fmax);
fprintf(stderr," %s: Source type : ",xargv[0]);
switch ( src->type ) {
diff --git a/fdelmodc3D/getRecTimes3D.c b/fdelmodc3D/getRecTimes3D.c
index fa007f88322b1cb038b7d2d820218d118ffe7835..4c19786a8028eb10c0bc3183fe3c31c72e3603e8 100644
--- a/fdelmodc3D/getRecTimes3D.c
+++ b/fdelmodc3D/getRecTimes3D.c
@@ -402,6 +402,7 @@ long getRecTimes3D(modPar mod, recPar rec, bndPar bnd, long itime, long isam,
}
}
if (rec.type.txx) rec_txx[irec*rec.nt+isam] = txx[iy*n1*n2+ix*n1+iz];
+ if (rec.type.tyy) rec_tyy[irec*rec.nt+isam] = tyy[iy*n1*n2+ix*n1+iz];
if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[iy*n1*n2+ix*n1+iz];
if (rec.type.txz) { /* time interpolation to be done */
if (rec.int_vz == 2 || rec.int_vx == 2) {
diff --git a/fdelmodc3D/readModel3D.c b/fdelmodc3D/readModel3D.c
index df93d7185f08428f42294fb63779a08a6e8daf9a..157f7a1790d4bcb1c8539aebf37385b170679e63 100644
--- a/fdelmodc3D/readModel3D.c
+++ b/fdelmodc3D/readModel3D.c
@@ -27,7 +27,7 @@
long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
- float *l2m, float *lam, float *muu, float *tss, float *tes, float *tep)
+ float *l2m, float *lam, float *muxz, float *tss, float *tes, float *tep)
{
FILE *fpcp, *fpcs, *fpro;
FILE *fpqp=NULL, *fpqs=NULL;
@@ -36,7 +36,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
long n1, n2, n3, ix, iy, iz, nz, ny, nx;
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, cs11, cs12, cs21, cs22, mul, mu, lamda2mu, lamda;
+ 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 a, b;
@@ -69,7 +69,7 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
ioPx=mod.ioPx;
ioPy=mod.ioPy;
ioPz=mod.ioPz;
- /* Txz, Txy, Tyz,: muu */
+ /* Txz, Txy, Tyz,: muxz */
ioTx=mod.ioTx;
ioTy=mod.ioTy;
ioTz=mod.ioTz;
@@ -231,23 +231,263 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
/* the edges of the model */
if (mod.ischeme>2) { /* Elastic Scheme */
+ iz = nz-1;
+ for (iy=0;iy<ny-1;iy++) {
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+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;
+ }
+ }
+
+ iy = ny-1;
+ 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+1)*nz+iz];
+ cs212 = cs2c*ro[iy*nx*nz+(ix+1)*nz+iz+1];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+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;
+ }
+ }
+
+ ix = nx-1;
+ for (iy=0;iy<ny-1;iy++) {
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+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;
+ }
+ }
+
iz = nz-1;
iy = ny-1;
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];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs21 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
- cs22 = cs2a*ro[iy*nx*nz+(ix+1)*nz+iz];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
lamda = lamda2mu - 2*mu;
@@ -260,58 +500,124 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
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;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix = nx-1;
iz = nz-1;
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];
- cs2b = cs[(iy+1)*nx*nz+ix*nz+iz]*cs[(iy+1)*nx*nz+ix*nz+iz];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2b*ro[iy*nx*nz+ix*nz+iz];
- cs21 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs22 = cs2b*ro[iy*nx*nz+ix*nz+iz];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+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]);
+ 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;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix = nx-1;
iy = ny-1;
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];
- cs2b = cs[iy*nx*nz+ix*nz+iz+1]*cs[iy*nx*nz+ix*nz+iz+1];
- cs11 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs12 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
- cs21 = cs2*ro[iy*nx*nz+ix*nz+iz];
- cs22 = cs2b*ro[iy*nx*nz+ix*nz+iz+1];
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
}
else {
mul = 0.0;
}
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
mu = cs2*ro[iy*nx*nz+ix*nz+iz];
lamda2mu = cp2*ro[iy*nx*nz+ix*nz+iz];
lamda = lamda2mu - 2*mu;
@@ -324,64 +630,100 @@ long readModel3D(modPar mod, bndPar bnd, float *rox, float *roy, float *roz,
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;
- muu[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
+ muxz[(iy+ioTy)*n2*n1+(ix+ioTx)*n1+iz+ioTz]=fac*mul;
}
ix=nx-1;
+ iy=ny-1;
iz=nz-1;
- cp2 = cp[ix*nz+iz]*cp[ix*nz+iz];
- cs2 = cs[ix*nz+iz]*cs[ix*nz+iz];
- mu = cs2*ro[ix*nz+iz];
- lamda2mu = cp2*ro[ix*nz+iz];
+ 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];
lamda = lamda2mu - 2*mu;
- bx = ro[ix*nz+iz];
- bz = ro[ix*nz+iz];
- rox[(ix+ioXx)*n1+iz+ioXz]=fac/bx;
- roz[(ix+ioZx)*n1+iz+ioZz]=fac/bz;
- l2m[(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
- lam[(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(ix+ioTx)*n1+iz+ioTz]=fac*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;
+
+ 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];
+ if (cs111 > 0.0) {
+ mul = 4.0/(1.0/cs111+1.0/cs112+1.0/cs211+1.0/cs212);
+ }
+ else {
+ mul = 0.0;
+ }
- for (ix=0;ix<nx-1;ix++) {
- for (iz=0;iz<nz-1;iz++) {
- cp2 = cp[ix*nz+iz]*cp[ix*nz+iz];
- cs2 = cs[ix*nz+iz]*cs[ix*nz+iz];
- cs2a = cs[(ix+1)*nz+iz]*cs[(ix+1)*nz+iz];
- cs2b = cs[ix*nz+iz+1]*cs[ix*nz+iz+1];
- cs2c = cs[(ix+1)*nz+iz+1]*cs[(ix+1)*nz+iz+1];
-
-/*
-Compute harmonic average of mul for accurate and stable fluid-solid interface
-see Finite-difference modeling of wave propagation in a fluid-solid configuration
-Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
-*/
-
- cs11 = cs2*ro[ix*nz+iz];
- cs12 = cs2b*ro[ix*nz+iz+1];
- cs21 = cs2a*ro[ix*nz+iz];
- cs22 = cs2c*ro[ix*nz+iz+1];
-// cpx = 0.5*(cp[ix*nz+iz]+cp[(ix+1)*nz+iz])
-// cpz = 0.5*(cp[ix*nz+iz]+cp[ix*nz+iz+1])
-
- if (cs11 > 0.0) {
- mul = 4.0/(1.0/cs11+1.0/cs12+1.0/cs21+1.0/cs22);
- }
- else {
- mul = 0.0;
- }
- mu = cs2*ro[ix*nz+iz];
- lamda2mu = cp2*ro[ix*nz+iz];
- lamda = lamda2mu - 2*mu; /* could also use mul to calculate lambda, but that might not be correct: question from Chaoshun Hu. Note use mu or mul as well on boundaries */
-
- bx = 0.5*(ro[ix*nz+iz]+ro[(ix+1)*nz+iz]);
- bz = 0.5*(ro[ix*nz+iz]+ro[ix*nz+iz+1]);
- rox[(ix+ioXx)*n1+iz+ioXz]=fac/bx;
- roz[(ix+ioZx)*n1+iz+ioZz]=fac/bz;
- l2m[(ix+ioPx)*n1+iz+ioPz]=fac*lamda2mu;
- lam[(ix+ioPx)*n1+iz+ioPz]=fac*lamda;
- muu[(ix+ioTx)*n1+iz+ioTz]=fac*mul;
- }
- }
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs121+1.0/cs112+1.0/cs122);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ /* 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];
+ // if (cs111 > 0.0) {
+ // mul = 4.0/(1.0/cs111+1.0/cs211+1.0/cs121+1.0/cs221);
+ // }
+ // else {
+ // mul = 0.0;
+ // }
+
+ mu = cs2*ro[iy*nx*nz+ix*nz+iz];
+ lamda2mu = cp2*ro[iy*nx*nz+ix*nz+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;
+ }
+ }
+ }
}
else { /* Acoustic Scheme */
@@ -607,7 +949,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ioTx+bnd.ntap;
iyo = mod.ioTy;
@@ -617,7 +959,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ixe*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ixe*n1+iz];
}
}
}
@@ -738,7 +1080,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ieTx-bnd.ntap;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -748,7 +1090,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+(ixo-1)*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+(ixo-1)*n1+iz];
}
}
}
@@ -880,7 +1222,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
if (bnd.lef==4 || bnd.lef==2) ixo -= bnd.ntap;
@@ -892,7 +1234,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iye*n2*n1+ix*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iye*n2*n1+ix*n1+iz];
}
}
}
@@ -1029,7 +1371,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
if (bnd.lef==4 || bnd.lef==2) ixo -= bnd.ntap;
@@ -1041,7 +1383,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[(iyo-1)*n2*n1+ix*n1+iz];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[(iyo-1)*n2*n1+ix*n1+iz];
}
}
}
@@ -1178,7 +1520,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu field */
+ /* muxz field */
ixo = mod.ioTx;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -1188,7 +1530,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ix*n1+ize];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+ize];
}
}
}
@@ -1320,7 +1662,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
}
- /* muu */
+ /* muxz */
ixo = mod.ioTx;
ixe = mod.ieTx;
iyo = mod.ioTy;
@@ -1330,7 +1672,7 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
for (iy=iyo; iy<iye; iy++) {
for (ix=ixo; ix<ixe; ix++) {
for (iz=izo; iz<ize; iz++) {
- muu[iy*n2*n1+ix*n1+iz] = muu[iy*n2*n1+ix*n1+izo-1];
+ muxz[iy*n2*n1+ix*n1+iz] = muxz[iy*n2*n1+ix*n1+izo-1];
}
}
}
@@ -1371,16 +1713,6 @@ Robbert van Vossen, Johan O. A. Robertsson, and Chris H. Chapman
}
-
- FILE *fptest;
- fptest = fopen("cp-test.bin","w");
- fwrite(cp, nz*ny*nx, sizeof(float), fptest);
- fclose(fptest);
-
- fptest = fopen("ro-test.bin","w");
- fwrite(ro, nz*ny*nx, sizeof(float), fptest);
- fclose(fptest);
-
free(cp);
free(ro);
free(cs);