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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit ba0080acfae1883be592065df53fe67aaf872f5a
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu Nov 2 15:27:08 2006 +1100

    3D implementation of gfs_youngs_normal()
    
    darcs-hash:20061102042708-d4795-87336ddc2e7b35d19f705fa5f1640b9aada1f22d.gz

diff --git a/src/vof.c b/src/vof.c
index 5b2f07f..960ff0f 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -426,44 +426,47 @@ 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])
+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;
+  gint x, y, z = 0;
   
-  f[1][1] = GFS_VARIABLE (cell, v->i);
+  f[1][1][1] = GFS_VARIABLE (cell, v->i);
   ftt_cell_pos (cell, &p);
-  for (x = -1; x <= 1; x++)
-    for (y = -1; y <= 1; y++)
-      if (x != 0 || y != 0) {
-	FttVector o;
-	o.x = p.x + h*x; o.y = p.y + h*y; o.z = 0.;
-	FttCell * neighbor = domain_and_boundary_locate (v->domain, o, level);
-	if (!neighbor) /* fixme: boundary conditions */
-	  f[x + 1][y + 1] = f[1][1];
-	else {
-	  guint l = ftt_cell_level (neighbor);
-	  FttVector m;
-	  gdouble alpha;
-	  if (l == level || !gfs_vof_plane (neighbor, v, &m, &alpha))
-	    f[x + 1][y + 1] = GFS_VARIABLE (neighbor, v->i);
+#if !FTT_2D
+  for (z = -1; z <= 1; z++)
+#endif
+    for (x = -1; x <= 1; x++)
+      for (y = -1; y <= 1; y++)
+	if (x != 0 || y != 0 || z != 0) {
+	  FttVector o;
+	  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) /* fixme: boundary conditions */
+	    f[x + 1][y + 1][z + 1] = f[1][1][1];
 	  else {
-	    FttComponent c;
-	    FttVector q;
-	    
-	    g_assert (l == level - 1);
-	    ftt_cell_pos (neighbor, &q);
-	    for (c = 0; c < FTT_DIMENSION; c++) {
-	      gdouble a = ((&o.x)[c] - (&q.x)[c])/h;
-	      g_assert (fabs (a) == 0.5);
-	      alpha -= (&m.x)[c]*(0.25 - a/2.);
+	    guint l = ftt_cell_level (neighbor);
+	    FttVector m;
+	    gdouble alpha;
+	    if (l == level || !gfs_vof_plane (neighbor, v, &m, &alpha))
+	      f[x + 1][y + 1][z + 1] = GFS_VARIABLE (neighbor, v->i);
+	    else {
+	      FttComponent c;
+	      FttVector q;
+	      
+	      g_assert (l == level - 1);
+	      ftt_cell_pos (neighbor, &q);
+	      for (c = 0; c < FTT_DIMENSION; c++) {
+		gdouble a = ((&o.x)[c] - (&q.x)[c])/h;
+		g_assert (fabs (a) == 0.5);
+		alpha -= (&m.x)[c]*(0.25 - a/2.);
+	      }
+	      f[x + 1][y + 1][z + 1] = gfs_plane_volume (&m, 2.*alpha);
 	    }
-	    f[x + 1][y + 1] = gfs_plane_volume (&m, 2.*alpha);
 	  }
 	}
-      }  
 }
 
 /**
@@ -477,37 +480,40 @@ static void stencil (FttCell * cell, GfsVariable * v, gdouble f[3][3])
  */
 void gfs_youngs_normal (FttCell * cell, GfsVariable * v, FttVector * n)
 {
-  gdouble f[3][3];
+  gdouble f[3][3][3];
 
   stencil (cell, v, f);
-  n->x = (2.*f[2][1] + f[2][2] + f[2][0] - f[0][2] - 2.*f[0][1] - f[0][0])/8.;
-  n->y = (2.*f[1][2] + f[2][2] + f[0][2] - f[2][0] - 2.*f[1][0] - f[0][0])/8.;
-  n->z = 0.;
-} 
-
-static void column_normal (gdouble f[3][3], FttVector * n)
-{
-  gdouble s1, s2;
-
-  s1 = f[2][0] + f[2][1] + f[2][2] - f[0][0] - f[0][1] - f[0][2];
-  s2 = f[0][2] + f[1][2] + f[2][2] - f[0][0] - f[1][0] - f[2][0];
-  if (fabs (s2) > fabs (s1)) {
-    n->x = s1; n->y = s2 < 0. ? -2. : 2.;
-  }
-  else {
-    n->y = s2; n->x = s1 < 0. ? - 2. : 2.;
-  }
+#if FTT_2D
+  n->x = (2.*f[2][1][1] + f[2][2][1] + f[2][0][1] - f[0][2][1] - 2.*f[0][1][1] - f[0][0][1])/8.;
+  n->y = (2.*f[1][2][1] + f[2][2][1] + f[0][2][1] - f[2][0][1] - 2.*f[1][0][1] - f[0][0][1])/8.;
   n->z = 0.;
+#else  /* 3D */
+  gdouble mm1 = f[0][0][0] + f[0][0][2] + f[0][2][0] + f[0][2][2] +
+    2.*(f[0][0][1] + f[0][2][1] + f[0][1][0] + f[0][1][2]) + 
+    4.*f[0][1][1];
+  gdouble mm2 = f[2][0][0] + f[2][0][2] + f[2][2][0] + f[2][2][2] + 
+    2.*(f[2][0][1] + f[2][2][1] + f[2][1][0] + f[2][1][2]) + 
+    4.*f[2][1][1];
+  n->x = (mm2 - mm1)/32.;
+    
+  mm1 = f[0][0][0] + f[0][0][2] + f[2][0][0] + f[2][0][2] + 
+    2.*(f[0][0][1] + f[2][0][1] + f[1][0][0] + f[1][0][2]) + 
+    4.*f[1][0][1];
+  mm2 = f[0][2][0] + f[0][2][2] + f[2][2][0] + f[2][2][2] + 
+    2.*(f[0][2][1] + f[2][2][1] + f[1][2][0] + f[1][2][2]) + 
+    4.*f[1][2][1];
+  n->y = (mm2 - mm1)/32.;
+                  
+  mm1 = f[0][0][0] + f[0][2][0] + f[2][0][0] + f[2][2][0] +
+    2.*(f[0][1][0] + f[2][1][0] + f[1][0][0] + f[1][2][0]) + 
+    4.*f[1][1][0];
+  mm2 = f[0][0][2] + f[0][2][2] + f[2][0][2] + f[2][2][2] + 
+    2.*(f[0][1][2] + f[2][1][2] + f[1][0][2] + f[1][2][2]) + 
+    4.*f[1][1][2];
+  n->z = (mm2 - mm1)/32.;
+#endif /* 3D */
 }
 
-void gfs_column_normal (FttCell * cell, GfsVariable * v, FttVector * n)
-{
-  gdouble f[3][3];
-
-  stencil (cell, v, f);
-  column_normal (f, n);
-} 
-
 static gdouble plane_volume_shifted (FttVector m, gdouble alpha, FttVector p[2])
 {
   FttComponent c;

-- 
Gerris Flow Solver