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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 523063d030b5aaf53dc2e914cef38f39767df998
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu Mar 1 07:54:14 2007 +1100

    VOF normal is computed using "myc" in 2D (thanks to Ruben Scardovelli)
    
    darcs-hash:20070228205414-d4795-4ccee5f52a314cba12ed25ad01168a0cd7c65ddf.gz

diff --git a/src/Makefile.am b/src/Makefile.am
index 224a48b..0654a04 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -87,6 +87,7 @@ SRC = \
 	ocean.c \
 	levelset.c \
 	myc.h \
+	myc2d.h \
 	$(GFS_HDS)
 
 libgfs3D_la_LDFLAGS = $(NO_UNDEFINED)\
diff --git a/src/myc2d.h b/src/myc2d.h
new file mode 100644
index 0000000..534a1bb
--- /dev/null
+++ b/src/myc2d.h
@@ -0,0 +1,68 @@
+#define NOT_ZERO 1.e-30
+
+/*-----------------------------------------------------*
+ *MYC - Mixed Youngs and Central Scheme (2D)           *
+ *-----------------------------------------------------*/
+static void mycs(double c[3][3],double mxy[2])
+{
+  int ix;
+  double c_t,c_b,c_r,c_l;
+  double mx0,my0,mx1,my1,mm1,mm2;
+  
+  /* top, bottom, right and left sums of c values */
+  c_t = c[0][2] + c[1][2] + c[2][2];
+  c_b = c[0][0] + c[1][0] + c[2][0];
+  c_r = c[2][0] + c[2][1] + c[2][2];
+  c_l = c[0][0] + c[0][1] + c[0][2];
+
+  /* consider two lines: sgn(my) Y =  mx0 X + alpha,
+     and: sgn(mx) X =  my0 Y + alpha */ 
+  mx0 = 0.5*(c_l-c_r);
+  my0 = 0.5*(c_b-c_t);
+
+  /* minimum coefficient between mx0 and my0 wins */
+  if (fabs(mx0) <= fabs(my0)) {
+    mm1 = fabs(my0) + NOT_ZERO; 
+    my0 = my0/mm1;
+    ix = 1;
+  }
+  else {
+    mm1 = fabs(mx0) + NOT_ZERO; 
+    mx0 = mx0/mm1;
+    ix = 0;
+  }
+
+  /* Youngs' normal to the interface */
+  mm1 = c[0][0] + 2.0*c[0][1] + c[0][2];
+  mm2 = c[2][0] + 2.0*c[2][1] + c[2][2];
+  mx1 = mm1 - mm2;
+  mm1 = c[0][0] + 2.0*c[1][0] + c[2][0];
+  mm2 = c[0][2] + 2.0*c[1][2] + c[2][2];
+  my1 = mm1 - mm2;
+
+  /* choose between the best central and Youngs' scheme */ 
+  if (ix) {
+    mm1 = fabs(my1) + NOT_ZERO; 
+    mm1 = fabs(mx1)/mm1;
+    if (mm1 > fabs(mx0)) {
+      mx0 = mx1;
+      my0 = my1;
+    }
+  }
+  else {
+    mm1 = fabs(mx1) + NOT_ZERO; 
+    mm1 = fabs(my1)/mm1;
+    if (mm1 > fabs(my0)) {
+      mx0 = mx1;
+      my0 = my1;
+    }
+  }
+	
+  /* normalize the set (mx0,my0): |mx0|+|my0|=1 and
+     write the two components of the normal vector  */
+  mm1 = fabs(mx0) + fabs(my0) + NOT_ZERO; 
+  mxy[0] = mx0/mm1;
+  mxy[1] = my0/mm1;
+  
+  return;
+}
diff --git a/src/vof.c b/src/vof.c
index 29a21a2..a5b57c2 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -473,14 +473,20 @@ static FttCell * domain_and_boundary_locate (GfsDomain * domain, FttVector p, gu
   return gfs_domain_boundary_locate (domain, p, level);
 }
 
