[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 97e1c87a4a46e9d4547c1680d892cc4f0a3fb936
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Wed May 23 12:42:52 2007 +1000

    Bug fix for marginal case in myc (3D version)
    
    darcs-hash:20070523024252-d4795-6ef8f9f5ec97dc7f9659254bedcd7478621dabf9.gz

diff --git a/src/myc.h b/src/myc.h
index 1cf061c..ea4b317 100644
--- a/src/myc.h
+++ b/src/myc.h
@@ -25,8 +25,7 @@ static void mycs(double c[3][3][3],double mxyz[3])
        c[0][1][1];
   m2 = c[2][1][0] + c[2][1][2] + c[2][0][1] + c[2][2][1] + 
        c[2][1][1];
-  t0 = fabs(m1-m2) + NOT_ZERO; 
-  m[0][0] = (m1-m2)/t0;
+  m[0][0] = m1 > m2 ? 1. : -1.;
 
   m1 = c[0][0][1]+ c[2][0][1]+ c[1][0][1];
   m2 = c[0][2][1]+ c[2][2][1]+ c[1][2][1];
@@ -46,8 +45,7 @@ static void mycs(double c[3][3][3],double mxyz[3])
        c[1][0][1];
   m2 = c[1][2][0] + c[1][2][2] + c[2][2][1] + c[0][2][1] +
        c[1][2][1];
-  t0 = fabs(m1-m2) + NOT_ZERO; 
-  m[1][1] = (m1-m2)/t0;
+  m[1][1] = m1 > m2 ? 1. : -1.;
 
   m1 = c[1][0][0]+ c[1][1][0]+ c[1][2][0];
   m2 = c[1][0][2]+ c[1][1][2]+ c[1][2][2];
@@ -68,21 +66,20 @@ static void mycs(double c[3][3][3],double mxyz[3])
        c[1][1][0];
   m2 = c[0][1][2] + c[2][1][2] + c[1][0][2] + c[1][2][2] +
        c[1][1][2];
-  t0 = fabs(m1-m2) + NOT_ZERO; 
-  m[2][2] = (m1-m2)/t0;
+  m[2][2] = m1 > m2 ? 1. : -1.;
 
   /* normalize each set (mx,my,mz): |mx|+|my|+|mz| = 1 */
-  t0 = fabs(m[0][0]) + fabs(m[0][1]) + fabs(m[0][2]) + NOT_ZERO;
+  t0 = fabs(m[0][0]) + fabs(m[0][1]) + fabs(m[0][2]);
   m[0][0] /= t0;
   m[0][1] /= t0;
   m[0][2] /= t0;
 
-  t0 = fabs(m[1][0]) + fabs(m[1][1]) + fabs(m[1][2]) + NOT_ZERO;
+  t0 = fabs(m[1][0]) + fabs(m[1][1]) + fabs(m[1][2]);
   m[1][0] /= t0;
   m[1][1] /= t0;
   m[1][2] /= t0;
 
-  t0 = fabs(m[2][0]) + fabs(m[2][1]) + fabs(m[2][2]) + NOT_ZERO;
+  t0 = fabs(m[2][0]) + fabs(m[2][1]) + fabs(m[2][2]);
   m[2][0] /= t0;
   m[2][1] /= t0;
   m[2][2] /= t0;
@@ -93,9 +90,11 @@ static void mycs(double c[3][3][3],double mxyz[3])
   t2 = fabs(m[2][2]);
 
   cn = 0;
-  if (t1 > t0 && t1 > t2)
+  if (t1 > t0) {
+    t0 = t1;
     cn = 1;
-  if (t2 > t0 && t2 > t1)
+  }
+  if (t2 > t0)
     cn = 2;
 
   /* Youngs-CIAM scheme */  
@@ -130,10 +129,15 @@ static void mycs(double c[3][3][3],double mxyz[3])
   m[3][2] /= t0;
 
   /* choose between the previous choice and Youngs-CIAM */
-  t0 = MAX(fabs(m[3][0]),fabs(m[3][1]));
-  t0 = MAX(t0,fabs(m[3][2]));
-
-  if (fabs(m[cn][cn]) > t0) 
+  t0 = fabs (m[3][0]);
+  t1 = fabs (m[3][1]);
+  t2 = fabs (m[3][2]);
+  if (t1 > t0)
+    t0 = t1;
+  if (t2 > t0)
+    t0 = t2;
+
+  if (fabs(m[cn][cn]) > t0)
     cn = 3;
 
   /* components of the normal vector */
diff --git a/src/myc2d.h b/src/myc2d.h
index 534a1bb..7caf540 100644
--- a/src/myc2d.h
+++ b/src/myc2d.h
@@ -22,13 +22,11 @@ static void mycs(double c[3][3],double mxy[2])
 
   /* minimum coefficient between mx0 and my0 wins */
   if (fabs(mx0) <= fabs(my0)) {
-    mm1 = fabs(my0) + NOT_ZERO; 
-    my0 = my0/mm1;
+    my0 = my0 > 0. ? 1. : -1.;
     ix = 1;
   }
   else {
-    mm1 = fabs(mx0) + NOT_ZERO; 
-    mx0 = mx0/mm1;
+    mx0 = mx0 > 0. ? 1. : -1.;
     ix = 0;
   }
 

-- 
Gerris Flow Solver