diff options
| author | Christian Reisswig <reisswig@tapir.caltech.edu> | 2011-06-06 16:40:52 -0700 |
|---|---|---|
| committer | Barry Wardell <barry.wardell@gmail.com> | 2011-12-14 18:26:22 +0000 |
| commit | 1801bc6589abec293dd80003740539dec71d74d3 (patch) | |
| tree | a9a195d88e2e64a62a65be148c71fda1985bbdf6 /Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc | |
| parent | af68e28f01bff4772eac13b4569de02bb11e222b (diff) | |
CC-ENO: Fixed indexing problem in first order failsafe interpolation.
Diffstat (limited to 'Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc')
| -rw-r--r-- | Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc | 31 |
1 files changed, 22 insertions, 9 deletions
diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc index bbd742f40..82ab24c48 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc +++ b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc @@ -289,7 +289,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (lV*rV < 0) + if (lV*rV <= 0) { // switch back to first order! res = 0; @@ -338,7 +338,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] < 0) + if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -389,6 +389,10 @@ namespace CarpetLib { } } + /*typedef coeffs1d<RT,1,di> coeffs; + for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) { + res += coeffs::get(i) * interp0<T,1> (p + i*d1); + }*/ return res; } @@ -419,7 +423,6 @@ namespace CarpetLib { { static_assert (dj==0 or dj==1, "dj must be 0 or 1"); typedef typename typeprops<T>::real RT; - typedef coeffs1d<RT,ORDER,dj> coeffs; T res = typeprops<T>::fromreal (0); // get function values needed for the stencil T f[coeffs1d<RT,ORDER,dj>::maxncoeffs]; @@ -446,7 +449,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (lV*rV < 0) + if (lV*rV <= 0) { // switch back to first order! res = 0; @@ -493,7 +496,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] < 0) + if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -543,6 +546,11 @@ namespace CarpetLib { res += coeffs1::get(i) * f[i-coeffs1d<RT,ORDER,dj>::minimin]; } } + + /*typedef coeffs1d<RT,1,dj> coeffs; + for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) { + res += coeffs::get(i) * interp1<T,1,di> (p + i*d2, d1); + }*/ return res; } @@ -575,7 +583,6 @@ namespace CarpetLib { { static_assert (dk==0 or dk==1, "dk must be 0 or 1"); typedef typename typeprops<T>::real RT; - typedef coeffs1d<RT,ORDER,dk> coeffs; T res = typeprops<T>::fromreal (0); // get function values needed for the stencil T f[coeffs1d<RT,ORDER,dk>::maxncoeffs]; @@ -602,7 +609,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (lV*rV < 0) + if (lV*rV <= 0) { // switch back to first order! res = 0; @@ -649,7 +656,7 @@ namespace CarpetLib { } // check that divided differences do not change sign: if so go back to first order! - if (V[0]*V[2] < 0) + if (V[0]*V[2] <= 0) { // switch back to first order! res = 0; @@ -690,7 +697,7 @@ namespace CarpetLib { } // check that result is reasonable! - if ((res - f[-coeffs1d<RT,ORDER,dj>::minimin-1+dk]) * (f[-coeffs1d<RT,ORDER,dk>::minimin+dk] - res) < 0) + if ((res - f[-coeffs1d<RT,ORDER,dk>::minimin-1+dk]) * (f[-coeffs1d<RT,ORDER,dk>::minimin+dk] - res) < 0) { res = 0; // switch back to first order @@ -700,6 +707,11 @@ namespace CarpetLib { } } + /*typedef coeffs1d<RT,1,dk> coeffs; + for (ptrdiff_t i=coeffs::imin; i<coeffs::imax; ++i) { + res += coeffs::get(i) * interp2<T,1,di,dj> (p + i*d3, d1, d2); + }*/ + return res; } @@ -1160,6 +1172,7 @@ namespace CarpetLib { #define TYPECASE(N,T) \ \ \ + \ template \ void \ prolongate_3d_cc_eno_rf2<T,2> (T const * restrict const src, \ |