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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 17218e8605aadc180433321c489278ddd6ce0e0a
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Tue Feb 21 03:18:27 2006 +1100

    Simplification of gfs_plane_volume()
    
    darcs-hash:20060220161827-d4795-1009c06583db1ad422a838c26a9cc5d4da362039.gz

diff --git a/src/solid.c b/src/solid.c
index 6ba524a..1f98a0f 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -627,18 +627,25 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s, In
 	sym[c] = FALSE;
       n += (&m.x)[c];
     }
-    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];
-  }
-  else { /* this is a "thin" cell */
-    p->thin++;
-    deal_with_thin_cell (cell, p);
+    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;
+      }
+    }
   }
   if (solid->a == 0.)
     GFS_VARIABLE (cell, p->status->i) = 1.;
diff --git a/src/vof.c b/src/vof.c
index c8c6113..d8d53b5 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -27,12 +27,11 @@
  * gfs_line_area:
  * @m: normal to the line.
  * @alpha: line constant.
- * @c1: width of the cell.
  *
- * Returns: the area of the fraction of a rectangular cell (c1,1)
- * lying under the line (@m, at alpha).
+ * Returns: the area of the fraction of a cell lying under the line
+ * (@m, at alpha).
  */
-gdouble gfs_line_area (FttVector * m, gdouble alpha, gdouble c1)
+gdouble gfs_line_area (FttVector * m, gdouble alpha)
 {
   FttVector n;
   gdouble a, v;
@@ -43,14 +42,14 @@ gdouble gfs_line_area (FttVector * m, gdouble alpha, gdouble c1)
   if (alpha <= 0.)
     return 0.;
 
-  if (alpha >= m->x*c1 + m->y || c1 == 0.)
-    return c1;
+  if (alpha >= m->x + m->y)
+    return 1.;
 
   n = *m; n.x += 1e-4; n.y += 1e-4;
 
   v = alpha*alpha;
 
-  a = alpha - n.x*c1;
+  a = alpha - n.x;
   if (a > 0.)
     v -= a*a;
 
@@ -117,12 +116,10 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
  * gfs_plane_volume:
  * @m: normal to the plane.
  * @alpha: plane constant.
- * @c1: width of the cell.
  *
- * Returns: the volume of a parallelepipedic cell (c1,1,1) lying under
- * the plane (@m, at alpha).
+ * Returns: the volume of a cell lying under the plane (@m, at alpha).
  */
-gdouble gfs_plane_volume (FttVector * m, gdouble alpha, gdouble c1)
+gdouble gfs_plane_volume (FttVector * m, gdouble alpha)
 {
   FttVector n;
   gdouble a, amax, v;
@@ -135,29 +132,23 @@ gdouble gfs_plane_volume (FttVector * m, gdouble alpha, gdouble c1)
   if (alpha <= 0.)
     return 0.;
 
-  if (alpha >= m->x*c1 + m->y + m->z || c1 == 0.)
-    return c1;
+  if (alpha >= m->x + m->y + m->z)
+    return 1.;
 
   n = *m; n.x += 1e-4; n.y += 1e-4; n.z += 1e-4;
-  amax = n.x*c1 + n.y + n.z;
+  amax = n.x + n.y + n.z;
 
   md = &n.x;
   v = alpha*alpha*alpha;
 
-  a = alpha - n.x*c1;
-  if (a > 0.)
-    v -= a*a*a;
-  for (j = 1; j < 3; j++) {
+  for (j = 0; j < 3; j++) {
     a = alpha - md[j];
     if (a > 0.)
       v -= a*a*a;
   }
 
   amax = alpha - amax;
-  a = amax + n.x*c1;
-  if (a > 0.)
-    v += a*a*a;
-  for (j = 1; j < 3; j++) {
+  for (j = 0; j < 3; j++) {
     a = amax + md[j];
     if (a > 0.)
       v += a*a*a;
@@ -470,15 +461,19 @@ static void gfs_cell_vof_advected_face_values (FttCell * cell,
     m.x /= 1. + u_right - u_left;
     alpha += m.x*u_left;
 
-    if (u_left < 0.)
-      GFS_STATE (cell)->f[left].v =
-	- gfs_plane_volume (&m, alpha - m.x*u_left, - u_left)/u_left;
+    if (u_left < 0.) {
+      m1 = m;
+      m1.x *= - u_left;
+      GFS_STATE (cell)->f[left].v = gfs_plane_volume (&m1, alpha + m1.x);
+    }
     else
       GFS_STATE (cell)->f[left].v = f;
 
-    if (u_right > 0.)
-      GFS_STATE (cell)->f[right].v =
-	gfs_plane_volume (&m, alpha - m.x, u_right)/u_right;
+    if (u_right > 0.) {
+      m1 = m;
+      m1.x *= u_right;
+      GFS_STATE (cell)->f[right].v = gfs_plane_volume (&m1, alpha - m.x);
+    }
     else
       GFS_STATE (cell)->f[right].v = f;
   }
@@ -667,7 +662,7 @@ void gfs_vof_coarse_fine (FttCell * parent, GfsVariable * v)
 	(&p.x)[c] = (&m.x)[c] < 0. ? (&p.x)[c] + 0.25 : 0.25 - (&p.x)[c];
 	alpha1 -= (&m1.x)[c]*(&p.x)[c];
       }
-      GFS_VARIABLE (child.c[i], v->i) = gfs_plane_volume (&m1, 2.*alpha1, 1.);
+      GFS_VARIABLE (child.c[i], v->i) = gfs_plane_volume (&m1, 2.*alpha1);
     }
   }
 }
diff --git a/src/vof.h b/src/vof.h
index abe952e..9c8a317 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -29,8 +29,7 @@ extern "C" {
 #define GFS_IS_FULL(f)             ((f) < 1e-6 || (f) > 1. - 1.e-6)
 
 gdouble gfs_line_area              (FttVector * m, 
-				    gdouble alpha, 
-				    gdouble c1);
+				    gdouble alpha);
 void    gfs_line_center            (FttVector * m, 
 				    gdouble alpha, 
 				    gdouble a, 
@@ -42,8 +41,7 @@ gdouble gfs_line_alpha             (FttVector * m,
 #  define gfs_plane_alpha          gfs_line_alpha
 #else /* 3D */
 gdouble gfs_plane_volume           (FttVector * m, 
-				    gdouble alpha, 
-				    gdouble c1);
+				    gdouble alpha);
 gdouble gfs_plane_alpha            (FttVector * m, 
 				    gdouble c);
 void    gfs_plane_center           (FttVector * m, 

-- 
Gerris Flow Solver