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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 1d88bced2f6346b0348878cb8ac6d9a920e013da
Author: Stephane Popinet <popinet at users.sourceforge.net>
Date:   Mon Nov 1 10:08:09 2004 +1100

    Added VOF calculation of 3D solid fraction (gerris--mainline--0.7--patch-17)
    
    gerris--mainline--0.7--patch-17
    Keywords:
    
    The center of mass calculation is not done yet.
    
    darcs-hash:20041031230809-aabb8-186b3c9e339eb7fa4b4d0ec76a44168825c3bdfb.gz

diff --git a/src/isocube.h b/src/isocube.h
index e40304b..0093dec 100644
--- a/src/isocube.h
+++ b/src/isocube.h
@@ -28,23 +28,6 @@ static FttVector vertex[8] = {
   {0.,0.,0.},{0.,0.,1.},{0.,1.,0.},{0.,1.,1.},
   {1.,0.,0.},{1.,0.,1.},{1.,1.,0.},{1.,1.,1.}
 };
-/* first index is the edge number, second index is the edge orientation 
-   (0 or 1), third index are the edges which this edge may connect to
-   in order and the corresponding face direction */
-static guint connect[12][2][4] = {
-  {{9, 1, 8, FTT_BOTTOM}, {4, 3, 7, FTT_BACK}},   /* 0 */
-  {{6, 2, 5, FTT_FRONT},  {8, 0, 9, FTT_BOTTOM}}, /* 1 */
-  {{10, 3, 11, FTT_TOP},  {5, 1, 6, FTT_FRONT}},  /* 2 */
-  {{7, 0, 4, FTT_BACK},   {11, 2, 10, FTT_TOP}},  /* 3 */
-  {{3, 7, 0, FTT_BACK},   {8, 5, 11, FTT_LEFT}},  /* 4 */
-  {{11, 4, 8, FTT_LEFT},  {1, 6, 2, FTT_FRONT}},  /* 5 */
-  {{2, 5, 1, FTT_FRONT},  {9, 7, 10, FTT_RIGHT}}, /* 6 */
-  {{10, 6, 9, FTT_RIGHT}, {0, 4, 3, FTT_BACK}},   /* 7 */
-  {{5, 11, 4, FTT_LEFT},  {0, 9, 1, FTT_BOTTOM}}, /* 8 */
-  {{1, 8, 0, FTT_BOTTOM}, {7, 10, 6, FTT_RIGHT}}, /* 9 */
-  {{6, 9, 7, FTT_RIGHT},  {3, 11, 2, FTT_TOP}},   /* 10 */
-  {{2, 10, 3, FTT_TOP},   {4, 8, 5, FTT_LEFT}}    /* 11 */
-};
 static guint face[6][4][2] = {
   {{7,0},{10,0},{6,1},{9,1}}, /* right */
   {{4,0},{11,0},{5,1},{8,1}}, /* left */
diff --git a/src/solid.c b/src/solid.c
index 3680130..60463a3 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -278,7 +278,7 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
   GfsSolidVector * solid = GFS_STATE (cell)->solid;
   gdouble h = ftt_cell_size (cell);
   CellCube cube;
-  FttVector o;
+  FttVector o, ca = {0., 0., 0.};
   guint i, n1 = 0;
   gint inside[8] = {0,0,0,0,0,0,0,0};
 
@@ -300,7 +300,13 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
       guint j = edge1[i][0], k = edge1[i][1];
 
       /* intersection vertex position is the average of all the n[i] intersections */
-      cube.x[i] /= cube.n[i]; 
+      cube.x[i] /= cube.n[i];
+
+      /* average of all intersections */
+      ca.x += (1. - cube.x[i])*cube.p[j].x + cube.x[i]*cube.p[k].x;
+      ca.y += (1. - cube.x[i])*cube.p[j].y + cube.x[i]*cube.p[k].y;
+      ca.z += (1. - cube.x[i])*cube.p[j].z + cube.x[i]*cube.p[k].z;
+
       g_assert (inside[j] == 0 || inside[j] == cube.inside[i]);
       g_assert (inside[k] == 0 || inside[k] == - cube.inside[i]);
       inside[j] = cube.inside[i];
@@ -320,18 +326,19 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
 
   if (!solid)
     GFS_STATE (cell)->solid = solid = g_malloc0 (sizeof (GfsSolidVector));
-  
-  for (i = 0; i < FTT_NEIGHBORS; i++) { /* for each face */
+
+  /* compute face fractions */
+  for (i = 0; i < FTT_NEIGHBORS; i++) {
     CellFace f;
-    guint j;
+    guint j, n2;
 
-    n1 = 0;
+    n2 = 0;
     for (j = 0; j < 4; j++) { /* initialise face i */
       guint e = face[i][j][0];
 
       f.p[j] = cube.p[face_v[i][j]];
       f.n[j] = cube.n[e];
-      if (f.n[j]) n1++;
+      if (f.n[j]) n2++;
       if (face[i][j][1]) {
 	f.x[j] = 1. - cube.x[e];
 	f.inside[j] = - cube.inside[e];
@@ -342,7 +349,7 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
       }
     }
 
-    switch (n1) {
+    switch (n2) {
     case 0: { /* the face is not cut */
       gint ins = 0;
 
@@ -367,8 +374,31 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
     }
     default:
       g_log (G_LOG_DOMAIN, G_LOG_LEVEL_ERROR,
-	     "the solid surface is probably not closed (n1 = %d)", n1);
+	     "the solid surface is probably not closed (n2 = %d)", n2);
+    }
+  }
+
+  /* now compute cell fraction, center of area, center of mass */
+  ca.x /= n1; ca.y /= n1; ca.z /= n1;
+  solid->ca = ca;
+  {
+    FttVector m;
+    gdouble alpha, n = 0.;
+    FttComponent c;
+
+    for (c = 0; c < FTT_DIMENSION; c++) {
+      (&ca.x)[c] = ((&ca.x)[c] - (&o.x)[c])/h;
+      (&m.x)[c] = solid->s[2*c + 1] - solid->s[2*c];
+      if ((&m.x)[c] < 0.) {
+	(&m.x)[c] = - (&m.x)[c];
+	(&ca.x)[c] = 1. - (&ca.x)[c];
+      }
+      (&m.x)[c] += 1e-6;
+      n += (&m.x)[c];
     }
+    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.);
   }
 }
 #endif /* 3D */

-- 
Gerris Flow Solver