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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit db21a034dc3df8f1e517d894f9aa8a434c55d66e
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu Dec 14 12:30:56 2006 +1100

    Density (and viscosity) is computed using face-interpolated volume fraction values
    
    This uses the new gfs_function_face_value(). Large density ratio flows
    can now be computed correctly. Needs further testing.
    
    darcs-hash:20061214013056-d4795-7a091ccfa81f6b5f2e774585845f4d6e82f22ef7.gz

diff --git a/src/poisson.c b/src/poisson.c
index 9139b89..3a5c862 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -331,29 +331,17 @@ static void poisson_coeff (FttCellFace * face, gdouble * lambda2)
   }
 }
 
-static void reset_alpha_coeff (FttCell * cell, gpointer * data)
-{
-  FttDirection d;
-  GfsFaceStateVector * f = GFS_STATE (cell)->f;
-  GfsFunction * alpha = data[0];
-  GfsVariable * a = data[1];
-  
-  for (d = 0; d < FTT_NEIGHBORS; d++)
-    f[d].v = 0.;
-  GFS_VARIABLE (cell, a->i) = gfs_function_value (alpha, cell);
-}
-
 static void poisson_alpha_coeff (FttCellFace * face,
 				 gpointer * data)
 {
   gdouble * lambda2 = data[0];
-  GfsVariable * alpha = data[1];
+  GfsFunction * alpha = data[1];
   gdouble v = lambda2[face->d/2];
   GfsStateVector * s = GFS_STATE (face->cell);
 
   if (GFS_IS_MIXED (face->cell))
     v *= s->solid->s[face->d];
-  v *= gfs_face_interpolated_value (face, alpha->i);
+  v *= gfs_function_face_value (alpha, face);
   s->f[face->d].v = v;
 
   switch (ftt_face_type (face)) {
@@ -418,29 +406,21 @@ void gfs_poisson_coefficients (GfsDomain * domain,
 
     lambda2[i] = lambda*lambda;
   }
-  if (alpha == NULL) {
-    gfs_domain_cell_traverse (domain,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_coeff, NULL);
+  gfs_domain_cell_traverse (domain,
+			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			    (FttCellTraverseFunc) reset_coeff, NULL);
+  if (alpha == NULL)
     gfs_domain_face_traverse (domain, FTT_XYZ, 
 			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			      (FttFaceTraverseFunc) poisson_coeff, lambda2);
-  }
   else {
     gpointer data[2];
 
-    data[0] = alpha;
-    data[1] = gfs_temporary_variable (domain);
-    gfs_domain_cell_traverse (domain,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_alpha_coeff, data);
-    gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, data[1]);
     data[0] = lambda2;
+    data[1] = alpha;
     gfs_domain_face_traverse (domain, FTT_XYZ, 
 			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttFaceTraverseFunc) poisson_alpha_coeff, 
-			      data);
-    gts_object_destroy (data[1]);
+			      (FttFaceTraverseFunc) poisson_alpha_coeff, data);
   }
   gfs_domain_cell_traverse (domain,
 			    FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
@@ -453,7 +433,8 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
   GfsStateVector * s = GFS_STATE (face->cell);
   GfsSourceTensionGeneric * t = data[1];
   gdouble v = lambda2[face->d/2]*t->sigma;
-  GfsVariable * alpha = data[2], * kappa = GFS_SOURCE_TENSION (data[1])->k;
+  GfsFunction * alpha = data[2];
+  GfsVariable * kappa = GFS_SOURCE_TENSION (data[1])->k;
   gdouble c1 = GFS_VARIABLE (face->cell, t->c->i);
   gdouble c2 = GFS_VARIABLE (face->neighbor, t->c->i);
   gdouble w1 = c1*(1. - c1);
@@ -477,7 +458,7 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
   }
 
   if (alpha)
-    v *= gfs_face_interpolated_value (face, alpha->i);
+    v *= gfs_function_face_value (alpha, face);
   if (GFS_IS_MIXED (face->cell))
     v *= s->solid->s[face->d];
   s->f[face->d].v = v;
@@ -521,28 +502,15 @@ void gfs_source_tension_coefficients (GfsSourceTension * s,
 
     lambda2[i] = lambda*lambda;
   }
-  if (alpha) {
-    data[0] = alpha;
-    data[1] = gfs_temporary_variable (domain);
-    gfs_domain_cell_traverse (domain,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_alpha_coeff, data);
-    gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, data[1]);
-    data[2] = data[1];
-  }
-  else {
-    gfs_domain_cell_traverse (domain,
-			      FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
-			      (FttCellTraverseFunc) reset_coeff, NULL);
-    data[2] = NULL;
-  }
+  gfs_domain_cell_traverse (domain,
+			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+			    (FttCellTraverseFunc) reset_coeff, NULL);
   data[0] = lambda2;
   data[1] = s;
+  data[2] = alpha;
   gfs_domain_face_traverse (domain, FTT_XYZ, 
 			    FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
 			    (FttFaceTraverseFunc) tension_coeff, data);
-  if (alpha)
-    gts_object_destroy (data[2]);
 }
 
 static void correct (FttCell * cell, gpointer * data)
diff --git a/src/source.c b/src/source.c
index 9cb28f4..23cc49c 100644
--- a/src/source.c
+++ b/src/source.c
@@ -462,8 +462,7 @@ static gboolean diffusion_event (GfsEvent * event, GfsSimulation * sim)
 
 static gdouble diffusion_face (GfsDiffusion * d, FttCellFace * f)
 {
-  return d->mu ? gfs_face_interpolated_value (f, d->mu->i) :
-    gfs_function_get_constant_value (d->val);
+  return gfs_function_face_value (d->val, f);
 }
 
 static gdouble diffusion_cell (GfsDiffusion * d, FttCell * cell)

-- 
Gerris Flow Solver