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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 397002761f3f2cdae316bea6d48f09bb19b5705a
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Sun Mar 5 18:59:59 2006 +1100

    Fix for special cases when computing solid fractions
    
    Solid fractions of cells cut more than once by the surface were not
    consistently computed by the approximate "VOF" technique. This was causing
    convergence problems in some cases.
    
    These cells (rare) are now treated using a simple approximation.
    
    darcs-hash:20060305075959-d4795-326b70715df4c283f783b0167e6b35b77ce02aa6.gz

diff --git a/src/isocube.h b/src/isocube.h
index 0093dec..2c3819f 100644
--- a/src/isocube.h
+++ b/src/isocube.h
@@ -44,3 +44,17 @@ static guint face_v[6][4] = {
   {1,5,7,3},/* front */
   {0,4,6,2} /* back */
 };
+static guint connect[12][2][3] = {
+  {{9, 1, 8}, {4, 3, 7}},   /* 0 */
+  {{6, 2, 5}, {8, 0, 9}},   /* 1 */
+  {{10, 3, 11}, {5, 1, 6}}, /* 2 */
+  {{7, 0, 4}, {11, 2, 10}}, /* 3 */
+  {{3, 7, 0}, {8, 5, 11}},  /* 4 */
+  {{11, 4, 8}, {1, 6, 2}},  /* 5 */
+  {{2, 5, 1}, {9, 7, 10}},  /* 6 */
+  {{10, 6, 9}, {0, 4, 3}},  /* 7 */
+  {{5, 11, 4}, {0, 9, 1}},  /* 8 */
+  {{1, 8, 0}, {7, 10, 6}},  /* 9 */
+  {{6, 9, 7}, {3, 11, 2}},  /* 10 */
+  {{2, 10, 3}, {4, 8, 5}}   /* 11 */
+};
diff --git a/src/solid.c b/src/solid.c
index e2ebde7..8399c45 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -340,6 +340,28 @@ static void cell_size (FttCell * cell, FttVector * h)
 #endif /* 3D */
 }
 
+/* Returns: the number of closed loops for the given isocube */
+static guint topology (CellCube * cube)
+{
+  guint l, nl = 0;
+
+  for (l = 0; l < 12; l++) {
+    guint nv = 0, e = l;
+    
+    while (cube->n[e]) {
+      guint m = 0, * ne = connect[e][cube->inside[e] > 0];
+
+      nv++;
+      cube->n[e] = 0;
+      while (m < 3 && !cube->n[e])
+	e = ne[m++];
+    }
+    if (nv > 2)
+      nl++;
+  }
+  return nl;
+}
+
 static void set_solid_fractions_from_surface3D (FttCell * cell, GtsSurface * s)
 {
   GfsSolidVector * solid = GFS_STATE (cell)->solid;
@@ -450,13 +472,10 @@ static void set_solid_fractions_from_surface3D (FttCell * cell, GtsSurface * s)
   }
 
   /* now compute cell fraction, center of area, center of mass */
-
-  /* fixme: the calculation of ca below is not strictly equivalent to
-     the true center of area of the piece of surface contained in the
-     cell if there are more than 4 intersection points (n1 > 4) */
   ca.x /= n1; ca.y /= n1; ca.z /= n1; 
   solid->ca = ca;
-  {
+  switch (topology (&cube)) {
+  case 1: {
     FttVector m;
     gdouble alpha, n = 0.;
     gboolean sym[FTT_DIMENSION];
@@ -474,25 +493,29 @@ static void set_solid_fractions_from_surface3D (FttCell * cell, GtsSurface * s)
 	sym[c] = FALSE;
       n += (&m.x)[c];
     }
-    if (n > 0.) {
-      m.x /= n; m.y /= n; m.z /= n;
-      alpha = m.x*ca.x + m.y*ca.y + m.z*ca.z;
-      solid->a = gfs_plane_volume (&m, alpha, 1.);
-      gfs_plane_center (&m, alpha, solid->a, &solid->cm);
-      for (c = 0; c < FTT_DIMENSION; c++)
-	(&solid->cm.x)[c] = (&o.x)[c] + 
-	  (sym[c] ? 1. - (&solid->cm.x)[c] : (&solid->cm.x)[c])*(&h.x)[c];
-    }
-    else { /* degenerate intersections */
-      solid->a = 0.;
-      for (i = 0; i < FTT_NEIGHBORS; i++)
-	solid->a += solid->s[i];
-      solid->a /= FTT_NEIGHBORS;
-      if (solid->a == 0. || solid->a == 1.) {
-	g_free (solid);
-	GFS_STATE (cell)->solid = NULL;
-      }
+    g_assert (n > 0.);
+    m.x /= n; m.y /= n; m.z /= n;
+    alpha = m.x*ca.x + m.y*ca.y + m.z*ca.z;
+    solid->a = gfs_plane_volume (&m, alpha, 1.);
+    gfs_plane_center (&m, alpha, solid->a, &solid->cm);
+    for (c = 0; c < FTT_DIMENSION; c++)
+      (&solid->cm.x)[c] = (&o.x)[c] + 
+	(sym[c] ? 1. - (&solid->cm.x)[c] : (&solid->cm.x)[c])*(&h.x)[c];
+    break;
+  }
+  case 2:
+    solid->cm = solid->ca;
+    solid->a = 0.;
+    for (i = 0; i < FTT_NEIGHBORS; i++)
+      solid->a += solid->s[i];
+    solid->a /= FTT_NEIGHBORS;
+    if (solid->a == 0. || solid->a == 1.) {
+      g_free (solid);
+      GFS_STATE (cell)->solid = NULL;
     }
+    break;
+  default:
+    g_assert_not_reached ();
   }
 }
 

-- 
Gerris Flow Solver