Finish trilinear case 4

PSORLab · Nov 6, 2020 · cbe7bc4 · cbe7bc4
1 parent 30ac33c
commit cbe7bc4
Showing 1 changed file with 38 additions and 35 deletions.
diff --git a/src/forward_operators/trilinear.jl b/src/forward_operators/trilinear.jl
@@ -335,24 +335,25 @@ function trilinear_case_4(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float6
         θcv3 = -0.5*(xyzUUL + xyzLLL - xyzULU - xyzLUU)/delX
         θcv4 = xyzLLL - θcv1*xL - θcv2*yL - θcv3*zL
 
-        # TODO: CHECK OTHER PAPER
-        #cv_b1 = -2.0*xyzUUU
-        #cv_b2 = -2.0*xyzULL
-        #cv_b3 = -2.0*xyzUUU
-        #cv_b4 = -2.0*xyzUUU
-        #cv_b5 = θcv4
-        #cv_b6 = TODO not in paper
-
-        #cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUU, yzLL, yzUL, yzUL, yzLU, yzLU
-        #cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUU, xzUL, xzLU, xzLL, -θcv1/delY, xzLU
-        #cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUU, xyUL, xyLU, xyLL, xyLL, θcv2/delZ
-
-        #cv1 = cv_ax1*x.cv + cv_ay1*y.cv + cv_az1*z.cv + cv_b1
-        #cv2 = cv_ax2*x.cv - cv_ay2*y.cc - cv_az2*z.cc + cv_b2
-        #cv3 = -cv_ax3*x.cc + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
-        #cv4 = -cv_ax4*x.cc - cv_ay4*y.cc - cv_az4*z.cc + cv_b4
-        #cv5 = -cv_ax5*x.cc + cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, -y.cc) - cv_az5*z.cc + cv_b5
-        #cv6 = -cv_ax6*x.cc + cv_ay6*y.cv + cv_az6*ifelse(cv_az6 > 0.0, z.cv, -z.cc) + cv_b6
+        cv_b1 = -2.0*xyzLUL
+        cv_b2 = -2.0*xyzLLU
+        cv_b3 = -2.0*xyzUUU
+        cv_b4 = -2.0*xyzULL
+        cv_b5 = θcv4
+        #cv_b6 = TODO? eq 6 not in paper
+
+        cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5 = yzUL, yzLU, yzUU, yzLL, θcv1  # TODO? eq 6 not in paper
+        cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5 = xzLL, xzLU, xzUU, xzUL, θcv2  # TODO? eq 6 not in paper
+        cv_az1, cv_az2, cv_az3, cv_az4, cv_az5 = xyLU, xyLL, xyUU, xyUL, θcv3  # TODO? eq 6 not in paper
+
+        cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
+        cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
+        cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
+        cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cc + cv_b4
+        cv5 = cv_ax5*ifelse(cv_ax5 > 0.0, x.cv, -x.cc) +
+              cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, -y.cc) +
+              cv_az5*ifelse(cv_az5 > 0.0, z.cv, -z.cc) + cv_b5
+        cv6 = -Inf # TODO? eq 6 not in paper
 
     elseif xyzUUL + xyzULU + xyzLUU >= xyzLLL + 2*xyzUUU
 
@@ -414,23 +415,25 @@ function trilinear_case_4(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float6
        θcc3 = 0.5*(xyzLLU + xyzUUU - xyzLUL - xyzULL)/delZ
        θcc4 = xyzUUU - θcc1*xU - θcc2*yU - θcc3*zU
 
-       #cc_b1 = -2.0*xyzUUL
-       #cc_b2 = -(xyzULU + xyzULL)
-       #cc_b3 = -(xyzLUU + xyzLLU)
-       #cc_b4 = -(xyzLUU + xyzLUL)
-       #cc_b5 = -(xyzULU + xyzLLU)
-       #cc_b6 = θcc*zU/delZ - xyzULL - xyzLUL + xyzUUU
-
-       #cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzUL, yzLL, yzUU, yzUU, yzLU, yzLL
-       #cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUL, xzUU, xzLU, xzLL, xzUU, xzLL
-       #cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUU, xyUL, xyLL, xyLU, xyLL, -θcc/delZ
-
-       #cc1 = -cc_ax1*x.cv - cc_ay1*y.cv + cc_az1*z.cc + cc_b1
-       #cc2 = -cc_ax2*x.cv + cc_ay2*y.cc + cc_az2*z.cv + cc_b2
-       #cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
-       #cc4 = cc_ax4*x.cc - cc_ay4*y.cv + cc_az4*z.cc + cc_b4
-       #cc5 = cc_ax5*x.cc + cc_ay5*y.cc + cc_az5*z.cv + cc_b5
-       #cc6 = -cc_ax6*x.cv - cc_ay6*y.cv + cc_az6*ifelse(cc_az6 > 0.0, z.cc, -z.cv) + cc_b6
+       cc_b1 = -2.0*xyzLLL
+       cc_b2 = -2.0*xyzUUL
+       cc_b3 = -2.0*xyzULU
+       cc_b4 = -2.0*xyzLUU
+       cc_b5 = θcc4
+       #cc_b6 = TODO? eq 6 not in paper
+
+       cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5 = yzLL, yzUL, yzLU, yzUU, θcc1
+       cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5 = xzLL, xzUL, xzUU, xzLU, θcc2
+       cc_az1, cc_az2, cc_az3, cc_az4, cc_az5 = xyLL, xyUU, xyUL, xyLU, θcc3
+
+       cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
+       cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
+       cc3 = cc_ax3*x.cv + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
+       cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
+       cc5 = cc_ax5*ifelse(cc_ax5 > 0.0, x.cc, -x.cv) +
+             cc_ay5*ifelse(cc_ay5 > 0.0, y.cc, -y.cv) +
+             cc_az5*ifelse(cc_az5 > 0.0, z.cc, -z.cv) + cc_b5
+       cc6 = Inf # TODO? eq 6 not in paper
 
    elseif (xyzULL + xyzLUL + xyzLLU) <= (xyzUUU + 2*xyzLLL)