diff --git a/AMSS_NCKU_source/fdderivs_c.C b/AMSS_NCKU_source/fdderivs_c.C index 4ae31d4..fc0694f 100644 --- a/AMSS_NCKU_source/fdderivs_c.C +++ b/AMSS_NCKU_source/fdderivs_c.C @@ -71,131 +71,106 @@ void fdderivs(const int ex[3], const double Fdxdz = F1o144 / (dX * dZ); const double Fdydz = F1o144 / (dY * dZ); - /* 只清零不被主循环覆盖的边界面 */ - { - /* 高边界：k0=ex3-1 */ - for (int j0 = 0; j0 < ex2; ++j0) - for (int i0 = 0; i0 < ex1; ++i0) { - const size_t p = idx_ex(i0, j0, ex3 - 1, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - /* 高边界：j0=ex2-1 */ - for (int k0 = 0; k0 < ex3 - 1; ++k0) - for (int i0 = 0; i0 < ex1; ++i0) { - const size_t p = idx_ex(i0, ex2 - 1, k0, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - /* 高边界：i0=ex1-1 */ - for (int k0 = 0; k0 < ex3 - 1; ++k0) - for (int j0 = 0; j0 < ex2 - 1; ++j0) { - const size_t p = idx_ex(ex1 - 1, j0, k0, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - - /* 低边界：当二阶模板也不可用时，对应 i0/j0/k0=0 面 */ - if (kminF == 1) { - for (int j0 = 0; j0 < ex2; ++j0) - for (int i0 = 0; i0 < ex1; ++i0) { - const size_t p = idx_ex(i0, j0, 0, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - } - if (jminF == 1) { - for (int k0 = 0; k0 < ex3; ++k0) - for (int i0 = 0; i0 < ex1; ++i0) { - const size_t p = idx_ex(i0, 0, k0, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - } - if (iminF == 1) { - for (int k0 = 0; k0 < ex3; ++k0) - for (int j0 = 0; j0 < ex2; ++j0) { - const size_t p = idx_ex(0, j0, k0, ex); - fxx[p]=ZEO; fyy[p]=ZEO; fzz[p]=ZEO; - fxy[p]=ZEO; fxz[p]=ZEO; fyz[p]=ZEO; - } - } + const size_t all = (size_t)ex1 * (size_t)ex2 * (size_t)ex3; + for (size_t p = 0; p < all; ++p) { + fxx[p] = ZEO; fxy[p] = ZEO; fxz[p] = ZEO; + fyy[p] = ZEO; fyz[p] = ZEO; fzz[p] = ZEO; } - /* - * 两段式： - * 1) 二阶可用区域先计算二阶模板 - * 2) 高阶可用区域再覆盖四阶模板 - */ - const int i2_lo = (iminF > 0) ? iminF : 0; - const int j2_lo = (jminF > 0) ? jminF : 0; - const int k2_lo = (kminF > 0) ? kminF : 0; - const int i2_hi = ex1 - 2; - const int j2_hi = ex2 - 2; - const int k2_hi = ex3 - 2; - - const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0; - const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0; - const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0; - const int i4_hi = ex1 - 3; - const int j4_hi = ex2 - 3; - const int k4_hi = ex3 - 3; - - /* - * Strategy A: - * Avoid redundant work in overlap of 2nd/4th-order regions. - * Only compute 2nd-order on shell points that are NOT overwritten by - * the 4th-order pass. - */ - const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi); - - if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) { - for (int k0 = k2_lo; k0 <= k2_hi; ++k0) { - const int kF = k0 + 1; - for (int j0 = j2_lo; j0 <= j2_hi; ++j0) { - const int jF = j0 + 1; - for (int i0 = i2_lo; i0 <= i2_hi; ++i0) { - if (has4 && - i0 >= i4_lo && i0 <= i4_hi && - j0 >= j4_lo && j0 <= j4_hi && - k0 >= k4_lo && k0 <= k4_hi) { - continue; - } - const int iF = i0 + 1; - const size_t p = idx_ex(i0, j0, k0, ex); + // Match Fortran (ghost_width=3, "for bam comparison") exactly: + // only compute when x/y/z all satisfy the same-order stencil at this point. + for (int k0 = 0; k0 <= ex3 - 2; ++k0) { + const int kF = k0 + 1; + for (int j0 = 0; j0 <= ex2 - 2; ++j0) { + const int jF = j0 + 1; + for (int i0 = 0; i0 <= ex1 - 2; ++i0) { + const int iF = i0 + 1; + const size_t p = idx_ex(i0, j0, k0, ex); + if ((iF + 2 <= imaxF && iF - 2 >= iminF) && + (jF + 2 <= jmaxF && jF - 2 >= jminF) && + (kF + 2 <= kmaxF && kF - 2 >= kminF)) { + fxx[p] = Fdxdx * ( + -fh[idx_fh_F_ord2(iF - 2, jF, kF, ex)] + + F16 * fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] - + F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - + fh[idx_fh_F_ord2(iF + 2, jF, kF, ex)] + + F16 * fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] + ); + fyy[p] = Fdydy * ( + -fh[idx_fh_F_ord2(iF, jF - 2, kF, ex)] + + F16 * fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] - + F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - + fh[idx_fh_F_ord2(iF, jF + 2, kF, ex)] + + F16 * fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] + ); + fzz[p] = Fdzdz * ( + -fh[idx_fh_F_ord2(iF, jF, kF - 2, ex)] + + F16 * fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] - + F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - + fh[idx_fh_F_ord2(iF, jF, kF + 2, ex)] + + F16 * fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] + ); + fxy[p] = Fdxdy * ( + (fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)] - fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)]) + - F8 * (fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)] - fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)]) + + F8 * (fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)] - fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)]) + - (fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)] - fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)]) + ); + fxz[p] = Fdxdz * ( + (fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)] - fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)]) + - F8 * (fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)] - fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)]) + + F8 * (fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)] - fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)]) + - (fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)] - F8 * fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)] + + F8 * fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)] - fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)]) + ); + fyz[p] = Fdydz * ( + (fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)] - F8 * fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)] + + F8 * fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)] - fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)]) + - F8 * (fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)] - F8 * fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)] + + F8 * fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)] - fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)]) + + F8 * (fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)] - F8 * fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)] + + F8 * fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)] - fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)]) + - (fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)] - F8 * fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)] + + F8 * fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)] - fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)]) + ); + } else if ((iF + 1 <= imaxF && iF - 1 >= iminF) && + (jF + 1 <= jmaxF && jF - 1 >= jminF) && + (kF + 1 <= kmaxF && kF - 1 >= kminF)) { fxx[p] = Sdxdx * ( fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] - TWO * fh[idx_fh_F_ord2(iF, jF, kF, ex)] + fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] ); - fyy[p] = Sdydy * ( fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] - TWO * fh[idx_fh_F_ord2(iF, jF, kF, ex)] + fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] ); - fzz[p] = Sdzdz * ( fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] - TWO * fh[idx_fh_F_ord2(iF, jF, kF, ex)] + fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] ); - fxy[p] = Sdxdy * ( fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)] - fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)] - fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)] + fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)] ); - fxz[p] = Sdxdz * ( fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)] - fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)] - fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)] + fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)] ); - fyz[p] = Sdydz * ( fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)] - fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)] - @@ -207,126 +182,5 @@ void fdderivs(const int ex[3], } } - if (has4) { - for (int k0 = k4_lo; k0 <= k4_hi; ++k0) { - const int kF = k0 + 1; - for (int j0 = j4_lo; j0 <= j4_hi; ++j0) { - const int jF = j0 + 1; - for (int i0 = i4_lo; i0 <= i4_hi; ++i0) { - const int iF = i0 + 1; - const size_t p = idx_ex(i0, j0, k0, ex); - - fxx[p] = Fdxdx * ( - -fh[idx_fh_F_ord2(iF - 2, jF, kF, ex)] + - F16 * fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] - - F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF, ex)] + - F16 * fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] - ); - - fyy[p] = Fdydy * ( - -fh[idx_fh_F_ord2(iF, jF - 2, kF, ex)] + - F16 * fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] - - F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF, ex)] + - F16 * fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] - ); - - fzz[p] = Fdzdz * ( - -fh[idx_fh_F_ord2(iF, jF, kF - 2, ex)] + - F16 * fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] - - F30 * fh[idx_fh_F_ord2(iF, jF, kF, ex)] - - fh[idx_fh_F_ord2(iF, jF, kF + 2, ex)] + - F16 * fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] - ); - - { - const double t_jm2 = - ( fh[idx_fh_F_ord2(iF - 2, jF - 2, kF, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF - 2, kF, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF - 2, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF - 2, kF, ex)] ); - - const double t_jm1 = - ( fh[idx_fh_F_ord2(iF - 2, jF - 1, kF, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF - 1, kF, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF - 1, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF - 1, kF, ex)] ); - - const double t_jp1 = - ( fh[idx_fh_F_ord2(iF - 2, jF + 1, kF, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF + 1, kF, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF + 1, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF + 1, kF, ex)] ); - - const double t_jp2 = - ( fh[idx_fh_F_ord2(iF - 2, jF + 2, kF, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF + 2, kF, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF + 2, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF + 2, kF, ex)] ); - - fxy[p] = Fdxdy * ( t_jm2 - F8 * t_jm1 + F8 * t_jp1 - t_jp2 ); - } - - { - const double t_km2 = - ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 2, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 2, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 2, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF - 2, ex)] ); - - const double t_km1 = - ( fh[idx_fh_F_ord2(iF - 2, jF, kF - 1, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF - 1, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF - 1, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF - 1, ex)] ); - - const double t_kp1 = - ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 1, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 1, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 1, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF + 1, ex)] ); - - const double t_kp2 = - ( fh[idx_fh_F_ord2(iF - 2, jF, kF + 2, ex)] - -F8*fh[idx_fh_F_ord2(iF - 1, jF, kF + 2, ex)] - +F8*fh[idx_fh_F_ord2(iF + 1, jF, kF + 2, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF + 2, ex)] ); - - fxz[p] = Fdxdz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 ); - } - - { - const double t_km2 = - ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 2, ex)] - -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 2, ex)] - +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 2, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF - 2, ex)] ); - - const double t_km1 = - ( fh[idx_fh_F_ord2(iF, jF - 2, kF - 1, ex)] - -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF - 1, ex)] - +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF - 1, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF - 1, ex)] ); - - const double t_kp1 = - ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 1, ex)] - -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 1, ex)] - +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 1, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF + 1, ex)] ); - - const double t_kp2 = - ( fh[idx_fh_F_ord2(iF, jF - 2, kF + 2, ex)] - -F8*fh[idx_fh_F_ord2(iF, jF - 1, kF + 2, ex)] - +F8*fh[idx_fh_F_ord2(iF, jF + 1, kF + 2, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF + 2, ex)] ); - - fyz[p] = Fdydz * ( t_km2 - F8 * t_km1 + F8 * t_kp1 - t_kp2 ); - } - } - } - } - } - // free(fh); } diff --git a/AMSS_NCKU_source/fderivs_c.C b/AMSS_NCKU_source/fderivs_c.C index 90f10a1..69ca48c 100644 --- a/AMSS_NCKU_source/fderivs_c.C +++ b/AMSS_NCKU_source/fderivs_c.C @@ -80,46 +80,48 @@ void fderivs(const int ex[3], fz[p] = ZEO; } - /* - * 两段式： - * 1) 先在二阶可用区域计算二阶模板 - * 2) 再在高阶可用区域覆盖为四阶模板 - * - * 与原 if/elseif 逻辑等价，但减少逐点分支判断。 - */ - const int i2_lo = (iminF > 0) ? iminF : 0; - const int j2_lo = (jminF > 0) ? jminF : 0; - const int k2_lo = (kminF > 0) ? kminF : 0; - const int i2_hi = ex1 - 2; - const int j2_hi = ex2 - 2; - const int k2_hi = ex3 - 2; - - const int i4_lo = (iminF + 1 > 0) ? (iminF + 1) : 0; - const int j4_lo = (jminF + 1 > 0) ? (jminF + 1) : 0; - const int k4_lo = (kminF + 1 > 0) ? (kminF + 1) : 0; - const int i4_hi = ex1 - 3; - const int j4_hi = ex2 - 3; - const int k4_hi = ex3 - 3; - - if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) { - for (int k0 = k2_lo; k0 <= k2_hi; ++k0) { - const int kF = k0 + 1; - for (int j0 = j2_lo; j0 <= j2_hi; ++j0) { - const int jF = j0 + 1; - for (int i0 = i2_lo; i0 <= i2_hi; ++i0) { - const int iF = i0 + 1; - const size_t p = idx_ex(i0, j0, k0, ex); + // Match Fortran (ghost_width=3, "for bam comparison") exactly: + // only compute when x/y/z all satisfy the same-order stencil at this point. + for (int k0 = 0; k0 <= ex3 - 2; ++k0) { + const int kF = k0 + 1; + for (int j0 = 0; j0 <= ex2 - 2; ++j0) { + const int jF = j0 + 1; + for (int i0 = 0; i0 <= ex1 - 2; ++i0) { + const int iF = i0 + 1; + const size_t p = idx_ex(i0, j0, k0, ex); + if ((iF + 2 <= imaxF && iF - 2 >= iminF) && + (jF + 2 <= jmaxF && jF - 2 >= jminF) && + (kF + 2 <= kmaxF && kF - 2 >= kminF)) { + fx[p] = d12dx * ( + fh[idx_fh_F_ord2(iF - 2, jF, kF, ex)] - + EIT * fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] + + EIT * fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] - + fh[idx_fh_F_ord2(iF + 2, jF, kF, ex)] + ); + fy[p] = d12dy * ( + fh[idx_fh_F_ord2(iF, jF - 2, kF, ex)] - + EIT * fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] + + EIT * fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] - + fh[idx_fh_F_ord2(iF, jF + 2, kF, ex)] + ); + fz[p] = d12dz * ( + fh[idx_fh_F_ord2(iF, jF, kF - 2, ex)] - + EIT * fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] + + EIT * fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] - + fh[idx_fh_F_ord2(iF, jF, kF + 2, ex)] + ); + } else if ((iF + 1 <= imaxF && iF - 1 >= iminF) && + (jF + 1 <= jmaxF && jF - 1 >= jminF) && + (kF + 1 <= kmaxF && kF - 1 >= kminF)) { fx[p] = d2dx * ( -fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] + fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] ); - fy[p] = d2dy * ( -fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] + fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] ); - fz[p] = d2dz * ( -fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] + fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] @@ -129,39 +131,5 @@ void fderivs(const int ex[3], } } - if (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi) { - for (int k0 = k4_lo; k0 <= k4_hi; ++k0) { - const int kF = k0 + 1; - for (int j0 = j4_lo; j0 <= j4_hi; ++j0) { - const int jF = j0 + 1; - for (int i0 = i4_lo; i0 <= i4_hi; ++i0) { - const int iF = i0 + 1; - const size_t p = idx_ex(i0, j0, k0, ex); - - fx[p] = d12dx * ( - fh[idx_fh_F_ord2(iF - 2, jF, kF, ex)] - - EIT * fh[idx_fh_F_ord2(iF - 1, jF, kF, ex)] + - EIT * fh[idx_fh_F_ord2(iF + 1, jF, kF, ex)] - - fh[idx_fh_F_ord2(iF + 2, jF, kF, ex)] - ); - - fy[p] = d12dy * ( - fh[idx_fh_F_ord2(iF, jF - 2, kF, ex)] - - EIT * fh[idx_fh_F_ord2(iF, jF - 1, kF, ex)] + - EIT * fh[idx_fh_F_ord2(iF, jF + 1, kF, ex)] - - fh[idx_fh_F_ord2(iF, jF + 2, kF, ex)] - ); - - fz[p] = d12dz * ( - fh[idx_fh_F_ord2(iF, jF, kF - 2, ex)] - - EIT * fh[idx_fh_F_ord2(iF, jF, kF - 1, ex)] + - EIT * fh[idx_fh_F_ord2(iF, jF, kF + 1, ex)] - - fh[idx_fh_F_ord2(iF, jF, kF + 2, ex)] - ); - } - } - } - } - // free(fh); }