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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit a1e1996231236f3c2b8dda9e9cf6af9268af9752
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Fri Jun 20 08:44:41 2008 +1000

    More accurate streamlines
    
    Uses new function gfs_mixed_cell_interpolate().
    
    darcs-hash:20080619224441-d4795-76cd5d88dfaf018378df7df29516277af53c9f4b.gz

diff --git a/src/fluid.c b/src/fluid.c
index 7030890..cc18d8f 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -1302,6 +1302,59 @@ void gfs_mixed_cell_gradient (FttCell * cell,
 }
 
 /**
+ * gfs_mixed_cell_interpolate:
+ * @cell: a mixed #FttCell.
+ * @p: a #FttVector.
+ * @v: a #GfsVariable.
+ *
+ * Returns: the value of variable @v interpolated at position @p within @cell.
+ */
+gdouble gfs_mixed_cell_interpolate (FttCell * cell,
+				    FttVector p,
+				    GfsVariable * v)
+{
+  FttCell * n[N_CELLS];
+  gdouble m[N_CELLS - 1][N_CELLS - 1], a[N_CELLS - 1];
+  gdouble v0, h;
+  FttVector * o;
+  guint i, j;
+
+  g_return_val_if_fail (cell != NULL, 0.);
+  g_return_val_if_fail (GFS_IS_MIXED (cell), 0.);
+  g_return_val_if_fail (v != NULL, 0.);
+
+  o = &GFS_STATE (cell)->solid->cm;
+  v0 = GFS_VALUE (cell, v);
+
+  if (v->surface_bc) {
+    (* GFS_SURFACE_GENERIC_BC_CLASS (GTS_OBJECT (v->surface_bc)->klass)->bc) (cell, v->surface_bc);
+    if (((cell)->flags & GFS_FLAG_DIRICHLET) != 0) {
+      o = &GFS_STATE (cell)->solid->ca;
+      v0 = GFS_STATE (cell)->solid->fv;
+    }
+  }
+  g_assert (cell_bilinear (cell, n, o, gfs_cell_cm, -1, m));
+
+  for (i = 0; i < N_CELLS - 1; i++) {
+    a[i] = 0.;
+    for (j = 0; j < N_CELLS - 1; j++)
+      a[i] += m[i][j]*(GFS_VALUE (n[j + 1], v) - v0);
+  }
+  
+  h = ftt_cell_size (cell);
+  p.x = (p.x - o->x)/h;
+  p.y = (p.y - o->y)/h;
+#if FTT_2D
+  return a[0]*p.x + a[1]*p.y + a[2]*p.x*p.y + v0;
+#else /* 3D */
+  p.z = (p.z - o->z)/h;
+  return (a[0]*p.x + a[1]*p.y + a[2]*p.z + 
+	  a[3]*p.x*p.y + a[4]*p.x*p.z + a[5]*p.y*p.z + 
+	  a[6]*p.x*p.y*p.z + v0);
+#endif /* 3D */
+}
+
+/**
  * gfs_cell_dirichlet_gradient_flux:
  * @cell: a #FttCell.
  * @v: a #GfsVariable index.
diff --git a/src/fluid.h b/src/fluid.h
index 75998f5..44d9c4f 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -152,6 +152,9 @@ gdouble               gfs_cell_dirichlet_gradient_flux (FttCell * cell,
 gdouble               gfs_cell_dirichlet_value         (FttCell * cell,
 							GfsVariable * v,
 							gint max_level);
+gdouble               gfs_mixed_cell_interpolate       (FttCell * cell,
+							FttVector p,
+							GfsVariable * v);
 
 void                  gfs_normal_divergence          (FttCell * cell,
 						      GfsVariable * v);
diff --git a/src/graphic.c b/src/graphic.c
index 4522dec..8859263 100644
--- a/src/graphic.c
+++ b/src/graphic.c
@@ -1668,6 +1668,21 @@ static void vorticity_vector (FttCell * cell, gpointer * data)
 }
 #endif /* 3D */
 
+static gdouble interpolated_velocity (FttCell * cell, FttVector p, GfsVariable ** U,
+				      gdouble direction,
+				      FttVector * u)
+{
+  FttComponent c;
+  gdouble nu = 0.;
+  gdouble (* interpolate) (FttCell *, FttVector, GfsVariable * v) = GFS_IS_MIXED (cell) ?
+    gfs_mixed_cell_interpolate : gfs_interpolate;
+  for (c = 0; c < FTT_DIMENSION; c++) {
+    (&u->x)[c] = direction* (*interpolate) (cell, p, U[c]);
+    nu += (&u->x)[c]*(&u->x)[c];
+  }
+  return nu;
+}
+
 static GList * grow_curve (GfsDomain * domain,
 			   GfsVariable ** U,
 			   FttVector p,
@@ -1745,10 +1760,7 @@ static GList * grow_curve (GfsDomain * domain,
       nstep = 0;
     }
 
-    for (c = 0; c < FTT_DIMENSION; c++) {
-      (&u.x)[c] = direction*gfs_interpolate (cell, p, U[c]);
-      nu += (&u.x)[c]*(&u.x)[c];
-    }
+    nu = interpolated_velocity (cell, p, U, direction, &u);
     if (nu > 0) {
       FttVector p1 = p;
       FttCell * cell1;
@@ -1756,14 +1768,10 @@ static GList * grow_curve (GfsDomain * domain,
       nu = 2.*sqrt (nu);
       for (c = 0; c < FTT_DIMENSION; c++)
 	(&p1.x)[c] += h*(&u.x)[c]/nu;
-      nu = 0.;
       cell1 = gfs_domain_locate (domain, p1, -1);
       if (!cell1)
 	break;
-      for (c = 0; c < FTT_DIMENSION; c++) {
-	(&u.x)[c] = direction*gfs_interpolate (cell1, p1, U[c]);
-	nu += (&u.x)[c]*(&u.x)[c];
-      }
+      nu = interpolated_velocity (cell1, p1, U, direction, &u);
     }
     else
       break;

-- 
Gerris Flow Solver