From 40b281c5f78b0f7e6cef7292708c32b17f654bf8 Mon Sep 17 00:00:00 2001
From: "joeri.brackenhoff" <joeri.brackenhoff@login0.ogun.local>
Date: Thu, 14 Mar 2019 10:47:56 -0300
Subject: [PATCH] 3D
---
fdelmodc3D/acoustic4_3D.c | 174 +++
fdelmodc3D/applySource3D.c | 365 ++++++
fdelmodc3D/boundaries3D.c | 2253 +++++++++++++++++++++++++++++++++
fdelmodc3D/fdelmodc3D.c | 32 +-
fdelmodc3D/getRecTimes3D.c | 324 +++++
fdelmodc3D/writeSrcRecPos3D.c | 160 +++
fdelmodc3D/writesufile3D.c | 162 +++
7 files changed, 3449 insertions(+), 21 deletions(-)
create mode 100644 fdelmodc3D/acoustic4_3D.c
create mode 100644 fdelmodc3D/applySource3D.c
create mode 100644 fdelmodc3D/boundaries3D.c
create mode 100644 fdelmodc3D/getRecTimes3D.c
create mode 100644 fdelmodc3D/writeSrcRecPos3D.c
create mode 100644 fdelmodc3D/writesufile3D.c
diff --git a/fdelmodc3D/acoustic4_3D.c b/fdelmodc3D/acoustic4_3D.c
new file mode 100644
index 0000000..453d59f
--- /dev/null
+++ b/fdelmodc3D/acoustic4_3D.c
@@ -0,0 +1,174 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+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 storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, 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 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 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 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:
+
+ The captial symbols T (=P) 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;
+ long ix, iy, iz;
+ long n1, n2;
+ long ioXx, ioXy, ioXz, ioYx, ioYy, ioYz, ioZx, ioZy, ioZz, ioPx, ioPy, ioPz;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ n1 = mod.naz;
+ n2 = mod.nay;
+
+ /* 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++) {
+ 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]) +
+ c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz]));
+ }
+ }
+ }
+
+ /* 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++) {
+ 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]) +
+ c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*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++) {
+#pragma ivdep
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+ 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]) +
+ c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2]));
+ }
+ }
+ }
+
+ /* boundary condition clears velocities on boundaries */
+ boundariesP3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
+
+ /* Add force source */
+ if (src.type > 5) {
+ 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);
+ }
+
+ /* this is needed because the P fields are not using tapered boundaries (bnd....=4) */
+ if (bnd.top==2) mod.ioPz += bnd.npml;
+ if (bnd.bot==2) mod.iePz -= bnd.npml;
+ if (bnd.fro==2) mod.ioPy += bnd.npml;
+ if (bnd.bac==2) mod.iePy -= bnd.npml;
+ if (bnd.lef==2) mod.ioPx += bnd.npml;
+ if (bnd.rig==2) mod.iePx -= bnd.npml;
+
+ /* 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++) {
+ 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]) +
+ c2*(vx[iy*n1*n2+(ix+2)*n1+iz] - vx[iy*n1*n2+(ix-1)*n1+iz]) +
+ c1*(vy[(iy+1)*n1*n2+ix*n1+iz] - vy[iy*n1*n2+ix*n1+iz]) +
+ 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]));
+ }
+ }
+ }
+ if (bnd.top==2) mod.ioPz -= bnd.npml;
+ if (bnd.bot==2) mod.iePz += bnd.npml;
+ if (bnd.fro==2) mod.ioPy -= bnd.npml;
+ if (bnd.bac==2) mod.iePy += bnd.npml;
+ if (bnd.lef==2) mod.ioPx -= bnd.npml;
+ if (bnd.rig==2) mod.iePx += bnd.npml;
+
+ /* Add stress source */
+ if (src.type < 6) {
+ 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);
+ }
+
+/* Free surface: calculate free surface conditions for stresses */
+
+ /* check if there are sources placed on the free surface */
+ storeSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
+
+ /* Free surface: calculate free surface conditions for stresses */
+ boundariesV3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
+
+ /* restore source positions on the edge */
+ reStoreSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
+
+ return 0;
+}
diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
new file mode 100644
index 0000000..80298b3
--- /dev/null
+++ b/fdelmodc3D/applySource3D.c
@@ -0,0 +1,365 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+void vmess(char *fmt, ...);
+
+#define c1 (9.0/8.0)
+#define c2 (-1.0/24.0)
+
+/*********************************************************************
+ *
+ * Add's the source amplitude(s) to the grid.
+ *
+ * For the acoustic schemes, the source-type must not be txx tzz or txz.
+ *
+ * AUTHOR:
+ * Jan Thorbecke (janth@xs4all.nl)
+ * The Netherlands
+ *
+ **********************************************************************/
+
+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;
+ static long first=1;
+
+ if (src.type==6) {
+ ibndz = mod.ioXz;
+ ibndy = mod.ioXy;
+ ibndx = mod.ioXx;
+ }
+ else if (src.type==7) {
+ ibndz = mod.ioZz;
+ ibndy = mod.ioZy;
+ ibndx = mod.ioZx;
+ }
+ else if (src.type==2) {
+ ibndz = mod.ioTz;
+ ibndy = mod.ioTy;
+ ibndx = mod.ioTx;
+ if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+ if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+ if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+ }
+ else {
+ ibndz = mod.ioPz;
+ ibndy = mod.ioPy;
+ ibndx = mod.ioPx;
+ if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+ if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+ if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+ }
+
+ n1 = mod.naz;
+ n2 = mod.nax;
+ dt = mod.dt;
+ sdx = 1.0/mod.dx;
+
+ /* special txz source activated? */
+
+ if ((bnd.top==1) && (src.type==2)) {
+ iz = izsrc + ibndz;
+ if (iz==ibndz) {
+ if (src.orient != 1) {
+ if (first) {
+ vmess("Only monopole Txz source allowed at surface. Reset to monopole");
+ first = 0;
+ }
+ src.orient=1;
+ }
+ }
+ }
+
+/*
+* for plane wave sources the sources are placed
+* around the central shot position
+* the first source position has an offset in x of is0
+*
+* itime = 0 corresponds with time=0
+* itime = 1 corresponds with time=dt
+* src[0] (the first sample) corresponds with time = 0
+*/
+
+ is0 = -1*floor((src.n-1)/2);
+#pragma omp for private (isrc, src_ampl, ix, iy, iz, time, id1, id2, id3, scl)
+ for (isrc=0; isrc<src.n; isrc++) {
+ src_ampl=0.0;
+ /* calculate the source position */
+ if (src.random || src.multiwav) {
+ ix = src.x[isrc] + ibndx;
+ iy = src.y[isrc] + ibndy;
+ iz = src.z[isrc] + ibndz;
+ }
+ else { /* plane wave and point sources */
+ ix = ixsrc + ibndx + is0 + isrc;
+ iy = iysrc + ibndy + is0 + isrc;
+ iz = izsrc + ibndz;
+ }
+ time = itime*dt - src.tbeg[isrc];
+ id1 = floor(time/dt);
+ id2 = id1+1;
+
+ /* delay not reached or no samples left in source wavelet? */
+ if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
+
+// fprintf(stderr,"isrc=%li ix=%li iz=%li src.x=%li src.z=%li\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
+
+ if (!src.multiwav) { /* only one wavelet for all sources */
+ src_ampl = src_nwav[0][id1]*(id2-time/dt) + src_nwav[0][id2]*(time/dt-id1);
+ }
+ else { /* multi-wavelet sources */
+ src_ampl = src_nwav[isrc][id1]*(id2-time/dt) + src_nwav[isrc][id2]*(time/dt-id1);
+ }
+
+ if (src_ampl==0.0) continue;
+ if ( ((ix-ibndx)<0) || ((ix-ibndx)>mod.nx) ) continue; /* source outside grid */
+ if ( ((iy-ibndy)<0) || ((iy-ibndy)>mod.ny) ) continue; /* source outside grid */
+
+ if (verbose>=4 && itime==0) {
+ vmess("Source %li positioned at grid ix=%li iy=%li iz=%li",isrc, ix, iy, iz);
+ }
+
+ /* cosine squared windowing to reduce edge effects on shot arrays */
+ if ( (src.n>1) && src.window) {
+ scl = 1.0;
+ if (isrc < src.window) {
+ scl = cos(0.5*M_PI*(src.window - isrc)/src.window);
+ }
+ else if (isrc > src.n-src.window+1) {
+ scl = cos(0.5*M_PI*(src.window - (src.n-isrc+1))/src.window);
+ }
+ src_ampl *= scl*scl;
+ }
+
+ /* source scaling factor to compensate for discretisation */
+
+ /* old amplitude setting does not obey reciprocity */
+ // src_ampl *= rox[ix*n1+iz]*l2m[ix*n1+iz]/(dt);
+
+/* in older version added factor 2.0 to be compliant with defined Green's functions in Marchenko algorithm */
+/* this is now set to 1.0 */
+ src_ampl *= (1.0/mod.dx)*l2m[iy*n1*n2+ix*n1+iz];
+
+ if (verbose>5) {
+ vmess("Source %li at grid [ix=%li,iy=%li,iz=%li] at itime %li has value %e",isrc, ix, iy, iz, itime, src_ampl);
+ }
+
+ /* Force source */
+
+ if (src.type == 6) {
+ vx[iy*n1*n2+ix*n1+iz] += src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+ /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+ //vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
+ }
+ else if (src.type == 7) {
+ vz[iy*n1*n2+ix*n1+iz] += src_ampl*roz[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+ /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+ /* stable implementation changes amplitude and more work is needed */
+ //vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+ //vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+ } /* src.type */
+
+
+ /* Stress source */
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+ /* Compressional source */
+ if (src.type == 1) {
+ if (src.orient != 1) src_ampl=src_ampl/mod.dx;
+
+ if (src.orient==1) { /* monopole */
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ }
+ else if (src.orient==2) { /* dipole +/- */
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+ }
+ else if (src.orient==3) { /* dipole - + */
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+ }
+ else if (src.orient==4) { /* dipole +/0/- */
+ if (iz > ibndz)
+ tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+ if (iz < mod.nz+ibndz-1)
+ tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+ }
+ else if (src.orient==5) { /* dipole + - */
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+ }
+ }
+ }
+ else { /* Elastic scheme */
+ /* Compressional source */
+ if (src.type == 1) {
+ if (src.orient==1) { /* monopole */
+ txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ }
+ else if (src.orient==2) { /* dipole +/- */
+ txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ txx[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+ }
+ else if (src.orient==3) { /* dipole - + */
+ txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ txx[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+ tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+ }
+ else if (src.orient==4) { /* dipole +/0/- */
+ if (iz > ibndz) {
+ txx[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+ }
+ if (iz < mod.nz+ibndz-1) {
+ txx[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+ }
+ }
+ else if (src.orient==5) { /* dipole + - */
+ txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ txx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+ tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+ }
+ }
+ else if (src.type == 2) {
+ /* Txz source */
+ if ((iz == ibndz) && bnd.top==1) {
+ txz[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl;
+ txz[iy*n1*n2+ix*n1+iz-1] += src_ampl;
+ }
+ else {
+ txz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ }
+ /* possible dipole orientations for a txz source */
+ if (src.orient == 2) { /* dipole +/- */
+ txz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+ }
+ else if (src.orient == 3) { /* dipole - + */
+ txz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+ }
+ else if (src.orient == 4) { /* dipole +/O/- */
+ /* correction: subtrace previous value to prevent z-1 values. */
+ txz[iy*n1*n2+ix*n1+iz] -= 2.0*src_ampl;
+ txz[iy*n1*n2+ix*n1+iz+1] += src_ampl;
+ }
+ else if (src.orient == 5) { /* dipole + - */
+ txz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+ }
+ }
+ /* Tzz source */
+ else if(src.type == 3) {
+ tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
+ }
+ /* Txx source */
+ else if(src.type == 4) {
+ txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ }
+
+/***********************************************************************
+* pure potential shear S source (experimental)
+* Curl S-pot = CURL(F) = dF_x/dz - dF_z/dx
+***********************************************************************/
+ else if(src.type == 5) {
+ src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+ if (src.orient == 3) src_ampl = -src_ampl;
+ /* first order derivatives */
+ vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vx[(iy-1)*n1*n2+ix*n1+iz-1] -= src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vy[iy*n1*n2+(ix-1)*n1+iz-1] -= src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vz[(iy-1)*n1*n2+(ix-1)*n1+iz] += src_ampl*sdx;
+
+ /* second order derivatives */
+ /*
+ vx[ix*n1+iz] += c1*src_ampl*sdx;
+ vx[ix*n1+iz-1] -= c1*src_ampl*sdx;
+ vx[ix*n1+iz+1] += c2*src_ampl*sdx;
+ vx[ix*n1+iz-2] -= c2*src_ampl*sdx;
+
+ vz[ix*n1+iz] -= c1*src_ampl*sdx;
+ vz[(ix-1)*n1+iz] += c1*src_ampl*sdx;
+ vz[(ix+1)*n1+iz] -= c2*src_ampl*sdx;
+ vz[(ix-2)*n1+iz] += c2*src_ampl*sdx;
+ */
+
+ /* determine second position of dipole */
+ if (src.orient == 2) { /* dipole +/- vertical */
+ iz += 1;
+ vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vx[(iy-1)*n1*n2+ix*n1+iz-1] += src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vy[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
+ }
+ else if (src.orient == 3) { /* dipole - + horizontal */
+ ix += 1;
+ vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vx[(iy-1)*n1*n2+ix*n1+iz-1] += src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vy[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
+ }
+ }
+/***********************************************************************
+* pure potential pressure P source (experimental)
+* Divergence P-pot = DIV(F) = dF_x/dx + dF_z/dz
+***********************************************************************/
+ else if(src.type == 8) {
+ src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
+ if (src.orient == 3) src_ampl = -src_ampl;
+ vx[iy*n1*n2+(ix+1)*n1+iz] += src_ampl*sdx;
+ vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vy[(iy+1)*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz+1] += src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ /* determine second position of dipole */
+ if (src.orient == 2) { /* dipole +/- */
+ iz += 1;
+ vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
+ vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz+1] -= src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ }
+ else if (src.orient == 3) { /* dipole - + */
+ ix += 1;
+ vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
+ vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
+ vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz+1] -= src_ampl*sdx;
+ vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
+ }
+ }
+ 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);
+
+ txx[iy*n1*n2+ix*n1+iz] -= Mxx*src_ampl;
+ tzz[iy*n1*n2+ix*n1+iz] -= Mzz*src_ampl;
+ txz[iy*n1*n2+ix*n1+iz] -= Mxz*src_ampl;
+ } /* src.type */
+ } /* ischeme */
+ } /* loop over isrc */
+
+ return 0;
+}
diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
new file mode 100644
index 0000000..3f4f555
--- /dev/null
+++ b/fdelmodc3D/boundaries3D.c
@@ -0,0 +1,2253 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+
+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)
+{
+/*********************************************************************
+
+ AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+
+***********************************************************************/
+
+ float c1, c2;
+ float dp, dvx, dvy, dvz;
+ long ix, iy, iz, ixs, iys, izs, ibnd, ib, ibx, iby, ibz;
+ long nx, ny, nz, n1, n2, n3;
+ long is0, isrc;
+ long ixo, ixe, iyo, iye, izo, ize;
+ long npml, ipml, pml;
+ float kappu, alphu, sigmax, R, a, m, fac, dx, dy, dt;
+ float dpx, dpy, dpz, *p;
+ static float *Vxpml, *Vypml, *Vzpml, *sigmu, *RA;
+ static long allocated=0;
+ float Jx, Jy, Jz, rho, d;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ nx = mod.nx;
+ ny = mod.ny;
+ nz = mod.nz;
+ n1 = mod.naz;
+ n2 = mod.nax;
+ n3 = mod.nay;
+ dx = mod.dx;
+ dy = mod.dy;
+ dt = mod.dt;
+ fac = dt/dx;
+ if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+ else pml=0;
+
+ ibnd = mod.iorder/2-1;
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+ if (bnd.top==1) { /* free surface at top */
+#pragma omp for private (ix,iy) nowait
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ iz = bnd.surface[ix];
+ vz[iy*n2*n1+ix*n1+iz] = vz[iy*n2*n1+ix*n1+iz+1];
+ vz[iy*n2*n1+ix*n1+iz-1] = vz[iy*n2*n1+ix*n1+iz+2];
+ }
+ }
+ }
+ }
+
+/************************************************************/
+/* rigid boundary condition clears velocities on boundaries */
+/************************************************************/
+
+ if (bnd.top==3) { /* rigid surface at top */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (iy=1; iy<=ny; iy++) {
+ for (ix=1; ix<=nx; ix++) {
+ vx[iy*n2*n1+ix*n1+ibnd] = 0.0;
+ vy[iy*n2*n1+ix*n1+ibnd] = 0.0;
+ vz[iy*n2*n1+ix*n1+ibnd] = -vz[iy*n2*n1+ix*n1+ibnd+1];
+ if (mod.iorder >= 4) vz[iy*n2*n1+ix*n1+ibnd-1] = -vz[iy*n2*n1+ix*n1+ibnd+2];
+ if (mod.iorder >= 6) vz[iy*n2*n1+ix*n1+ibnd-2] = -vz[iy*n2*n1+ix*n1+ibnd+3];
+ }
+ }
+ }
+ if (bnd.rig==3) { /* rigid surface at right */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (iy=1; iy<=ny; iy++) {
+ for (iz=1; iz<=nz; iz++) {
+ vz[iy*n2*n1+(nx+ibnd-1)*n1+iz] = 0.0;
+ vy[iy*n2*n1+(nx+ibnd-1)*n1+iz] = 0.0;
+ vx[iy*n2*n1+(nx+ibnd)*n1+iz] = -vx[iy*n2*n1+(nx+ibnd-1)*n1+iz];
+ if (mod.iorder == 4) vx[iy*n2*n1+(nx+2)*n1+iz] = -vx[iy*n2*n1+(nx-1)*n1+iz];
+ if (mod.iorder == 6) {
+ vx[iy*n2*n1+(nx+1)*n1+iz] = -vx[iy*n2*n1+(nx)*n1+iz];
+ vx[iy*n2*n1+(nx+3)*n1+iz] = -vx[iy*n2*n1+(nx-2)*n1+iz];
+ }
+ }
+ }
+ }if (bnd.bac==3) { /* rigid surface at back */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (ix=1; ix<=nx; ix++) {
+ for (iz=1; iz<=nz; iz++) {
+ vz[(ny+ibnd-1)*n2*n1+ix*n1+iz] = 0.0;
+ vx[(ny+ibnd-1)*n2*n1+ix*n1+iz] = 0.0;
+ vy[(ny+ibnd)*n2*n1+ix*n1+iz] = -vy[(ny+ibnd-1)*n2*n1+ix*n1+iz];
+ if (mod.iorder == 4) vy[(ny+2)*n2*n1+ix*n1+iz] = -vy[(ny-1)*n2*n1+iy*n1+iz];
+ if (mod.iorder == 6) {
+ vy[(ny+1)*n2*n1+ix*n1+iz] = -vy[ny*n2*n1+ix*n1+iz];
+ vy[(ny+3)*n2*n1+ix*n1+iz] = -vy[(ny-2)*n2*n1+ix*n1+iz];
+ }
+ }
+ }
+ }
+ if (bnd.bot==3) { /* rigid surface at bottom */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (iy=1; iy<=ny; iy++) {
+ for (ix=1; ix<=nx; ix++) {
+ vx[iy*n2*n1+ix*n1+nz+ibnd-1] = 0.0;
+ vy[iy*n2*n1+ix*n1+nz+ibnd-1] = 0.0;
+ vz[iy*n2*n1+ix*n1+nz+ibnd] = -vz[iy*n2*n1+ix*n1+nz+ibnd-1];
+ if (mod.iorder == 4) vz[iy*n2*n1+ix*n1+nz+2] = -vz[iy*n2*n1+ix*n1+nz-1];
+ if (mod.iorder == 6) {
+ vz[iy*n2*n1+ix*n1+nz+1] = -vz[iy*n2*n1+ix*n1+nz];
+ vz[iy*n2*n1+ix*n1+nz+3] = -vz[iy*n2*n1+ix*n1+nz-2];
+ }
+ }
+ }
+ }
+ if (bnd.lef==3) { /* rigid surface at left */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (iy=1; iy<=ny; iy++) {
+ for (iz=1; iz<=nz; iz++) {
+ vz[iy*n2*n1+ibnd*n1+iz] = 0.0;
+ vy[iy*n2*n1+ibnd*n1+iz] = 0.0;
+ vx[iy*n2*n1+ibnd*n1+iz] = -vx[iy*n2*n1+(ibnd+1)*n1+iz];
+ if (mod.iorder == 4) vx[iy*n2*n1+0*n1+iz] = -vx[iy*n2*n1+3*n1+iz];
+ if (mod.iorder == 6) {
+ vx[iy*n2*n1+1*n1+iz] = -vx[iy*n2*n1+4*n1+iz];
+ vx[iy*n2*n1+0*n1+iz] = -vx[iy*n2*n1+5*n1+iz];
+ }
+ }
+ }
+ }
+ if (bnd.fro==3) { /* rigid surface at front */
+#pragma omp for private (ix, iy, iz) nowait
+#pragma ivdep
+ for (ix=1; ix<=nx; ix++) {
+ for (iz=1; iz<=nz; iz++) {
+ vz[ibnd*n2*n1+ix*n1+iz] = 0.0;
+ vx[ibnd*n2*n1+ix*n1+iz] = 0.0;
+ vy[ibnd*n2*n1+ix*n1+iz] = -vy[(ibnd+1)*n2*n1+ix*n1+iz];
+ if (mod.iorder == 4) vy[0*n2*n1+ix*n1+iz] = -vy[3*n2*n1+ix*n1+iz];
+ if (mod.iorder == 6) {
+ vy[1*n2*n1+ix*n1+iz] = -vy[4*n2*n1+ix*n1+iz];
+ vy[0*n2*n1+ix*n1+iz] = -vy[5*n2*n1+ix*n1+iz];
+ }
+ }
+ }
+ }
+
+/************************************************************/
+/* Tapered boundaries for both elastic and acoustic schemes */
+/* compute all field values in tapered areas */
+/************************************************************/
+
+ /*********/
+ /* Top */
+ /*********/
+ if (bnd.top==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ 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++) {
+ 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[ibz-iz];
+ }
+ }
+ }
+ /* right top corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieXx;
+ ixe = ixo+bnd.ntap;
+ 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++) {
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ /* right top front corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* right top back corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ioXz-bnd.ntap;
+ ize = mod.ioXz;
+
+ /* left top corner */
+ if (bnd.lef==4) {
+ ixo = mod.ioXx-bnd.ntap;
+ ixe = mod.ioXx;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ /* left top front corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* left top back corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+ /* front top corner */
+ if (bnd.fro==4) {
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* Back top corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ioXy+bnd.ntap;
+ ibz = (bnd.ntap+izo-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+(ibz-iz)];
+ }
+ }
+ }
+ }
+
+ /* Vy field */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz-bnd.ntap;
+ ize = mod.ioYz;
+
+ 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++) {
+ 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[ib-iz];
+ }
+ }
+ }
+ /* right top corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieYx;
+ ixe = ixo+bnd.ntap;
+ 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++) {
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ /* right top front corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* right top back corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz-bnd.ntap;
+ ize = mod.ioYz;
+
+ /* left top corner */
+ if (bnd.lef==4) {
+ ixo = mod.ioYx-bnd.ntap;
+ ixe = mod.ioYx;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ /* left top front corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* left top back corner */
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+ /* front top corner */
+ if (bnd.fro==4) {
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* Back top corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ioXy+bnd.ntap;
+ ibz = (bnd.ntap+izo-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+(ibz-iz)];
+ }
+ }
+ }
+ }
+
+
+ /* Vz field */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ 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[ib-iz];
+ }
+ }
+ }
+ /* right top corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieZx;
+ ixe = ixo+bnd.ntap;
+ ibz = (bnd.ntap+izo-1);
+ ibx = (ixo);
+#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+(ibz-iz)];
+ }
+ }
+ }
+ /* right top 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* right top 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ioZz-bnd.ntap;
+ ize = mod.ioZz;
+
+ /* left top corner */
+ if (bnd.lef==4) {
+ ixo = mod.ioZx-bnd.ntap;
+ ixe = mod.ioZx;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ /* left top 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+(ibx-ix)*bnd.ntap+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* left top 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+ /* front top corner */
+ if (bnd.fro==4) {
+ ixo = mod.ioXx;
+ ixe = mod.ieXx;
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ 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
+ 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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ /* Back top corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ioXy+bnd.ntap;
+ ibz = (bnd.ntap+izo-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+(ibz-iz)];
+ }
+ }
+ }
+ }
+ }
+
+ }
+
+ /*********/
+ /* Bottom */
+ /*********/
+ if (bnd.bot==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ 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 corner */
+ if (bnd.rig==4) {
+ ixo = mod.ieXx;
+ 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++) {
+ 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 */
+ 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 */
+ 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 corner */
+ 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 */
+ 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 */
+ 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)];
+ }
+ }
+ }
+ }
+ }
+ /* front bottom corner */
+ 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)];
+ }
+ }
+ }
+ }
+ /* Back bottom corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ioXy+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 */
+ ixo = mod.ioYx;
+ ixe = mod.ieYx;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz;
+ ize = mod.ioYz+bnd.ntap;
+
+ 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++) {
+ 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 corner */
+ 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 */
+ 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 */
+ 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.ioYz;
+ ize = mod.ioYz+bnd.ntap;
+
+ /* left bottom corner */
+ 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 */
+ 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 */
+ 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)];
+ }
+ }
+ }
+ }
+ }
+ /* front bottom corner */
+ 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++) {
+ 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)];
+ }
+ }
+ }
+ }
+ /* Back bottom corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ioXy+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, iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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.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,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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[(ix-ibx)*bnd.ntap+(iz-ibz)];
+ }
+ }
+ /* right bottom front corner */
+ 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++) {
+ 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.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++) {
+ 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)];
+ }
+ }
+ }
+ }
+ }
+ /* 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,iz)
+ for (ix=ixo; ix<ixe; ix++) {
+#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.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++) {
+ 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.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*n1*n2+ix*n1+iz-1]) +
+ c2*(tzz[iy*n1*n2+ix*n1+iz+1] - tzz[iy*n1*n2+ix*n1+iz-2]));
+
+ vy[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.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++) {
+ 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.ieXy;
+ iye = mod.ioXy+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)];
+ }
+ }
+ }
+ }
+
+ }
+ }
+
+ /*********/
+ /* Left */
+ /*********/
+ if (bnd.lef==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ 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 front corner */
+ 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 back corner */
+ 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 */
+ 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 front corner */
+ 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 back corner */
+ 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 */
+ 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 front corner */
+ 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 back corner */
+ 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)];
+ }
+ }
+ }
+ }
+
+ }
+
+ }
+
+ /*********/
+ /* Right */
+ /*********/
+ if (bnd.rig==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ ixo = mod.ieXx;
+ ixe = mod.ieXx+bnd.ntap;
+ iyo = mod.ioXy;
+ iye = mod.ieXy;
+ izo = mod.ioXz;
+ ize = mod.ieXz;
+
+ ib = (ixe-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[ix-ib];
+ }
+ }
+ }
+ /* left front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioXy-bnd.ntap;
+ iye = mod.ioXy;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+ /* left back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieXy;
+ iye = mod.ieXy+bnd.ntap;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+
+ /* Vy field */
+ ixo = mod.ieYx;
+ ixe = mod.ieYx+bnd.ntap;
+ iyo = mod.ioYy;
+ iye = mod.ieYy;
+ izo = mod.ioYz;
+ ize = mod.ieYz;
+
+ ib = (ixe-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.tapx[ix-ib];
+ }
+ }
+ }
+ /* left front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioYy-bnd.ntap;
+ iye = mod.ioYy;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+ /* left back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieYy;
+ iye = mod.ieYy+bnd.ntap;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+
+ /* Vz field */
+ ixo = mod.ieZx;
+ ixe = mod.ieZx+bnd.ntap;
+ iyo = mod.ioZy;
+ iye = mod.ieZy;
+ izo = mod.ioZz;
+ ize = mod.ieZz;
+
+ ib = (ixe-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[ix-ib];
+ }
+ }
+ }
+ /* left front corner */
+ if (bnd.fro==4) {
+ iyo = mod.ioZy-bnd.ntap;
+ iye = mod.ioZy;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+ /* left back corner */
+ if (bnd.bac==4) {
+ iyo = mod.ieZy;
+ iye = mod.ieZy+bnd.ntap;
+ ibx = (ixo);
+ 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+(ix-ibx)];
+ }
+ }
+ }
+ }
+
+ }
+
+ }
+
+ /*********/
+ /* Front */
+ /*********/
+ if (bnd.fro==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ 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*(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.tapy[iby-iy];
+ }
+ }
+ }
+
+ /* Vy field */
+ 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*(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[iby-iy];
+ }
+ }
+ }
+
+ /* Vz field */
+ 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*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.tapy[iby-iy];
+ }
+ }
+ }
+ }
+
+ }
+
+ /********/
+ /* Back */
+ /********/
+ if (bnd.bac==4) {
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+
+ /* Vx field */
+ 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*(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.tapy[iy-iby];
+ }
+ }
+ }
+
+ /* Vy field */
+ 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*(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[iy-iby];
+ }
+ }
+ }
+
+ /* Vz field */
+ ixo = mod.ioZx;
+ ixe = mod.ieZx;
+ iyo = mod.ioZy;
+ 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*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.tapy[iy-iby];
+ }
+ }
+ }
+ }
+
+ }
+
+ if ( (npml != 0) && (itime==mod.nt-1) && pml) {
+#pragma omp master
+{
+ if (allocated) {
+ free(Vxpml);
+ free(Vzpml);
+ free(sigmu);
+ free(RA);
+ allocated=0;
+ }
+}
+ }
+
+ return 0;
+}
+
+long boundariesV3D(modPar mod, bndPar bnd, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose)
+{
+/*********************************************************************
+
+ AUTHOR:
+ Jan Thorbecke (janth@xs4all.nl)
+ The Netherlands
+
+***********************************************************************/
+
+ float c1, c2;
+ float dp, dvx, dvy, dvz;
+ long ix, iy, iz, ixs, iys, izs, izp, ib;
+ long nx, ny, nz, n1, n2, n3;
+ long is0, isrc;
+ long ixo, ixe, iyo, iye, izo, ize;
+ long npml, ipml, ipml2, pml;
+ float kappu, alphu, sigmax, R, a, m, fac, dx, dy, dt;
+ float *p;
+ static float *Pxpml, *Pypml, *Pzpml, *sigmu, *RA;
+ static long allocated=0;
+ float Jx, Jy, Jz, rho, d;
+
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ nx = mod.nx;
+ ny = mod.ny;
+ nz = mod.nz;
+ n1 = mod.naz;
+ n2 = mod.nax;
+ n3 = mod.nay;
+ dx = mod.dx;
+ dy = mod.dy;
+ dt = mod.dt;
+ fac = dt/dx;
+ if ( (bnd.top==2) || (bnd.bot==2) || (bnd.lef==2) || (bnd.rig==2) ) pml=1;
+ else pml=0;
+
+/************************************************************/
+/* PML boundaries for acoustic schemes */
+/* compute all field values in tapered areas */
+/************************************************************/
+
+ npml=bnd.npml; /* lenght of pml in grid-points */
+ if ( (npml != 0) && (itime==0) && pml) {
+#pragma omp master
+{
+ if (allocated) {
+ free(Pxpml);
+ free(Pypml);
+ free(Pzpml);
+ free(sigmu);
+ free(RA);
+ }
+ Pxpml = (float *)calloc(2*n1*n3*npml,sizeof(float));
+ Pypml = (float *)calloc(2*n2*n1*npml,sizeof(float));
+ Pzpml = (float *)calloc(2*n3*n2*npml,sizeof(float));
+ sigmu = (float *)calloc(npml,sizeof(float));
+ RA = (float *)calloc(npml,sizeof(float));
+ allocated = 1;
+
+ /* calculate sigmu and RA only once with fixed velocity Cp */
+ m=bnd.m; /* scaling order */
+ R=bnd.R; /* the theoretical reflection coefficient after discretization */
+ kappu = 1.0; /* auxiliary attenuation coefficient for small angles */
+ alphu=0.0; /* auxiliary attenuation coefficient for low frequencies */
+ d = (npml-1)*dx; /* depth of pml */
+ /* sigmu attenuation factor representing the loss in the PML depends on the grid position in the PML */
+
+ sigmax = ((3.0*mod.cp_min)/(2.0*d))*log(1.0/R);
+ for (ib=0; ib<npml; ib++) { /* ib=0 interface between PML and interior */
+ a = (float) (ib/(npml-1.0));
+ sigmu[ib] = sigmax*pow(a,m);
+ RA[ib] = (1.0)/(1.0+0.5*dt*sigmu[ib]);
+ }
+}
+ }
+
+#pragma omp barrier
+ if (mod.ischeme == 1 && pml) { /* Acoustic scheme PML's */
+ p = tzz; /* Tzz array pointer points to P-field */
+
+ if (bnd.top==2) mod.ioPz += bnd.npml;
+ if (bnd.bot==2) mod.iePz -= bnd.npml;
+ if (bnd.lef==2) mod.ioPx += bnd.npml;
+ if (bnd.rig==2) mod.iePx -= bnd.npml;
+ if (bnd.fro==2) mod.ioPy += bnd.npml;
+ if (bnd.bac==2) mod.iePy -= bnd.npml;
+
+ /* PML top P */
+ if (bnd.top == 2) {
+ /* PML top P-Vz-component */
+#pragma omp for private (ix, iy, iz, dvx, dvy, dvz, Jz, ipml)
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ ipml2 = npml-1;
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ ipml = npml-1;
+ for (iz=mod.ioPz-npml; iz<mod.ioPz; 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]);
+ Jz = RA[ipml2]*RA[ipml]*dvz - RA[ipml2]*RA[ipml]*dt*Pzpml[ix*npml+ipml];
+ Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz+dvx);
+ ipml--;
+ }
+ }
+ }
+ }
+
+ /* PML left P */
+ if (bnd.lef == 2) {
+ /* PML left P-Vx-component */
+#pragma omp for private (ix, iz, dvx, dvz, Jx, ipml)
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ ipml = npml-1;
+ for (ix=mod.ioPx-npml; ix<mod.ioPx; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml];
+ Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx+dvz);
+ ipml--;
+ }
+ }
+ }
+
+ /* PML corner left-top P */
+ if (bnd.lef == 2 && bnd.top == 2) {
+ /* PML left P-Vx-component */
+#pragma omp for private (ix, iz, dvx, Jx, ipml)
+ for (iz=mod.ioPz-npml; iz<mod.ioPz; iz++) {
+ ipml = npml-1;
+ for (ix=mod.ioPx-npml; ix<mod.ioPx; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml];
+ Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx);
+ ipml--;
+ }
+ }
+ /* PML top P-Vz-component */
+#pragma omp for private (ix, iz, dvz, Jz, ipml)
+ for (ix=mod.ioPx-npml; ix<mod.ioPx; ix++) {
+ ipml = npml-1;
+ for (iz=mod.ioPz-npml; iz<mod.ioPz; iz++) {
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[ix*npml+ipml];
+ Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz);
+ ipml--;
+ }
+ }
+ }
+
+ /* PML right P */
+ if (bnd.rig == 2) {
+ /* PML right P Vx-component */
+#pragma omp for private (ix, iz, dvx, dvz, Jx, ipml)
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ ipml = 0;
+ for (ix=mod.iePx; ix<mod.iePx+npml; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml];
+ Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx+dvz);
+ ipml++;
+ }
+ }
+ }
+
+ /* PML corner right-top P */
+ if (bnd.rig == 2 && bnd.top == 2) {
+ /* PML right P Vx-component */
+#pragma omp for private (ix, iz, dvx, Jx, ipml)
+ for (iz=mod.ioPz-npml; iz<mod.ioPz; iz++) {
+ ipml = 0;
+ for (ix=mod.iePx; ix<mod.iePx+npml; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml];
+ Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx);
+ ipml++;
+ }
+ }
+ /* PML top P-Vz-component */
+#pragma omp for private (ix, iz, dvz, Jz, ipml)
+ for (ix=mod.iePx; ix<mod.iePx+npml; ix++) {
+ ipml = npml-1;
+ for (iz=mod.ioPz-npml; iz<mod.ioPz; iz++) {
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[ix*npml+ipml];
+ Pzpml[ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz);
+ ipml--;
+ }
+ }
+ }
+
+ /* PML bottom P */
+ if (bnd.bot == 2) {
+ /* PML bottom P Vz-component */
+#pragma omp for private (ix, iz, dvx, dvz, Jz, ipml)
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ ipml = 0;
+ for (iz=mod.iePz; iz<mod.iePz+npml; iz++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml];
+ Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz+dvx);
+ ipml++;
+ }
+ }
+ }
+
+ /* PML corner bottom-right P */
+ if (bnd.bot == 2 && bnd.rig == 2) {
+ /* PML bottom P Vz-component */
+#pragma omp for private (ix, iz, dvz, Jz, ipml)
+ for (ix=mod.iePx; ix<mod.iePx+npml; ix++) {
+ ipml = 0;
+ for (iz=mod.iePz; iz<mod.iePz+npml; iz++) {
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml];
+ Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz);
+ ipml++;
+ }
+ }
+ /* PML right P Vx-component */
+#pragma omp for private (ix, iz, dvx, Jx, ipml)
+ for (iz=mod.iePz; iz<mod.iePz+npml; iz++) {
+ ipml = 0;
+ for (ix=mod.iePx; ix<mod.iePx+npml; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[n1*npml+iz*npml+ipml];
+ Pxpml[n1*npml+iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx);
+ //p[ix*n1+iz] -= l2m[ix*n1+iz]*(dvx);
+ ipml++;
+ }
+ }
+ }
+
+ /* PML corner left-bottom P */
+ if (bnd.bot == 2 && bnd.lef == 2) {
+ /* PML bottom P Vz-component */
+#pragma omp for private (ix, iz, dvz, Jz, ipml)
+ for (ix=mod.ioPx-npml; ix<mod.ioPx; ix++) {
+ ipml = 0;
+ for (iz=mod.iePz; iz<mod.iePz+npml; iz++) {
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ Jz = RA[ipml]*dvz - RA[ipml]*dt*Pzpml[n2*npml+ix*npml+ipml];
+ Pzpml[n2*npml+ix*npml+ipml] += sigmu[ipml]*Jz;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jz);
+ ipml++;
+ }
+ }
+ /* PML left P Vx-component */
+#pragma omp for private (ix, iz, dvx, Jx, ipml)
+ for (iz=mod.iePz; iz<mod.iePz+npml; iz++) {
+ ipml = npml-1;
+ for (ix=mod.ioPx-npml; ix<mod.ioPx; ix++) {
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ Jx = RA[ipml]*dvx - RA[ipml]*dt*Pxpml[iz*npml+ipml];
+ Pxpml[iz*npml+ipml] += sigmu[ipml]*Jx;
+ p[ix*n1+iz] -= l2m[ix*n1+iz]*(Jx);
+ ipml--;
+ }
+ }
+ }
+ if (bnd.top==2) mod.ioPz -= bnd.npml;
+ if (bnd.bot==2) mod.iePz += bnd.npml;
+ if (bnd.lef==2) mod.ioPx -= bnd.npml;
+ if (bnd.rig==2) mod.iePx += bnd.npml;
+
+ } /* end acoustic PML */
+
+
+
+/****************************************************************/
+/* Free surface: calculate free surface conditions for stresses */
+/****************************************************************/
+
+
+ ixo = mod.ioPx;
+ ixe = mod.iePx;
+ iyo = mod.ioPy;
+ iye = mod.iePy;
+ izo = mod.ioPz;
+ ize = mod.iePz;
+
+ if (mod.ischeme <= 2) { /* Acoustic scheme */
+ if (bnd.top==1) { /* free surface at top */
+#pragma omp for private (ix,iy) nowait
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ iz = bnd.surface[iy*n2+ix];
+ tzz[iy*n2*n1+ix*n1+iz] = 0.0;
+ //vz[ix*n1+iz] = -vz[ix*n1+iz+1];
+ //vz[ix*n1+iz-1] = -vz[ix*n1+iz+2];
+ }
+ }
+ }
+ if (bnd.rig==1) { /* free surface at right */
+#pragma omp for private (iy,iz) nowait
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ tzz[iy*n1*n2+(mod.iePx-1)*n1+iz] = 0.0;
+ }
+ }
+ }
+ if (bnd.fro==1) { /* free surface at front */
+#pragma omp for private (ix,iz) nowait
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ tzz[(mod.ioPy-1)*n1*n2+ix*n1+iz] = 0.0;
+ }
+ }
+ }
+ if (bnd.bot==1) { /* free surface at bottom */
+#pragma omp for private (ix,iy) nowait
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ tzz[iy*n1*n2+ix*n1+mod.iePz-1] = 0.0;
+ }
+ }
+ }
+ if (bnd.lef==1) { /* free surface at left */
+#pragma omp for private (iy,iz) nowait
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ tzz[iy*n1*n2+(mod.ioPx-1)*n1+iz] = 0.0;
+ }
+ }
+ }
+ if (bnd.bac==1) { /* free surface at back */
+#pragma omp for private (ix,iz) nowait
+ for (iz=mod.ioPz; iz<mod.iePz; iz++) {
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ tzz[(mod.iePy-1)*n1*n2+ix*n1+iz] = 0.0;
+ }
+ }
+ }
+ }
+
+ if ( (npml != 0) && (itime==mod.nt-1) && pml) {
+#pragma omp master
+{
+ if (allocated) {
+ free(Pxpml);
+ free(Pypml);
+ free(Pzpml);
+ free(sigmu);
+ free(RA);
+ allocated=0;
+ }
+}
+ }
+
+ return 0;
+}
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index a79ede1..8ba30a7 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -589,51 +589,41 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
// vx, vz, tzz, rox, roz, l2m, verbose);
// break;
case -1 : /* Acoustic dissipative media FD kernel */
- acoustic4_qr(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, verbose);
+ vmess("Acoustic dissipative not yet available");
break;
case 1 : /* Acoustic FD kernel */
if (mod.iorder==2) {
- acoustic2(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, verbose);
+ vmess("Acoustic order 2 not yet available");
}
else if (mod.iorder==4) {
if (mod.sh) {
- acousticSH4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, verbose);
+ vmess("SH order 4 not yet available");
}
else {
- acoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, verbose);
+ acoustic4_3D(mod, src, wav, bnd, it, ixsrc, iysrc, izsrc, src_nwav,
+ vx, vy, vz, tzz, rox, roy, roz, l2m, verbose);
}
}
else if (mod.iorder==6) {
- acoustic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, verbose);
+ vmess("Acoustic order 6 not yet available");
}
break;
case 2 : /* Visco-Acoustic FD kernel */
- viscoacoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, rox, roz, l2m, tss, tep, q, verbose);
+ vmess("Visco-Acoustic order 4 not yet available");
break;
case 3 : /* Elastic FD kernel */
if (mod.iorder==4) {
- elastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+ vmess("Elastic order 4 not yet available");
}
else if (mod.iorder==6) {
- elastic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+ vmess("Elastic order 6 not yet available");
}
break;
case 4 : /* Visco-Elastic FD kernel */
- viscoelastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul,
- tss, tep, tes, r, q, p, verbose);
+ vmess("Visco-Elastic order 4 not yet available");
break;
case 5 : /* Elastic FD kernel with S-velocity set to zero*/
- elastic4dc(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
- vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
+ vmess("DC-Elastic order 4 not yet available");
break;
}
diff --git a/fdelmodc3D/getRecTimes3D.c b/fdelmodc3D/getRecTimes3D.c
new file mode 100644
index 0000000..ea2c276
--- /dev/null
+++ b/fdelmodc3D/getRecTimes3D.c
@@ -0,0 +1,324 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"fdelmodc3D.h"
+#include"par.h"
+
+/**
+* Stores the wavefield at the receiver positions.
+*
+* On a staggered grid the fields are all on different positions,
+* to compensate for that the rec.int_vx and rec.int_vz options
+* can be set.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+int getRecTimes3D(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *l2m, float *rox, float *roy, float *roz, float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_txy, float *rec_tyz, float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, int verbose)
+{
+ int n1, n2, ibndx, ibndy, ibndz;
+ int irec, ix, iy, iz, ix2, iy2, iz2, ix1, iy1, iz1;
+ float dvx, dvy, dvz, rdz, rdy, rdx;
+ float C000, C100, C010, C001, C110, C101, C011, C111;
+ float *vz_t, c1, c2, lroz, field;
+
+ ibndx = mod.ioPx;
+ ibndy = mod.ioPy;
+ ibndz = mod.ioPz;
+ if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
+ if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
+ if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
+ n1 = mod.naz;
+ n2 = mod.nax;
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+
+ if (!rec.n) return 0;
+
+/***********************************************************************
+* velocity or txz or potential registrations issues:
+* rec_x and rec_z are related to actual txx/tzz/p positions.
+* offsets from virtual boundaries must be taken into account.
+*
+* vx velocities have one sample less in x-direction
+* vz velocities have one sample less in z-direction
+* txz stresses have one sample less in z-direction and x-direction
+*
+* Note, in the acoustic scheme P is stored in the Tzz array.
+***********************************************************************/
+
+ for (irec=0; irec<rec.n; irec++) {
+ iz = rec.z[irec]+ibndz;
+ iy = rec.y[irec]+ibndy;
+ ix = rec.x[irec]+ibndx;
+ iz1 = iz-1;
+ iy1 = iy-1;
+ ix1 = ix-1;
+ iz2 = iz+1;
+ iy2 = iy+1;
+ ix2 = ix+1;
+ /* interpolation to precise (not necessary on a grid point) position */
+ if ( rec.int_p==3 ) {
+
+ iz = (int)floorf(rec.zr[irec]/mod.dz)+ibndz;
+ iy = (int)floorf(rec.yr[irec]/mod.dy)+ibndy;
+ ix = (int)floorf(rec.xr[irec]/mod.dx)+ibndx;
+ rdz = (rec.zr[irec] - (iz-ibndz)*mod.dz)/mod.dz;
+ rdy = (rec.yr[irec] - (iy-ibndy)*mod.dy)/mod.dy;
+ rdx = (rec.xr[irec] - (ix-ibndx)*mod.dx)/mod.dx;
+ iz1 = iz-1;
+ iy1 = iy-1;
+ ix1 = ix-1;
+ iz2 = iz+1;
+ iy2 = iy+1;
+ ix2 = ix+1;
+
+ /*
+ // Interpolate according to Dirk Kraaijpool's scheme
+ // Reference: "Seismic ray fields and ray field maps : theory and algorithms" ,
+ // PhD thesis Utrecht University,Faculty of Geosciences, 2003)
+
+ C00 = tzz[ix*n1+iz] + 0.5*((tzz[(ix+1)*n1+iz] +tzz[(ix-1)*n1+iz]+
+ tzz[(ix )*n1+iz+1] +tzz[(ix )*n1+iz-1])/(2.0*mod.dx));
+ C10 = tzz[(ix+1)*n1+iz] + 0.5*((tzz[(ix+2)*n1+iz] +tzz[(ix )*n1+iz]+
+ tzz[(ix+1)*n1+iz+1] +tzz[(ix+1)*n1+iz-1])/(2.0*mod.dz));
+ C01 = tzz[ix*n1+iz+1] + 0.5*((tzz[(ix+1)*n1+iz+1] +tzz[(ix-1)*n1+iz+1]+
+ tzz[(ix)*n1+iz+2] +tzz[(ix )*n1+iz])/(2.0*mod.dx));
+ C11 = tzz[(ix+1)*n1+iz+1]+ 0.5*((tzz[(ix+2)*n1+iz+1] +tzz[(ix )*n1+iz+1]+
+ tzz[(ix+1)*n1+iz+2] +tzz[(ix+1)*n1+iz])/(2.0*mod.dz));
+ */
+
+ if (rec.type.p){
+ /* bi-linear interpolation */
+ C000 = tzz[iy*n1*n2+ix*n1+iz];
+ C100 = tzz[iy*n1*n2+(ix+1)*n1+iz];
+ C010 = tzz[iy*n1*n2+ix*n1+iz+1];
+ C001 = tzz[(iy+1)*n1*n2+ix*n1+iz];
+ C110 = tzz[iy*n1*n2+(ix+1)*n1+iz+1];
+ C101 = tzz[(iy+1)*n1*n2+ix*n1+iz+1];
+ C011 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz];
+ C111 = tzz[(iy+1)*n1*n2+(ix+1)*n1+iz+1];
+ rec_p[irec*rec.nt+isam] = C000*(1.0-rdx)*(1.0-rdy)*(1.0-rdz) +
+ C100*rdx*(1.0-rdy)*(1.0-rdz) +
+ C010*(1.0-rdx)*(1.0-rdy)*rdz +
+ C11*rdx*rdz;
+ }
+ if (rec.type.txx) {
+ C00 = txx[ix*n1+iz];
+ C10 = txx[(ix+1)*n1+iz];
+ C01 = txx[ix*n1+iz+1];
+ C11 = txx[(ix+1)*n1+iz+1];
+ rec_txx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.tzz) {
+ C00 = tzz[ix*n1+iz];
+ C10 = tzz[(ix+1)*n1+iz];
+ C01 = tzz[ix*n1+iz+1];
+ C11 = tzz[(ix+1)*n1+iz+1];
+ rec_tzz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.txz) {
+ C00 = txz[ix2*n1+iz2];
+ C10 = txz[(ix2+1)*n1+iz2];
+ C01 = txz[ix2*n1+iz2+1];
+ C11 = txz[(ix2+1)*n1+iz2+1];
+ rec_txz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.pp) {
+ C00 = (vx[ix2*n1+iz]-vx[ix*n1+iz] +
+ vz[ix*n1+iz2]-vz[ix*n1+iz])/mod.dx;
+ C10 = (vx[(ix2+1)*n1+iz]-vx[(ix+1)*n1+iz] +
+ vz[(ix+1)*n1+iz2]-vz[(ix+1)*n1+iz])/mod.dx;
+ C01 = (vx[ix2*n1+iz+1]-vx[ix*n1+iz+1] +
+ vz[ix*n1+iz2+1]-vz[ix*n1+iz+1])/mod.dx;
+ C11 = (vx[(ix2+1)*n1+iz+1]-vx[(ix+1)*n1+iz+1] +
+ vz[(ix+1)*n1+iz2+1]-vz[(ix+1)*n1+iz+1])/mod.dx;
+ rec_pp[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.ss) {
+ C00 = (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
+ (vz[ix2*n1+iz2]-vz[ix*n1+iz2]))/mod.dx;
+ C10 = (vx[(ix2+1)*n1+iz2]-vx[(ix2+1)*n1+iz] -
+ (vz[(ix2+1)*n1+iz2]-vz[(ix+1)*n1+iz2]))/mod.dx;
+ C01 = (vx[ix2*n1+iz2+1]-vx[ix2*n1+iz+1] -
+ (vz[ix2*n1+iz2+1]-vz[ix*n1+iz2+1]))/mod.dx;;
+ C11 = (vx[(ix2+1)*n1+iz2+1]-vx[(ix2+1)*n1+iz+1] -
+ (vz[(ix2+1)*n1+iz2+1]-vz[(ix+1)*n1+iz2+1]))/mod.dx;
+ rec_ss[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.vz) {
+ C00 = vz[ix*n1+iz2];
+ C10 = vz[(ix+1)*n1+iz2];
+ C01 = vz[ix*n1+iz2+1];
+ C11 = vz[(ix+1)*n1+iz2+1];
+ rec_vz[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+ if (rec.type.vx) {
+ C00 = vx[ix2*n1+iz];
+ C10 = vx[(ix2+1)*n1+iz];
+ C01 = vx[ix2*n1+iz+1];
+ C11 = vx[(ix2+1)*n1+iz+1];
+ rec_vx[irec*rec.nt+isam] = C00*(1.0-rdx)*(1.0-rdz) + C10*rdx*(1.0-rdz) +
+ C01*(1.0-rdx)*rdz + C11*rdx*rdz;
+ }
+
+ }
+ else { /* read values directly from the grid points */
+ if (verbose>=4 && isam==0) {
+ vmess("Receiver %d read at gridpoint ix=%d iz=%d",irec, ix, iz);
+ }
+ /* interpolation of receivers to same time step is only done for acoustic scheme */
+ if (rec.type.p) {
+ if (rec.int_p == 1) {
+ if (mod.ischeme == 1) { /* interpolate Tzz times -1/2 Dt backward to Vz times */
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ field = tzz[ix*n1+iz] + 0.5*l2m[ix*n1+iz]*(dvx+dvz);
+ dvx = c1*(vx[(ix+1)*n1+iz1] - vx[ix*n1+iz1]) +
+ c2*(vx[(ix+2)*n1+iz1] - vx[(ix-1)*n1+iz1]);
+ dvz = c1*(vz[ix*n1+iz1+1] - vz[ix*n1+iz1]) +
+ c2*(vz[ix*n1+iz1+2] - vz[ix*n1+iz1-1]);
+ field += tzz[ix*n1+iz1] + 0.5*l2m[ix*n1+iz1]*(dvx+dvz);
+ rec_p[irec*rec.nt+isam] = 0.5*field;
+ }
+ else {
+ rec_p[irec*rec.nt+isam] = 0.5*(tzz[ix*n1+iz1]+tzz[ix*n1+iz]);
+ }
+ }
+ else if (rec.int_p == 2) {
+ if (mod.ischeme == 1) { /* interpolate Tzz times -1/2 Dt backward to Vx times */
+ dvx = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
+ dvz = c1*(vz[ix*n1+iz+1] - vz[ix*n1+iz]) +
+ c2*(vz[ix*n1+iz+2] - vz[ix*n1+iz-1]);
+ field = tzz[ix*n1+iz] + 0.5*l2m[ix*n1+iz]*(dvx+dvz);
+ dvx = c1*(vx[(ix1+1)*n1+iz] - vx[ix1*n1+iz]) +
+ c2*(vx[(ix+2)*n1+iz] - vx[(ix1-1)*n1+iz]);
+ dvz = c1*(vz[ix1*n1+iz+1] - vz[ix1*n1+iz]) +
+ c2*(vz[ix1*n1+iz+2] - vz[ix1*n1+iz-1]);
+ field += tzz[ix1*n1+iz] + 0.5*l2m[ix1*n1+iz]*(dvx+dvz);
+ rec_p[irec*rec.nt+isam] = 0.5*field;
+ }
+ else {
+ rec_p[irec*rec.nt+isam] = 0.5*(tzz[ix1*n1+iz]+tzz[ix*n1+iz]);
+ }
+ }
+ else {
+ rec_p[irec*rec.nt+isam] = tzz[ix*n1+iz];
+ }
+ }
+ if (rec.type.txx) rec_txx[irec*rec.nt+isam] = txx[ix*n1+iz];
+ if (rec.type.tzz) rec_tzz[irec*rec.nt+isam] = tzz[ix*n1+iz];
+ if (rec.type.txz) { /* time interpolation to be done */
+ if (rec.int_vz == 2 || rec.int_vx == 2) {
+ rec_txz[irec*rec.nt+isam] = 0.25*(
+ txz[ix*n1+iz2]+txz[ix2*n1+iz2]+
+ txz[ix*n1+iz]+txz[ix2*n1+iz]);
+ }
+ else {
+ rec_txz[irec*rec.nt+isam] = txz[ix2*n1+iz2];
+ }
+ }
+ if (rec.type.pp) {
+ rec_pp[irec*rec.nt+isam] = (vx[ix2*n1+iz]-vx[ix*n1+iz] +
+ vz[ix*n1+iz2]-vz[ix*n1+iz])/mod.dx;
+ }
+ if (rec.type.ss) {
+ rec_ss[irec*rec.nt+isam] = (vx[ix2*n1+iz2]-vx[ix2*n1+iz] -
+ (vz[ix2*n1+iz2]-vz[ix*n1+iz2]))/mod.dx;
+ }
+ if (rec.type.vz) {
+/* interpolate vz to vx position to the right and above of vz */
+ if (rec.int_vz == 1) {
+ rec_vz[irec*rec.nt+isam] = 0.25*(
+ vz[ix*n1+iz2]+vz[ix1*n1+iz2]+
+ vz[ix*n1+iz] +vz[ix1*n1+iz]);
+ }
+/* interpolate vz to Txx/Tzz position by taking the mean of 2 values */
+ else if (rec.int_vz == 2) {
+ if (mod.ischeme == 1) { /* interpolate Vz times +1/2 Dt forward to P times */
+ field = vz[ix*n1+iz] - 0.5*roz[ix*n1+iz]*(
+ c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1]) +
+ c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2]));
+ field += vz[ix*n1+iz2] - 0.5*roz[ix*n1+iz2]*(
+ c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
+ c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
+ rec_vz[irec*rec.nt+isam] = 0.5*field;
+ }
+ else {
+ rec_vz[irec*rec.nt+isam] = 0.5*(vz[ix*n1+iz2]+vz[ix*n1+iz]);
+ }
+ }
+ else {
+ rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz2];
+ //rec_vz[irec*rec.nt+isam] = vz[ix*n1+iz];
+ //fprintf(stderr,"isam=%d vz[%d]=%e vz[%d]=%e vz[%d]=%e \n",isam, iz-1,vz[ix*n1+iz-1],iz,vz[ix*n1+iz], iz+1, vz[ix*n1+iz+1]);
+ }
+ }
+ if (rec.type.vx) {
+/* interpolate vx to vz position to the left and below of vx */
+ if (rec.int_vx == 1) {
+ rec_vx[irec*rec.nt+isam] = 0.25*(
+ vx[ix2*n1+iz]+vx[ix2*n1+iz1]+
+ vx[ix*n1+iz]+vx[ix*n1+iz1]);
+ }
+/* interpolate vx to Txx/Tzz position by taking the mean of 2 values */
+ else if (rec.int_vx == 2) {
+ if (mod.ischeme == 1) { /* interpolate Vx times +1/2 Dt forward to P times */
+ field = vx[ix*n1+iz] - 0.5*rox[ix*n1+iz]*(
+ c1*(tzz[ix*n1+iz] - tzz[(ix-1)*n1+iz]) +
+ c2*(tzz[(ix+1)*n1+iz] - tzz[(ix-2)*n1+iz]));
+ field += vx[ix2*n1+iz] - 0.5*rox[ix2*n1+iz]*(
+ c1*(tzz[ix2*n1+iz] - tzz[(ix2-1)*n1+iz]) +
+ c2*(tzz[(ix2+1)*n1+iz] - tzz[(ix2-2)*n1+iz]));
+ rec_vx[irec*rec.nt+isam] = 0.5*field;
+ }
+ else {
+ rec_vx[irec*rec.nt+isam] = 0.5*(vx[ix2*n1+iz]+vx[ix*n1+iz]);
+ }
+ }
+ else {
+ rec_vx[irec*rec.nt+isam] = vx[ix2*n1+iz];
+ }
+ }
+ }
+
+ } /* end of irec loop */
+
+ /* store all x-values on z-level for P Vz for up-down decomposition */
+ if (rec.type.ud) {
+ iz = rec.z[0]+ibndz;
+ iz2 = iz+1;
+ vz_t = (float *)calloc(2*mod.nax,sizeof(float));
+ /* P and Vz are staggered in time and need to correct for this */
+ /* -1- compute Vz at next time step and average with current time step */
+ lroz = mod.dt/(mod.dx*rec.rho);
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+ vz_t[ix] = vz[ix*n1+iz] - lroz*(
+ c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1]) +
+ c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2]));
+ vz_t[mod.nax+ix] = vz[ix*n1+iz2] - lroz*(
+ c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
+ c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
+ }
+ for (ix=0; ix<mod.nax; ix++) {
+ /* -2- compute average in time and depth to get Vz at same depth and time as P */
+ rec_udvz[ix*rec.nt+isam] = 0.25*(vz[ix*n1+iz2]+vz[ix*n1+iz]+vz_t[mod.nax+ix]+vz_t[ix]);
+ rec_udp[ix*rec.nt+isam] = tzz[ix*n1+iz];
+ }
+ free(vz_t);
+ }
+
+ return 0;
+}
diff --git a/fdelmodc3D/writeSrcRecPos3D.c b/fdelmodc3D/writeSrcRecPos3D.c
new file mode 100644
index 0000000..92b20b6
--- /dev/null
+++ b/fdelmodc3D/writeSrcRecPos3D.c
@@ -0,0 +1,160 @@
+#include<stdlib.h>
+#include<stdio.h>
+#include<math.h>
+#include<assert.h>
+#include"par.h"
+#include"fdelmodc3D.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* Writes the source and receiver positions into a gridded file,
+* which has the same size as the input gridded model files.
+* Source positions have a value +1 and receivers -1.
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2);
+
+long writeSrcRecPos3D(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
+{
+ FILE *fp;
+ float *dum, sub_x0, sub_y0, sub_z0, dx, dy, dz;
+ long is, nx, ny, nz, is0, ish, ix, iy, iz, ndot, idx, idy, idz;
+ char tmpname[1024];
+
+ ndot = 2;
+ nx = mod->nx;
+ ny = mod->ny;
+ nz = mod->nz;
+ dx = mod->dx;
+ dy = mod->dy;
+ dz = mod->dz;
+ sub_x0 = mod->x0;
+ sub_y0 = mod->y0;
+ sub_z0 = mod->z0;
+
+ /* write velocity field with positions of the sources */
+ dum = (float *)calloc(nx*ny*nz, sizeof(float));
+ vmess("Positions: shot=%li src=%li rec=%li", shot->n, src->n, rec->n);
+ /* source positions for random shots */
+ if (src->random) {
+ sprintf(tmpname,"SrcPositions%li.txt",src->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<src->n; is++) {
+ for (idy=0; idy<=ndot; idy++) {
+ for (idx=0; idx<=ndot; idx++) {
+ for (idz=0; idz<=ndot; idz++) {
+ dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
+ }
+ }
+ }
+ fprintf(fp, "%f %f %f\n", src->z[is]*dz+sub_z0, src->y[is]*dy+sub_y0, src->x[is]*dx+sub_x0);
+ }
+ fclose(fp);
+ }
+ /* source positions for single shot sources with plane waves */
+ else if (src->plane) {
+ is0 = -1*floor((src->n-1)/2);
+ sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+ fp = fopen(tmpname, "w+");
+ for (ish=0; ish<shot->n; ish++) {
+ for (is=0; is<src->n; is++) {
+ ix = shot->x[ish] + 1 + is0 + is;
+ iy = shot->y[ish] + 1 + is0 + is;
+ iz = shot->z[ish] + 1;
+ dum[iy*nx*nz+ix*nz+iz] = 1.0;
+
+ dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
+ dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
+ dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
+ dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+ dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
+ dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
+ fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
+ }
+ }
+ fclose(fp);
+ }
+ else if (src->multiwav) {
+ /* source positions for single shot sources with multiple wavelets */
+ sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+ fp = fopen(tmpname, "w+");
+ for (ish=0; ish<shot->n; ish++) {
+ for (is=0; is<src->n; is++) {
+ ix = src->x[is];
+ iy = src->y[is];
+ iz = src->z[is];
+ dum[iy*nx*nz+ix*nz+iz] = 1.0;
+
+ dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
+ dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
+ dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
+ dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
+ dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
+ dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
+ fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
+ }
+ }
+ fclose(fp);
+ }
+ else {
+ sprintf(tmpname,"SrcPositions%li.txt",shot->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<shot->n; is++) {
+ for (idy=0; idy<=ndot; idy++) {
+ for (idx=0; idx<=ndot; idx++) {
+ for (idz=0; idz<=ndot; idz++) {
+ dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
+ dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
+ }
+ }
+ }
+ fprintf(fp, "%f %f %f\n", shot->z[is]*dz+sub_z0, shot->y[is]*dy+sub_y0, shot->x[is]*dx+sub_x0);
+ }
+ fclose(fp);
+ }
+
+ /* receiver positions */
+ sprintf(tmpname,"RcvPositions%li.txt",rec->n);
+ fp = fopen(tmpname, "w+");
+ for (is=0; is<rec->n; is++) {
+ dum[rec->y[is]*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+ dum[(MAX(0,rec->y[is]-1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+ dum[(MIN(ny-1,rec->y[is]+1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
+ dum[rec->y[is]*nx*nz+(MAX(0,rec->x[is]-1))*nz+rec->z[is]] = -1.0;
+ dum[rec->y[is]*nx*nz+(MIN(nx-1,rec->x[is]+1))*nz+rec->z[is]] = -1.0;
+ dum[rec->y[is]*nx*nz+rec->x[is]*nz+MAX(0,rec->z[is]-1)] = -1.0;
+ dum[rec->y[is]*nx*nz+rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
+
+ if (rec->int_vx==3) {
+ fprintf(fp, "(%f, %f, %f)\n", rec->xr[is]*dx+sub_x0, rec->yr[is]*dy+sub_y0, rec->zr[is]*dz+sub_z0);
+ }
+ else {
+ fprintf(fp, "(%f, %f, %f)\n", rec->x[is]*dx+sub_x0, rec->y[is]*dy+sub_y0, rec->z[is]*dz+sub_z0);
+ }
+ }
+ fclose(fp);
+ writesufile3D("SrcRecPositions.su", dum, nz, nx*ny, sub_z0, sub_x0, dz, dx);
+ free(dum);
+
+ return 0;
+}
diff --git a/fdelmodc3D/writesufile3D.c b/fdelmodc3D/writesufile3D.c
new file mode 100644
index 0000000..6f2af26
--- /dev/null
+++ b/fdelmodc3D/writesufile3D.c
@@ -0,0 +1,162 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <string.h>
+#include "par.h"
+#include "fdelmodc3D.h"
+#include "SUsegy.h"
+#include "segy.h"
+
+/**
+* Writes an 2D array to a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+#define TRCBYTES 240
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define ISODD(n) ((n) & 01)
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+ FILE *file_out;
+ size_t nwrite, itrace;
+ long ns;
+ segy *hdr;
+
+/* Read in parameters */
+
+
+
+ if (n1 > USHRT_MAX) {
+ vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
+ }
+ ns = MIN(n1,USHRT_MAX);
+
+ file_out = fopen( filename, "w+" );
+ assert( file_out );
+
+ hdr = (segy *)calloc(1,TRCBYTES);
+ hdr->ns = ns;
+ hdr->dt = NINT(1000000*(d1));
+ hdr->d1 = d1;
+ hdr->d2 = d2;
+ hdr->f1 = f1;
+ hdr->f2 = f2;
+ hdr->fldr = 1;
+ hdr->trwf = n2;
+
+ for (itrace=0; itrace<n2; itrace++) {
+ hdr->tracl = itrace+1;
+ nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
+ assert (nwrite == TRCBYTES);
+ nwrite = fwrite( &data[itrace*n1], sizeof(float), ns, file_out );
+ assert (nwrite == ns);
+ }
+ fclose(file_out);
+ free(hdr);
+
+ return 0;
+}
+
+/**
+* Writes an 2D array to a SU file
+* special routine for src_nwav array which has a different number of samples for each shot
+*
+**/
+
+long writesufilesrcnwav3D(char *filename, float **src_nwav, wavPar wav, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+ FILE *file_out;
+ size_t nwrite, itrace;
+ float *trace;
+ long ns;
+ segy *hdr;
+
+/* Read in parameters */
+
+ if (n1 > USHRT_MAX) {
+ vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
+ }
+ ns = MIN(n1,USHRT_MAX);
+
+ file_out = fopen( filename, "w+" );
+ assert( file_out );
+
+ trace = (float *)malloc(n1*sizeof(float));
+ hdr = (segy *)calloc(1,TRCBYTES);
+ hdr->ns = ns;
+ hdr->dt = NINT(1000000*(d1));
+ hdr->d1 = d1;
+ hdr->d2 = d2;
+ hdr->f1 = f1;
+ hdr->f2 = f2;
+ hdr->fldr = 1;
+ hdr->trwf = n2;
+
+ for (itrace=0; itrace<n2; itrace++) {
+ hdr->tracl = itrace+1;
+ nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
+ assert (nwrite == TRCBYTES);
+ memset(trace, 0, n1*sizeof(float));
+ memcpy(trace, &src_nwav[itrace][0], wav.nsamp[itrace]*sizeof(float));
+ nwrite = fwrite( &trace[0], sizeof(float), ns, file_out );
+ assert (nwrite == ns);
+ }
+ fclose(file_out);
+ free(hdr);
+ free(trace);
+
+ return 0;
+}
+
+/**
+* Writes an 2D array to a SU file
+* special routine which used segyhdrs which have ns defined as integer (32 bit)
+* to handle more than 2^16 samples per trace.
+*
+**/
+
+long writeSUfile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
+{
+ FILE *file_out;
+ size_t nwrite, itrace;
+ SUsegy *SUhdr;
+ char *ptr;
+
+/* Read in parameters */
+
+ ptr = strstr(filename, " ");
+ *ptr = '\0';
+
+ file_out = fopen( filename, "w+" );
+ assert( file_out );
+
+ SUhdr = (SUsegy *)calloc(1,TRCBYTES);
+ SUhdr->ns = n1;
+ SUhdr->dt = NINT(1000000*(d1));
+ SUhdr->d1 = d1;
+ SUhdr->d2 = d2;
+ SUhdr->f1 = f1;
+ SUhdr->f2 = f2;
+ SUhdr->fldr = 1;
+ SUhdr->trwf = n2;
+
+ for (itrace=0; itrace<n2; itrace++) {
+ SUhdr->tracl = itrace+1;
+ nwrite = fwrite( SUhdr, 1, TRCBYTES, file_out );
+ assert (nwrite == TRCBYTES);
+ nwrite = fwrite( &data[itrace*n1], sizeof(float), n1, file_out );
+ assert (nwrite == n1);
+ }
+ fclose(file_out);
+ free(SUhdr);
+
+ return 0;
+}
+
--
GitLab