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

Stephane Popinet popinet at users.sf.net
Fri May 15 02:54:27 UTC 2009


The following commit has been merged in the upstream branch:
commit f2dc10aa5a62af7aeff759cb68f3ce06d16ad6e9
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Thu May 10 09:59:46 2007 +1000

    Paraboloid fitting for curvature calculation in 3D
    
    darcs-hash:20070509235946-d4795-e799265378bc55564ad96eaf4fc2fc7c31f51d87.gz

diff --git a/src/solid.c b/src/solid.c
index 759cf3e..01fb64e 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -1340,9 +1340,11 @@ void gfs_solid_coarse_fine (FttCell * parent)
 
       ftt_cell_pos (child.c[i], &p);
       gfs_plane_center (&m, 2.*alpha1, a, &s->cm);
-      for (c = 0; c < FTT_DIMENSION; c++)
+      gfs_plane_area_center (&m, 2.*alpha1, &s->ca);
+      for (c = 0; c < FTT_DIMENSION; c++) {
 	(&s->cm.x)[c] = (&p.x)[c] + h*((&s->cm.x)[c] - 0.5);
-      g_assert (gfs_vof_plane_center (child.c[i], &m, 2.*alpha1, &s->ca));
+	(&s->ca.x)[c] = (&p.x)[c] + h*((&s->ca.x)[c] - 0.5);
+      }
 
       FttDirection d;
       FttCellNeighbors n;
diff --git a/src/tension.c b/src/tension.c
index 1d28cb2..6086278 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -383,7 +383,7 @@ static void curvature (FttCell * cell, GfsVariable * v)
   gdouble f = GFS_VARIABLE (cell, t->i);
 
   if (GFS_IS_FULL (f))
-    GFS_VARIABLE (cell, v->i) = 0.;
+    GFS_VARIABLE (cell, v->i) = G_MAXDOUBLE;
   else
     GFS_VARIABLE (cell, v->i) = gfs_height_curvature (cell, GFS_VARIABLE_TRACER_VOF (t));
 }
@@ -445,6 +445,7 @@ static void variable_curvature_from_fraction (GfsEvent * event, GfsSimulation *
   gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
 			    (FttCellTraverseFunc) GFS_VARIABLE1 (event)->fine_coarse, event);
   gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, GFS_VARIABLE1 (event));
+  variable_curvature_diffuse (event, sim);
 
   gfs_domain_timer_stop (domain, "variable_curvature");
 }
diff --git a/src/vof.c b/src/vof.c
index a43a839..b3ebf3c 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -18,6 +18,7 @@
  */
 
 #include <math.h>
+#include <stdlib.h>
 #include "vof.h"
 #include "variable.h"
 #include "graphic.h"
@@ -177,6 +178,68 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
     p->y = 1. - p->y;
 }
 
