[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:54:26 UTC 2009
The following commit has been merged in the upstream branch:
commit f981db354868ee390dd35e0d93d02d0b7d05bee6
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Apr 30 20:15:19 2007 +1000
More robust implementation of gfs_height_curvature()
Tries all directions of integration (based on normal orientation) rather than
only the first guess.
darcs-hash:20070430101519-d4795-de9ff583878ff55c3524fe153459544ac6c8ff08.gz
diff --git a/src/vof.c b/src/vof.c
index 1bcae33..e6309dd 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -1337,54 +1337,45 @@ static guint local_height (FttVector * p,
return n;
}
-static FttComponent orientation (FttVector * m)
+static void orientation (FttVector * m, FttComponent * c)
{
- gdouble max = fabs (m->x);
- FttComponent c = FTT_X, i;
- for (i = 1; i < FTT_DIMENSION; i++)
- if (fabs ((&m->x)[i]) > max) {
- max = fabs ((&m->x)[i]);
- c = i;
- }
- return c;
+ FttComponent i, j;
+ for (i = 0; i < FTT_DIMENSION; i++)
+ c[i] = i;
+ for (i = 0; i < FTT_DIMENSION - 1; i++)
+ for (j = 0; j < FTT_DIMENSION - 1 - i; j++)
+ if (fabs ((&m->x)[c[j + 1]]) > fabs ((&m->x)[c[j]])) {
+ FttComponent tmp = c[j];
+ c[j] = c[j + 1];
+ c[j + 1] = tmp;
+ }
}
-/**
- * gfs_height_curvature:
- * @cell: a #FttCell containing an interface.
- * @v: a #GfsVariableTracerVOF.
- *
- * 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.
- */
-gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
+static gboolean curvature_along_direction (FttCell * cell,
+ GfsVariableTracerVOF * t,
+ FttComponent c,
+ gdouble * kappa)
{
- g_return_val_if_fail (cell != NULL, 0.);
- g_return_val_if_fail (t != NULL, 0.);
-
GfsVariable * v = GFS_VARIABLE1 (t);
- g_return_val_if_fail (!GFS_IS_FULL (GFS_VARIABLE (cell, v->i)), 0.);
FttVector m;
- FttComponent c;
- for (c = 0; c < FTT_DIMENSION; c++)
- (&m.x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
- c = orientation (&m);
+ FttComponent i;
+ for (i = 0; i < FTT_DIMENSION; i++)
+ (&m.x)[i] = GFS_VARIABLE (cell, t->m[i]->i);
FttVector p;
ftt_cell_pos (cell, &p);
guint level = ftt_cell_level (cell);
+ gdouble size = ftt_level_size (level);
+
gdouble H;
if (!local_height (&p, &p, level, v, c, &H))
- return G_MAXDOUBLE;
- if (H < -0.5 || H > 0.5)
- return G_MAXDOUBLE;
+ return FALSE;
+ if (H < -0.5 || H > 0.5) {
+ *kappa = G_MAXDOUBLE;
+ return TRUE;
+ }
- gdouble size = ftt_level_size (level);
#ifdef FTT_2D
FttComponent cp = FTT_ORTHOGONAL_COMPONENT (c);
gdouble h[2];
@@ -1393,26 +1384,19 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
(&q.x)[cp] += size;
(&q.x)[c] -= slope*size;
- // fprintf (stderr, "\n (%g %g) ", q.x, q.y);
- if (!local_height (&q, &p, level, v, c, &h[0])) {
- g_warning ("Failed to compute local height at (%g,%g)", q.x, q.y);
- return G_MAXDOUBLE;
- }
+ if (!local_height (&q, &p, level, v, c, &h[0]))
+ return FALSE;
q = p;
(&q.x)[cp] -= size;
(&q.x)[c] += slope*size;
- // fprintf (stderr, "\n (%g %g) ", q.x, q.y);
- if (!local_height (&q, &p, level, v, c, &h[1])) {
- g_warning ("Failed to compute local height at (%g,%g)", q.x, q.y);
- return G_MAXDOUBLE;
- }
+ if (!local_height (&q, &p, level, v, c, &h[1]))
+ return FALSE;
gdouble hxx = h[0] - 2*H + h[1];
gdouble hx = (h[0] - h[1])/2.;
gdouble dnm = 1. + hx*hx;
- // fprintf (stderr, " %g\n", hxx/(size*sqrt (dnm*dnm*dnm)));
- return hxx/(size*sqrt (dnm*dnm*dnm));
+ *kappa = hxx/(size*sqrt (dnm*dnm*dnm));
#else /* 3D */
static FttComponent or[3][2] = { { FTT_Y, FTT_Z }, { FTT_X, FTT_Z }, { FTT_X, FTT_Y } };
gdouble h[3][3];
@@ -1427,11 +1411,8 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
(&q.x)[or[c][0]] += size*x;
(&q.x)[or[c][1]] += size*y;
(&q.x)[c] -= slope*size;
- guint n = local_height (&q, &p, level, v, c, &h[x + 1][y + 1]);
- if (!n) {
- g_warning ("Failed to compute local height at (%g,%g,%g)", q.x, q.y, q.z);
- return G_MAXDOUBLE;
- }
+ if (!local_height (&q, &p, level, v, c, &h[x + 1][y + 1]))
+ return FALSE;
}
gdouble hxx = h[2][1] - 2.*H + h[0][1];
@@ -1440,6 +1421,45 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
gdouble hy = (h[1][2] - h[1][0])/2.;
gdouble hxy = (h[2][2] + h[0][0] - h[2][0] - h[0][2])/2.;
gdouble dnm = 1. + hx*hx + hy*hy;
- return (hxx + hyy + hxx*hy*hy + hyy*hx*hx - hxy*hx*hy)/(size*sqrt (dnm*dnm*dnm));
+ *kappa = (hxx + hyy + hxx*hy*hy + hyy*hx*hx - hxy*hx*hy)/(size*sqrt (dnm*dnm*dnm));
#endif /* 3D */
+ return TRUE;
+}
+
+/**
+ * gfs_height_curvature:
+ * @cell: a #FttCell containing an interface.
+ * @v: a #GfsVariableTracerVOF.
+ *
+ * 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.
+ */
+gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
+{
+ g_return_val_if_fail (cell != NULL, 0.);
+ g_return_val_if_fail (t != NULL, 0.);
+
+ GfsVariable * v = GFS_VARIABLE1 (t);
+ g_return_val_if_fail (!GFS_IS_FULL (GFS_VARIABLE (cell, v->i)), 0.);
+
+ FttVector m;
+ FttComponent c;
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&m.x)[c] = GFS_VARIABLE (cell, t->m[c]->i);
+
+ FttComponent try[FTT_DIMENSION];
+ orientation (&m, try); /* sort directions according to normal */
+ gdouble kappa = G_MAXDOUBLE;
+ for (c = 0; c < FTT_DIMENSION; c++) /* try each direction */
+ if (curvature_along_direction (cell, t, try[c], &kappa))
+ return kappa;
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ g_warning ("Failed to compute curvature at (%g,%g,%g) level: %d",
+ p.x, p.y, p.z, ftt_cell_level (cell));
+ return G_MAXDOUBLE;
}
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list