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

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 6c73e25fee85e3c3b5cf814bb1f2113f003694bb
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Tue Sep 9 14:52:46 2008 +1000

    Axisymmetric surface tension implementation
    
    darcs-hash:20080909045246-d4795-cf6cabb34ca1180069dc0a4fcfb72f16b5b42554.gz

diff --git a/configure.in b/configure.in
index 246d521..bd7b862 100644
--- a/configure.in
+++ b/configure.in
@@ -14,7 +14,7 @@ dnl are available for $ac_help expansion (don't we all *love* autoconf?)
 # Remember to update rpm/gerris.spec when changing the version number.
 GFS_MAJOR_VERSION=1
 GFS_MINOR_VERSION=3
-GFS_MICRO_VERSION=0
+GFS_MICRO_VERSION=1
 GFS_INTERFACE_AGE=0
 GFS_BINARY_AGE=0
 GFS_VERSION=$GFS_MAJOR_VERSION.$GFS_MINOR_VERSION.$GFS_MICRO_VERSION
diff --git a/src/poisson.c b/src/poisson.c
index 2968cc6..d12e19a 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -472,11 +472,10 @@ void gfs_poisson_coefficients (GfsDomain * domain,
 static void tension_coeff (FttCellFace * face, gpointer * data)
 {
   gdouble * lambda2 = data[0];
-  GfsStateVector * s = GFS_STATE (face->cell);
   GfsSourceTensionGeneric * t = data[1];
-  gdouble v = lambda2[face->d/2]*t->sigma;
-  gdouble alpha = data[2] ? gfs_function_face_value (data[2], face) : 1.;
   GfsVariable * kappa = GFS_SOURCE_TENSION (data[1])->k;
+  gdouble alpha = data[2] ? gfs_function_face_value (data[2], face) : 1.;
+  gdouble v = lambda2[face->d/2]*alpha*gfs_domain_face_fraction (kappa->domain, face)*t->sigma;
   gdouble k1 = GFS_VARIABLE (face->cell, kappa->i);
   gdouble k2 = GFS_VARIABLE (face->neighbor, kappa->i);
 #if 0
@@ -511,11 +510,7 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
 	   "Please check your definition.",
 	   alpha, p.x, p.y, p.z);
   }
-  v *= alpha;
-  /* fixme: mapping */
-  if (GFS_IS_MIXED (face->cell))
-    v *= s->solid->s[face->d];
-  s->f[face->d].v = v;
+  GFS_STATE (face->cell)->f[face->d].v = v;
 
   switch (ftt_face_type (face)) {
   case FTT_FINE_FINE:
diff --git a/src/tension.c b/src/tension.c
index c675498..6196eb9 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -640,6 +640,9 @@ static void variable_curvature_from_distance (GfsEvent * event, GfsSimulation *
 
   gfs_domain_timer_start (domain, "variable_curvature");
 
+  if (GFS_IS_AXI (sim))
+    g_assert_not_implemented ();
+
   for (c = 0; c < FTT_DIMENSION + 1; c++) {
     n[c] = gfs_temporary_variable (domain);
     gfs_variable_set_vector (n[c], c);
diff --git a/src/vof.c b/src/vof.c
index d48b5cb..b33a3ef 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -1877,6 +1877,19 @@ static gboolean curvature_along_direction (FttCell * cell,
     *kappa = hxx/(size*sqrt (dnm*dnm*dnm));
     if (kmax)
       *kmax = fabs (*kappa);
+    if (GFS_IS_AXI (v->domain)) {
+      gdouble nr, r = p.y;
+      if (c == FTT_X)
+	nr = hx;
+      else {
+	r += (d == FTT_TOP ? - H*size : H*size);
+	nr = (d == FTT_TOP ? 1. : -1.);
+      }
+      gdouble kaxi = nr/(sqrt(dnm)*r);
+      *kappa += kaxi;
+      if (kmax)
+	*kmax = MAX (*kmax, fabs (kaxi));
+    }
   }
 #else  /* 3D */  
   static FttComponent or[3][2] = { { FTT_Y, FTT_Z }, { FTT_X, FTT_Z }, { FTT_X, FTT_Y } };
@@ -2113,6 +2126,18 @@ static gdouble parabola_fit_curvature (ParabolaFit * p, gdouble kappamax,
   return kappa;
 }
 
+#if FTT_2D
+static void parabola_fit_axi_curvature (const ParabolaFit * p, gdouble r, 
+					gdouble * kappa, gdouble * kmax)
+{
+  gdouble nr = (p->m[0]*p->a[1] + p->m[1])/sqrt (1. + p->a[1]*p->a[1]);
+  gdouble kaxi = - nr/r;
+  *kappa += kaxi;
+  if (kmax)
+    *kmax = MAX (*kmax, fabs (kaxi));
+}
+#endif /* 2D */
+
 static void parabola_fit_destroy (ParabolaFit * p)
 {
 #if (FTT_2D || PARABOLA_SIMPLER)
@@ -2207,6 +2232,10 @@ gdouble gfs_fit_curvature (FttCell * cell, GfsVariableTracerVOF * t, gdouble * k
   gdouble kappa = parabola_fit_curvature (&fit, 2., kmax)/h;
   if (kmax)
     *kmax /= h;
+#if FTT_2D
+  if (GFS_IS_AXI (v->domain))
+    parabola_fit_axi_curvature (&fit, fc.y*h + p.y, &kappa, kmax);
+#endif
   parabola_fit_destroy (&fit);
   return kappa;
 }
@@ -2319,6 +2348,10 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t, gdouble
   kappa = parabola_fit_curvature (&fit, 2., kmax)/h;
   if (kmax)
     *kmax /= h;
+#if FTT_2D
+  if (GFS_IS_AXI (v->domain))
+    parabola_fit_axi_curvature (&fit, fc.y*h + p.y, &kappa, kmax);
+#endif
   parabola_fit_destroy (&fit);
   return kappa;
 }

-- 
Gerris Flow Solver