+/**
+ * gfs_line_area_center:
+ * @m: normal to the line.
+ * @alpha: line constant.
+ * @p: a #FttVector.
+ *
+ * Fills @p with the position of the center of area of the fraction of
+ * a square cell lying under the line (@m, at alpha).
+ */
+void gfs_line_area_center (FttVector * m, gdouble alpha, FttVector * p)
+{
+  FttVector n;
+
+  g_return_if_fail (m != NULL);
+  g_return_if_fail (p != NULL);
+
+  n = *m;
+  if (n.x < 0.) {
+    alpha -= n.x;
+    n.x = - n.x;
+  }
+  if (n.y < 0.) {
+    alpha -= n.y;
+    n.y = - n.y;
+  }
+
+  p->z = 0.;
+  if (alpha <= 0. || alpha >= n.x + n.y) {
+    p->x = p->y = 0.;
+    return;
+  }
+
+  n.x += 1e-4; n.y += 1e-4;
+
+  p->x = p->y = 0.;
+
+  if (alpha >= n.x) {
+    p->x += 1.;
+    p->y += (alpha - n.x)/n.y;
+  }
+  else
+    p->x += alpha/n.x;
+
+  if (alpha >= n.y) {
+    p->y += 1.;
+    p->x += (alpha - n.y)/n.x;
+  }
+  else
+    p->y += alpha/n.y;
+  
+  p->x /= 2.;
+  p->y /= 2.;
+
+  g_assert (p->x >= 0. && p->x <= 1.);
+  g_assert (p->y >= 0. && p->y <= 1.);
+
+  if (m->x < 0.)
+    p->x = 1. - p->x;
+  if (m->y < 0.)
+    p->y = 1. - p->y;
+}
+
 #if (!FTT_2D)
 /**
  * gfs_plane_volume:
@@ -408,6 +471,128 @@ void gfs_plane_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
   if (m->y < 0.) p->y = 1. - p->y;
   if (m->z < 0.) p->z = 1. - p->z;
 }
+
+/**
+ * gfs_plane_area_center:
+ * @m: normal to the plane.
+ * @alpha: plane constant.
+ * @p: a #FttVector.
+ *
+ * Fills @p with the position of the center of area of the fraction of
+ * a cubic cell lying under the plane (@m, at alpha).
+ */
+void gfs_plane_area_center (FttVector * m, gdouble alpha, FttVector * p)
+{
+  g_return_if_fail (m != NULL);
+  g_return_if_fail (p != NULL);
+
+  if (fabs (m->x) < 1e-4) {
+    FttVector n, q;
+    n.x = m->y;
+    n.y = m->z;
+    gfs_line_area_center (&n, alpha, &q);
+    p->x = 0.5;
+    p->y = q.x;
+    p->z = q.y;
+    return;
+  }
+  if (fabs (m->y) < 1e-4) {
+    FttVector n, q;
+    n.x = m->z;
+    n.y = m->x;
+    gfs_line_area_center (&n, alpha, &q);
+    p->x = q.y;
+    p->y = 0.5;
+    p->z = q.x;
+    return;
+  }
+  if (fabs (m->z) < 1e-4) {
+    gfs_line_area_center (m, alpha, p);
+    p->z = 0.5;
+    return;
+  }
+
+  FttVector n = *m;
+  if (n.x < 0.) {
+    alpha -= n.x;
+    n.x = - n.x;
+  }
+  if (n.y < 0.) {
+    alpha -= n.y;
+    n.y = - n.y;
+  }  
+  if (n.z < 0.) {
+    alpha -= n.z;
+    n.z = - n.z;
+  }
+
+  gdouble amax = n.x + n.y + n.z;
+  if (alpha <= 0. || alpha >= amax) {
+    p->x = p->y = p->z = 0.;
+    return;
+  }
+
+  gdouble area = alpha*alpha;
+  p->x = p->y = p->z = area*alpha;
+
+  gdouble b = alpha - n.x;
+  if (b > 0.) {
+    area -= b*b;
+    p->x -= b*b*(2.*n.x + alpha);
+    p->y -= b*b*b;
+    p->z -= b*b*b;
+  }
+  b = alpha - n.y;
+  if (b > 0.) {
+    area -= b*b;
+    p->y -= b*b*(2.*n.y + alpha);
+    p->x -= b*b*b;
+    p->z -= b*b*b;
+  }
+  b = alpha - n.z;
+  if (b > 0.) {
+    area -= b*b;
+    p->z -= b*b*(2.*n.z + alpha);
+    p->x -= b*b*b;
+    p->y -= b*b*b;
+  }
+
+  amax = alpha - amax;
+  b = amax + n.x;
+  if (b > 0.) {
+    area += b*b;
+    p->y += b*b*(2.*n.y + alpha - n.z);
+    p->z += b*b*(2.*n.z + alpha - n.y);
+    p->x += b*b*b;
+  }
+  b = amax + n.y;
+  if (b > 0.) {
+    area += b*b;
+    p->x += b*b*(2.*n.x + alpha - n.z);
+    p->z += b*b*(2.*n.z + alpha - n.x);
+    p->y += b*b*b;
+  }
+  b = amax + n.z;
+  if (b > 0.) {
+    area += b*b;
+    p->x += b*b*(2.*n.x + alpha - n.y);
+    p->y += b*b*(2.*n.y + alpha - n.x);
+    p->z += b*b*b;
+  }
+
+  area *= 3.;
+  p->x /= area*n.x;
+  p->y /= area*n.y;
+  p->z /= area*n.z;
+
+  g_assert (p->x >= 0. && p->x <= 1.);
+  g_assert (p->y >= 0. && p->y <= 1.);
+  g_assert (p->z >= 0. && p->z <= 1.);
+
+  if (m->x < 0.) p->x = 1. - p->x;
+  if (m->y < 0.) p->y = 1. - p->y;
+  if (m->z < 0.) p->z = 1. - p->z;
+}
 #endif /* 3D */
 
 /**
@@ -452,6 +637,54 @@ void gfs_youngs_gradient (FttCell * cell, GfsVariable * v, FttVector * g)
 }
 
 /**
+ * gfs_vof_plane_interpolate:
+ * @cell: a #FttCell containing location @p.
+ * @p: the center of the virtual cell.
+ * @level: the level of the virtual cell.
+ * @t: a #GfsVariableTracerVOF.
+ * @m: a #FttVector.
+ *
+ * Computes the equation @m.x = alpha of the volume fraction plane of
+ * a virtual cell at @level centered on @p.
+ *
+ * Returns: alpha for the virtual cell.
+ */
+gdouble gfs_vof_plane_interpolate (FttCell * cell,
+				   FttVector * p,
+				   guint level,
+				   GfsVariableTracerVOF * t,
+				   FttVector * m)
+{
+  guint l = ftt_cell_level (cell);
+
+  g_return_val_if_fail (cell != NULL, 0.);
+  g_return_val_if_fail (l <= level, 0.);
+  g_return_val_if_fail (t != NULL, 0.);
+  g_return_val_if_fail (m != NULL, 0.);
+
+  GfsVariable * v = GFS_VARIABLE1 (t);
+  gdouble f = GFS_VARIABLE (cell, v->i);
+  g_return_val_if_fail (!GFS_IS_FULL (f), 0.);
+  FttComponent c;
+  for (c = 0; c < FTT_DIMENSION; c++)
+    (&m->x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
+
+  gdouble alpha = GFS_VARIABLE (cell, t->alpha->i);
+  if (l < level) {
+    gdouble h = ftt_level_size (level);
+    gdouble H = ftt_cell_size (cell);
+    FttVector q;
+    
+    ftt_cell_pos (cell, &q);
+    alpha *= H;
+    for (c = 0; c < FTT_DIMENSION; c++)
+      alpha -= (&m->x)[c]*((&p->x)[c] - h/2. - (&q.x)[c] + H/2);
+    alpha /= h;
+  }
+  return alpha;
+}
+
+/**
  * gfs_vof_interpolate:
  * @cell: a #FttCell containing location @p.
  * @p: the center of the virtual cell.
@@ -479,19 +712,9 @@ gdouble gfs_vof_interpolate (FttCell * cell,
   if (l == level || GFS_IS_FULL (f))
     return f;
   else {
-    gdouble alpha = GFS_VARIABLE (cell, t->alpha->i);
-    gdouble h = ftt_level_size (level);
-    gdouble H = ftt_cell_size (cell);
-    FttComponent c;
-    FttVector m, q;
-    
-    for (c = 0; c < FTT_DIMENSION; c++)
-      (&m.x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
-    ftt_cell_pos (cell, &q);
-    alpha *= H;
-    for (c = 0; c < FTT_DIMENSION; c++)
-      alpha -= (&m.x)[c]*((&p->x)[c] - h/2. - (&q.x)[c] + H/2);
-    return gfs_plane_volume (&m, alpha/h);
+    FttVector m;
+    gdouble alpha = gfs_vof_plane_interpolate (cell, p, level, t, &m);
+    return gfs_plane_volume (&m, alpha);
   }
 }
 
@@ -1112,30 +1335,40 @@ gdouble gfs_vof_face_value (const FttCellFace * face, GfsVariableTracerVOF * t)
 }
 
 /**
- * gfs_vof_plane_facet:
+ * gfs_vof_facet:
  * @cell: a #FttCell.
- * @m: the plane normal.
- * @alpha: the plane position.
+ * @t: a #GfsVariableTracerVOF.
  * @p: a #FttVector array (of size 2 in 2D and 6 in 3D)
+ * @m: a #FttVector.
  *
  * Fills @p with the coordinates of points defining the
- * VOF-reconstructed interface facet defined by @m and @alpha.
+ * VOF-reconstructed interface facet defined by @t.
+ *
+ * Fills @m with the normal to the interface.
  *
  * Returns: the number of points defining the facet.
  */
-guint gfs_vof_plane_facet (FttCell * cell,
-			   FttVector * m,
-			   gdouble alpha,
-			   FttVector * p)
+guint gfs_vof_facet (FttCell * cell,
+		     GfsVariableTracerVOF * t,
+		     FttVector * p,
+		     FttVector * m)
 {
   g_return_val_if_fail (cell != NULL, 0);
-  g_return_val_if_fail (m != NULL, 0);
+  g_return_val_if_fail (t != NULL, 0);
   g_return_val_if_fail (p != NULL, 0);
+  g_return_val_if_fail (m != NULL, 0);
+
+  if (GFS_IS_FULL (GFS_VARIABLE (cell, GFS_VARIABLE1 (t)->i)))
+    return 0;
 
   guint n = 0;
   FttVector q;
   ftt_cell_pos (cell, &q);
   gdouble h = ftt_cell_size (cell);
+  FttComponent c;
+  for (c = 0; c < FTT_DIMENSION; c++)
+    (&m->x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
+  gdouble alpha = GFS_VARIABLE (cell, t->alpha->i);
 
 #if FTT_2D
   gdouble x, y;
@@ -1167,7 +1400,7 @@ guint gfs_vof_plane_facet (FttCell * cell,
   g_assert (n <= 2);
 #else /* 3D */
   gdouble max = fabs (m->x);
-  FttComponent c = FTT_X;
+  c = FTT_X;
   if (fabs (m->y) > max) {
     max = fabs (m->y);
     c = FTT_Y;
@@ -1187,78 +1420,13 @@ guint gfs_vof_plane_facet (FttCell * cell,
 }
 
 /**
- * gfs_vof_facet:
- * @cell: a #FttCell.
- * @t: a #GfsVariableTracerVOF.
- * @p: a #FttVector array (of size 2 in 2D and 6 in 3D)
- * @m: a #FttVector.
- *
- * Fills @p with the coordinates of points defining the
- * VOF-reconstructed interface facet defined by @t.
- *
- * Fills @m with the normal to the interface.
- *
- * Returns: the number of points defining the facet.
- */
-guint gfs_vof_facet (FttCell * cell,
-		     GfsVariableTracerVOF * t,
-		     FttVector * p,
-		     FttVector * m)
-{
-  g_return_val_if_fail (cell != NULL, 0);
-  g_return_val_if_fail (t != NULL, 0);
-  g_return_val_if_fail (p != NULL, 0);
-  g_return_val_if_fail (m != NULL, 0);
-
-  if (GFS_IS_FULL (GFS_VARIABLE (cell, GFS_VARIABLE1 (t)->i)))
-    return 0;
-  else {
-    FttComponent c;
-    for (c = 0; c < FTT_DIMENSION; c++)
-      (&m->x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
-    return gfs_vof_plane_facet (cell, m, GFS_VARIABLE (cell, t->alpha->i), p);
-  }
-}
-
-/**
- * gfs_vof_plane_center:
- * @cell: a #FttCell.
- * @m: the plane normal.
- * @alpha: the plane position.
- * @p: a #FttVector.
- *
- * Fills @p with the (approximate) coordinates of the center
- * of mass of the VOF-reconstructed interface facet defined by @m and @alpha.
- *
- * Returns: %TRUE if the cell contains the interface, %FALSE otherwise.
- */
-gboolean gfs_vof_plane_center (FttCell * cell, FttVector * m, gdouble alpha, FttVector * p)
-{
-  g_return_val_if_fail (cell != NULL, FALSE);
-  g_return_val_if_fail (m != NULL, FALSE);
-  g_return_val_if_fail (p != NULL, FALSE);
-
-  FttVector q[6];
-  guint i, nv = gfs_vof_plane_facet (cell, m, alpha, q);
-  if (nv > 0) {
-    p->x = p->y = p->z = 0.;
-    for (i = 0; i < nv; i++) {
-      p->x += q[i].x; p->y += q[i].y; p->z += q[i].z;
-    }
-    p->x /= nv; p->y /= nv; p->z /= nv;
-    return TRUE;
-  }
-  return FALSE;
-}
-
-/**
  * gfs_vof_center:
  * @cell: a #FttCell.
  * @t: a #GfsVariableTracerVOF.
  * @p: a #FttVector.
  *
- * Fills @p with the (approximate) coordinates of the center
- * of mass of the VOF-reconstructed interface facet defined by @t.
+ * Fills @p with the coordinates of the center of mass of the
+ * VOF-reconstructed interface facet defined by @t.
  *
  * Returns: %TRUE if the cell contains the interface, %FALSE otherwise.
  */
@@ -1268,17 +1436,19 @@ gboolean gfs_vof_center (FttCell * cell, GfsVariableTracerVOF * t, FttVector * p
   g_return_val_if_fail (t != NULL, FALSE);
   g_return_val_if_fail (p != NULL, 0);
 
-  FttVector m, q[6];
-  guint i, nv = gfs_vof_facet (cell, t, q, &m);
-  if (nv > 0) {
-    p->x = p->y = p->z = 0.;
-    for (i = 0; i < nv; i++) {
-      p->x += q[i].x; p->y += q[i].y; p->z += q[i].z;
-    }
-    p->x /= nv; p->y /= nv; p->z /= nv;
-    return TRUE;
-  }
-  return FALSE;
+  if (GFS_IS_FULL (GFS_VARIABLE (cell, GFS_VARIABLE1 (t)->i)))
+    return FALSE;
+
+  FttVector m, o;
+  FttComponent c;
+  for (c = 0; c < FTT_DIMENSION; c++)
+    (&m.x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
+  gfs_plane_area_center (&m, GFS_VARIABLE (cell, t->alpha->i), p);
+  ftt_cell_pos (cell, &o);
+  gdouble h = ftt_cell_size (cell);
+  for (c = 0; c < FTT_DIMENSION; c++)
+    (&p->x)[c] = (&o.x)[c] + h*((&p->x)[c] - 0.5);
+  return TRUE;
 }
 
 static gdouble fraction (FttVector * p,
@@ -1374,12 +1544,6 @@ static gboolean curvature_along_direction (FttCell * cell,
     found_all_heights = FALSE;
   else
     (*n)++;
-#if 0
-  if (H < -0.5 || H > 0.5) {
-    *kappa = G_MAXDOUBLE;
-    return TRUE;
-  }
-#endif
       
 #ifdef FTT_2D
   FttComponent cp = FTT_ORTHOGONAL_COMPONENT (c);
@@ -1446,38 +1610,98 @@ static gboolean curvature_along_direction (FttCell * cell,
 
 typedef struct {
   gdouble x2y, x4;
-  GtsVector o, m;
+  GtsVector o;
+#if FTT_2D
+  GtsVector m;
   gdouble a, b;
+#else /* 3D */
+  GtsVector4 M[3];
+  GtsVector t[3], rhs, a;
+#endif /* 3D */
 } ParabolaFit;
 
 static void parabola_fit_init (ParabolaFit * p, FttVector * o, FttVector * m)
 {
   p->o[0] = o->x; p->o[1] = o->y; p->o[2] = o->z;
+#if FTT_2D
   p->m[0] = m->x; p->m[1] = m->y; p->m[2] = m->z;
   gts_vector_normalize (p->m);
   p->x2y = p->x4 = p->b = 0.;
-}
-
-static void parabola_fit_destroy (ParabolaFit * p)
-{
+#else /* 3D */
+  gdouble max;
+  GtsVector nx = {0., 0., 0.}, ny, nz;
+  guint d = 0;
+
+  nz[0] = m->x; nz[1] = m->y; nz[2] = m->z;
+  gts_vector_normalize (nz);
+  max = nz[0]*nz[0];
+  /* build a vector orthogonal to nz */
+  if (nz[1]*nz[1] > max) { max = nz[1]*nz[1]; d = 1; }
+  if (nz[2]*nz[2] > max) d = 2;
+  switch (d) {
+  case 0: nx[0] = - nz[2]/nz[0]; nx[2] = 1.0; break;
+  case 1: nx[1] = - nz[2]/nz[1]; nx[2] = 1.0; break;
+  case 2: nx[2] = - nz[0]/nz[2]; nx[0] = 1.0; break;
+  }
+  gts_vector_normalize (nx);
+
+  /* build a second vector orthogonal to nx and nz */
+  gts_vector_cross (ny, nz, nx);
+
+  /* transformation matrix from (i,j,k) to (nx, ny, nz) */
+  p->t[0][0] = nx[0]; p->t[0][1] = nx[1]; p->t[0][2] = nx[2];
+  p->t[1][0] = ny[0]; p->t[1][1] = ny[1]; p->t[1][2] = ny[2];
+  p->t[2][0] = nz[0]; p->t[2][1] = nz[1]; p->t[2][2] = nz[2];
+
+  p->M[0][0] = 0.;
+  p->M[1][0] = p->M[1][1] = 0.;
+  p->M[2][0] = p->M[2][1] = p->M[2][2] = 0.;
+  p->rhs[0] = p->rhs[1] = p->rhs[2] = 0.;
+#endif /* 3D */
 }
 
 static void parabola_fit_add (ParabolaFit * p, GtsVector m, gdouble w)
 {
+#if FTT_2D
   gdouble x = m[0] - p->o[0];
   gdouble y = m[1] - p->o[1];
   gdouble x1 = p->m[1]*x - p->m[0]*y;
   gdouble y1 = p->m[0]*x + p->m[1]*y;
   p->x2y += w*x1*x1*y1;
   p->x4 += w*x1*x1*x1*x1;
+#else /* 3D */
+  gdouble x1 = m[0] - p->o[0];
+  gdouble y1 = m[1] - p->o[1];
+  gdouble z1 = m[2] - p->o[2];
+  gdouble x = p->t[0][0]*x1 + p->t[0][1]*y1 + p->t[0][2]*z1;
+  gdouble y = p->t[1][0]*x1 + p->t[1][1]*y1 + p->t[1][2]*z1;
+  gdouble z = p->t[2][0]*x1 + p->t[2][1]*y1 + p->t[2][2]*z1;
+  p->M[0][0] += w*x*x*x*x;
+  p->M[1][0] += w*x*x*y*y; p->M[1][1] += w*y*y*y*y;
+  p->M[2][0] += w*x*x*x*y; p->M[2][1] += w*x*y*y*y; p->M[2][2] += w*x*x*y*y;
+  p->rhs[0] += w*z*x*x; p->rhs[1] += w*z*y*y; p->rhs[2] += w*z*x*y;
+#endif /* 3D */
 }
 
 static void parabola_fit_solve (ParabolaFit * p)
 {
+#if FTT_2D
   if (p->x4 > 0.)
     p->a = p->x2y/p->x4;
   else /* this is a degenerate/isolated interface fragment */
     p->a = 0.;
+#else /* 3D */
+  p->M[0][1] = p->M[1][0]; p->M[0][2] = p->M[2][0];
+  p->M[1][2] = p->M[2][1];  
+  GtsMatrix * M = gts_matrix3_inverse ((GtsMatrix *) p->M);
+  if (M) {
+    p->a[0] = M[0][0]*p->rhs[0] + M[0][1]*p->rhs[1] + M[0][2]*p->rhs[2];
+    p->a[1] = M[1][0]*p->rhs[0] + M[1][1]*p->rhs[1] + M[1][2]*p->rhs[2];
+    gts_matrix_destroy (M);
+  }
+  else /* this is a degenerate/isolated interface fragment */
+    p->a[0] = p->a[1] = 0.;
+#endif /* 3D */
 }
 
 static gdouble parabola_fit_curvature (ParabolaFit * p, gdouble kappamax)
@@ -1486,18 +1710,28 @@ static gdouble parabola_fit_curvature (ParabolaFit * p, gdouble kappamax)
 #if FTT_2D
   kappa = 2.*p->a;
 #else
-  g_assert_not_implemented ();
+  kappa = 2.*(p->a[0] + p->a[1]);
 #endif
   if (fabs (kappa) > kappamax)
     return kappa > 0. ? kappamax : - kappamax;
   return kappa;
 }
 
-static void add_vof_center (ParabolaFit * fit, FttCell * cell, GfsVariable * v, gdouble w)
+static void add_vof_center (FttCell * cell, FttVector * p, guint level,
+			    GfsVariableTracerVOF * t,
+			    ParabolaFit * fit, gdouble w)
 {
-  FttVector c;
-  if (gfs_vof_center (cell, GFS_VARIABLE_TRACER_VOF (v), &c))
+  gdouble f = GFS_VARIABLE (cell, GFS_VARIABLE1 (t)->i);
+  if (!GFS_IS_FULL (f)) {
+    FttVector m, c;
+    gdouble alpha = gfs_vof_plane_interpolate (cell, p, level, t, &m);
+    gfs_plane_area_center (&m, alpha, &c);
+    gdouble h = ftt_level_size (level);
+    FttComponent i;
+    for (i = 0; i < FTT_DIMENSION; i++)
+      (&c.x)[i] = (&p->x)[i] + h*((&c.x)[i] - 0.5);
     parabola_fit_add (fit, &c.x, w);
+  }
 }
 
 static void fit_from_fractions (FttCell * cell, GfsVariable * v, ParabolaFit * fit)
@@ -1506,23 +1740,22 @@ static void fit_from_fractions (FttCell * cell, GfsVariable * v, ParabolaFit * f
   guint level = ftt_cell_level (cell);
   gint x, y, z = 0;
   FttVector p;
-  #define SMALL 0.1
-  static gdouble w[3][3] = {
-    { SMALL, 1., SMALL },
-    {  1.,   0.,  1.   },
-    { SMALL, 1., SMALL }
-  };
+  static gdouble w[4] = {1., 1., 0.1, 0.01 };
   
   ftt_cell_pos (cell, &p);
-  for (x = -1; x <= 1; x++)
-    for (y = -1; y <= 1; y++)
-      if (x != 0 || y != 0) {
-	FttVector o;
-	o.x = p.x + h*x; o.y = p.y + h*y; o.z = p.z + h*z;
-	FttCell * neighbor = domain_and_boundary_locate (v->domain, o, level);
-	if (neighbor)
-	  add_vof_center (fit, neighbor, v, w[x + 1][y + 1]);
-      }
+#if !FTT_2D
+  for (z = -1; z <= 1; z++)
+#endif
+    for (x = -1; x <= 1; x++)
+      for (y = -1; y <= 1; y++)
+	if (x != 0 || y != 0 || z != 0) {
+	  FttVector o;
+	  o.x = p.x + h*x; o.y = p.y + h*y; o.z = p.z + h*z;
+	  FttCell * neighbor = domain_and_boundary_locate (v->domain, o, level);
+	  if (neighbor)
+	    add_vof_center (neighbor, &o, level, GFS_VARIABLE_TRACER_VOF (v),
+			    fit, w[abs (x) + abs (y) + abs (z)]);
+	}
 }
 
 static gdouble gfs_curvature_fit (FttCell * cell, GfsVariableTracerVOF * t)
@@ -1545,7 +1778,6 @@ static gdouble gfs_curvature_fit (FttCell * cell, GfsVariableTracerVOF * t)
   fit_from_fractions (cell, GFS_VARIABLE1 (t), &fit);
   parabola_fit_solve (&fit);
   gdouble kappa = parabola_fit_curvature (&fit, 2./ftt_cell_size (cell));
-  parabola_fit_destroy (&fit);
   return kappa;
 }
 
@@ -1577,9 +1809,7 @@ static void orientation (FttVector * m, FttComponent * c)
  * An implementation of the Height-Function (HF) method generalised to
  * adaptive meshes.
  *
- * Returns: the curvature of the interface contained in @cell or
- * G_MAXDOUBLE if the curvature cannot be computed using the HF
- * method.
+ * Returns: the curvature of the interface contained in @cell.
  */
 gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
 {
@@ -1597,7 +1827,7 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
   FttComponent try[FTT_DIMENSION];
   orientation (&m, try); /* sort directions according to normal */
 
-  gdouble kappa = G_MAXDOUBLE;
+  gdouble kappa = 0.;
   GtsVector interface[FTT_DIMENSION*NI];
   guint n = 0;
   for (c = 0; c < FTT_DIMENSION; c++) /* try each direction */
@@ -1609,23 +1839,20 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
 
   if (n < NI)
     return gfs_curvature_fit (cell, t);
-  else {
-    gdouble h = ftt_cell_size (cell);
-    ParabolaFit fit;
-    guint j;
-    
-    FttVector p, fc;
-    ftt_cell_pos (cell, &p);
-    gfs_vof_center (cell, t, &fc);
-    fc.x = (fc.x - p.x)/h;
-    fc.y = (fc.y - p.y)/h;
-    fc.z = (fc.z - p.z)/h;
-    parabola_fit_init (&fit, &fc, &m);
-    for (j = 0; j < n; j++)
-      parabola_fit_add (&fit, interface[j], 1.);
-    parabola_fit_solve (&fit);
-    kappa = parabola_fit_curvature (&fit, 2.)/h;
-    parabola_fit_destroy (&fit);
-  }
-  return kappa;
+
+  gdouble h = ftt_cell_size (cell);
+  ParabolaFit fit;
+  guint j;
+  
+  FttVector p, fc;
+  ftt_cell_pos (cell, &p);
+  gfs_vof_center (cell, t, &fc);
+  fc.x = (fc.x - p.x)/h;
+  fc.y = (fc.y - p.y)/h;
+  fc.z = (fc.z - p.z)/h;
+  parabola_fit_init (&fit, &fc, &m);
+  for (j = 0; j < n; j++)
+    parabola_fit_add (&fit, interface[j], 1.);
+  parabola_fit_solve (&fit);
+  return parabola_fit_curvature (&fit, 2.)/h;
 }
diff --git a/src/vof.h b/src/vof.h
index 68ed68e..dbad090 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -35,12 +35,16 @@ void    gfs_line_center            (FttVector * m,
 				    gdouble alpha, 
 				    gdouble a, 
 				    FttVector * p);
+void    gfs_line_area_center       (FttVector * m, 
+				    gdouble alpha, 
+				    FttVector * p);
 gdouble gfs_line_alpha             (FttVector * m, 
 				    gdouble c);
 #if FTT_2D
 #  define gfs_plane_volume         gfs_line_area
 #  define gfs_plane_alpha          gfs_line_alpha
 #  define gfs_plane_center         gfs_line_center
+#  define gfs_plane_area_center     gfs_line_area_center
 #else /* 3D */
 gdouble gfs_plane_volume           (FttVector * m, 
 				    gdouble alpha);
@@ -50,6 +54,9 @@ void    gfs_plane_center           (FttVector * m,
 				    gdouble alpha, 
 				    gdouble a,
 				    FttVector * p);
+void    gfs_plane_area_center      (FttVector * m, 
+				    gdouble alpha, 
+				    FttVector * p);
 #endif /* 3D */
 void    gfs_youngs_gradient        (FttCell * cell, 
 				    GfsVariable * v,
@@ -79,21 +86,18 @@ void     gfs_tracer_vof_advection  (GfsDomain * domain,
 				    GfsAdvectionParams * par);
 gdouble  gfs_vof_face_value        (const FttCellFace * face, 
 				    GfsVariableTracerVOF * t);
-guint    gfs_vof_plane_facet       (FttCell * cell,
-				    FttVector * m,
-				    gdouble alpha,
-				    FttVector * p);
 guint    gfs_vof_facet             (FttCell * cell,
 				    GfsVariableTracerVOF * t,
 				    FttVector * p,
 				    FttVector * m);
-gboolean gfs_vof_plane_center      (FttCell * cell, 
-				    FttVector * m, 
-				    gdouble alpha, 
-				    FttVector * p);
 gboolean gfs_vof_center            (FttCell * cell,
 				    GfsVariableTracerVOF * t,
 				    FttVector * p);
+gdouble  gfs_vof_plane_interpolate (FttCell * cell,
+				    FttVector * p,
+				    guint level,
+				    GfsVariableTracerVOF * t,
+				    FttVector * m);
 gdouble  gfs_vof_interpolate       (FttCell * cell,
 				    FttVector * p,
 				    guint level,

-- 
Gerris Flow Solver



More information about the debian-science-commits mailing list