-static void stencil (FttCell * cell, GfsVariable * v, gdouble f[3][3][3])
+#if FTT_2D
+# define F(x,y,z) f[x][y]
+#else
+# define F(x,y,z) f[x][y][z]
+#endif
+
+static void stencil (FttCell * cell, GfsVariable * v, gdouble F(3,3,3))
 {
   gdouble h = ftt_cell_size (cell);
   guint level = ftt_cell_level (cell);
   FttVector p;
   gint x, y, z = 0;
   
-  f[1][1][1] = GFS_VARIABLE (cell, v->i);
+  F(1,1,1) = GFS_VARIABLE (cell, v->i);
   ftt_cell_pos (cell, &p);
 #if !FTT_2D
   for (z = -1; z <= 1; z++)
@@ -492,20 +498,20 @@ static void stencil (FttCell * cell, GfsVariable * v, gdouble f[3][3][3])
 	  o.x = p.x + h*x; o.y = p.y + h*y; o.z = p.z + h*z;
 	  FttCell * neighbor = domain_and_boundary_locate (v->domain, o, level);
 	  if (neighbor)
-	    f[x + 1][y + 1][z + 1] = 
+	    F(x + 1, y + 1, z + 1) =
 	      gfs_vof_interpolate (neighbor, &o, level, GFS_VARIABLE_TRACER_VOF (v));
 	  else
-	    f[x + 1][y + 1][z + 1] = -1.;
+	    F(x + 1, y + 1, z + 1) = -1.;
 	}
   /* boundary conditions (symmetry) */
 #if FTT_2D
   for (x = 0; x <= 2; x++) {
-    if (f[x][0][1] < 0.) f[x][0][1] = f[x][1][1];
-    if (f[x][2][1] < 0.) f[x][2][1] = f[x][1][1];
+    if (f[x][0] < 0.) f[x][0] = f[x][1];
+    if (f[x][2] < 0.) f[x][2] = f[x][1];
   }
   for (y = 0; y <= 2; y++) {
-    if (f[0][y][1] < 0.) f[0][y][1] = f[1][y][1];
-    if (f[2][y][1] < 0.) f[2][y][1] = f[1][y][1];
+    if (f[0][y] < 0.) f[0][y] = f[1][y];
+    if (f[2][y] < 0.) f[2][y] = f[1][y];
   }
 #else /* 3D */
   for (x = 0; x <= 2; x++)
@@ -528,12 +534,12 @@ static void stencil (FttCell * cell, GfsVariable * v, gdouble f[3][3][3])
 
 static void youngs_normal (FttCell * cell, GfsVariable * v, FttVector * n)
 {
-  gdouble f[3][3][3];
+  gdouble F(3,3,3);
 
   stencil (cell, v, f);
 #if FTT_2D
-  n->x = (f[0][2][1] + 2.*f[0][1][1] + f[0][0][1] - 2.*f[2][1][1] - f[2][2][1] - f[2][0][1])/8.;
-  n->y = (f[2][0][1] + 2.*f[1][0][1] + f[0][0][1] - 2.*f[1][2][1] - f[2][2][1] - f[0][2][1])/8.;
+  n->x = (f[0][2] + 2.*f[0][1] + f[0][0] - 2.*f[2][1] - f[2][2] - f[2][0])/8.;
+  n->y = (f[2][0] + 2.*f[1][0] + f[0][0] - 2.*f[1][2] - f[2][2] - f[0][2])/8.;
   n->z = 0.;
 #else  /* 3D */
   gdouble mm1 = f[0][0][0] + f[0][0][2] + f[0][2][0] + f[0][2][2] +
@@ -562,21 +568,21 @@ static void youngs_normal (FttCell * cell, GfsVariable * v, FttVector * n)
 #endif /* 3D */
 }
 
-#include "myc.h"
+#if FTT_2D
+# include "myc2d.h"
+#else
+# include "myc.h"
+#endif
 
 static void myc_normal (FttCell * cell, GfsVariable * v, FttVector * n)
 {
-  gdouble f[3][3][3];
+  gdouble F(3,3,3);  
 
   stencil (cell, v, f);
+  mycs (f, &n->x);
 #if FTT_2D
-  /* fixme: this is still Youngs */
-  n->x = (f[0][2][1] + 2.*f[0][1][1] + f[0][0][1] - 2.*f[2][1][1] - f[2][2][1] - f[2][0][1])/8.;
-  n->y = (f[2][0][1] + 2.*f[1][0][1] + f[0][0][1] - 2.*f[1][2][1] - f[2][2][1] - f[0][2][1])/8.;
   n->z = 0.;
-#else  /* 3D */
-  mycs (f, &n->x);
-#endif /* 3D */
+#endif
 }
 
 static void vof_plane (FttCell * cell, GfsVariable * v)

-- 
Gerris Flow Solver