[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 b627903e5306998be3ca6266213e69be5eb0af67
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed May 23 10:41:23 2007 +1000
VariableCurvature optionally computes the maximum curvature
darcs-hash:20070523004123-d4795-3d7ec7c7a342dfd993b66a79543873efc38f65e5.gz
diff --git a/src/tension.c b/src/tension.c
index 4d8fa03..1aaded8 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -327,6 +327,43 @@ GfsSourceGenericClass * gfs_source_tension_class (void)
/* GfsVariableCurvature: object */
+static void variable_curvature_destroy (GtsObject * o)
+{
+ if (GFS_VARIABLE_CURVATURE (o)->kmax)
+ gts_object_destroy (GTS_OBJECT (GFS_VARIABLE_CURVATURE (o)->kmax));
+
+ (* GTS_OBJECT_CLASS (gfs_variable_curvature_class ())->parent_class->destroy) (o);
+}
+
+static void curvature_coarse_fine (FttCell * parent, GfsVariable * v)
+{
+ FttCellChildren child;
+ guint n;
+
+ ftt_cell_children (parent, &child);
+ for (n = 0; n < FTT_CELLS; n++)
+ if (child.c[n])
+ GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
+}
+
+static void curvature_fine_coarse (FttCell * parent, GfsVariable * v)
+{
+ FttCellChildren child;
+ gdouble val = 0., sa = 0.;
+ guint i;
+
+ ftt_cell_children (parent, &child);
+ for (i = 0; i < FTT_CELLS; i++)
+ if (child.c[i] && GFS_VARIABLE (child.c[i], v->i) < G_MAXDOUBLE) {
+ val += GFS_VARIABLE (child.c[i], v->i);
+ sa += 1.;
+ }
+ if (sa > 0.)
+ GFS_VARIABLE (parent, v->i) = val/sa;
+ else
+ GFS_VARIABLE (parent, v->i) = G_MAXDOUBLE;
+}
+
static void variable_curvature_read (GtsObject ** o, GtsFile * fp)
{
GfsVariableCurvature * v = GFS_VARIABLE_CURVATURE (*o);
@@ -351,21 +388,30 @@ static void variable_curvature_read (GtsObject ** o, GtsFile * fp)
if (!GFS_IS_VARIABLE_TRACER_VOF (v->f)) {
gts_file_error (fp, "variable `%s' is not a VOF tracer", fp->token->str);
return;
- }
+ }
GFS_VARIABLE1 (v)->description = g_strjoin (" ",
"Curvature of the interface defined by tracer",
v->f->name, NULL);
+ gts_file_next_token (fp);
+ if (fp->type == GTS_STRING) {
+ if ((v->kmax = gfs_variable_from_name (domain->variables, fp->token->str)) ||
+ (v->kmax = gfs_domain_add_variable (domain, fp->token->str, "Maximum curvature"))) {
+ v->kmax->coarse_fine = curvature_coarse_fine;
+ v->kmax->fine_coarse = curvature_fine_coarse;
+ gts_file_next_token (fp);
+ }
+ }
}
- else if (GFS_IS_VARIABLE_DISTANCE (v->f))
+ else if (GFS_IS_VARIABLE_DISTANCE (v->f)) {
GFS_VARIABLE1 (v)->description = g_strjoin (" ",
"Curvature of the interface defined by distance",
- v->f->name, NULL);
+ v->f->name, NULL);
+ gts_file_next_token (fp);
+ }
else {
- GFS_VARIABLE1 (v)->description = NULL;
gts_file_error (fp, "variable `%s' is neither a tracer nor a distance", fp->token->str);
return;
}
- gts_file_next_token (fp);
}
static void variable_curvature_write (GtsObject * o, FILE * fp)
@@ -375,17 +421,30 @@ static void variable_curvature_write (GtsObject * o, FILE * fp)
(* GTS_OBJECT_CLASS (gfs_variable_curvature_class ())->parent_class->write) (o, fp);
fprintf (fp, " %s", v->f->name);
+ if (v->kmax)
+ fprintf (fp, " %s", v->kmax->name);
}
static void height_curvature (FttCell * cell, GfsVariable * v)
{
GfsVariable * t = GFS_VARIABLE_CURVATURE (v)->f;
+ GfsVariable * kmax = GFS_VARIABLE_CURVATURE (v)->kmax;
gdouble f = GFS_VARIABLE (cell, t->i);
- if (GFS_IS_FULL (f))
+ if (GFS_IS_FULL (f)) {
GFS_VARIABLE (cell, v->i) = G_MAXDOUBLE;
- else
- GFS_VARIABLE (cell, v->i) = gfs_height_curvature (cell, GFS_VARIABLE_TRACER_VOF (t));
+ if (kmax)
+ GFS_VARIABLE (cell, kmax->i) = G_MAXDOUBLE;
+ }
+ else {
+ if (kmax) {
+ gdouble k;
+ GFS_VARIABLE (cell, v->i) = gfs_height_curvature (cell, GFS_VARIABLE_TRACER_VOF (t), &k);
+ GFS_VARIABLE (cell, kmax->i) = k;
+ }
+ else
+ GFS_VARIABLE (cell, v->i) = gfs_height_curvature (cell, GFS_VARIABLE_TRACER_VOF (t), NULL);
+ }
}
static void fit_curvature (FttCell * cell, GfsVariable * v)
@@ -393,8 +452,16 @@ static void fit_curvature (FttCell * cell, GfsVariable * v)
GfsVariable * t = GFS_VARIABLE_CURVATURE (v)->f;
gdouble f = GFS_VARIABLE (cell, t->i);
- if (!GFS_IS_FULL (f) && GFS_VARIABLE (cell, v->i) == G_MAXDOUBLE)
- GFS_VARIABLE (cell, v->i) = gfs_fit_curvature (cell, GFS_VARIABLE_TRACER_VOF (t));
+ if (!GFS_IS_FULL (f) && GFS_VARIABLE (cell, v->i) == G_MAXDOUBLE) {
+ GfsVariable * kmax = GFS_VARIABLE_CURVATURE (v)->kmax;
+ if (kmax) {
+ gdouble k;
+ GFS_VARIABLE (cell, v->i) = gfs_fit_curvature (cell, GFS_VARIABLE_TRACER_VOF (t), &k);
+ GFS_VARIABLE (cell, kmax->i) = k;
+ }
+ else
+ GFS_VARIABLE (cell, v->i) = gfs_fit_curvature (cell, GFS_VARIABLE_TRACER_VOF (t), NULL);
+ }
}
typedef struct {
@@ -423,11 +490,11 @@ static void diffuse (FttCell * cell, DiffuseParms * p)
}
}
-static void variable_curvature_diffuse (GfsEvent * event, GfsSimulation * sim, guint n)
+static void variable_curvature_diffuse (GfsVariable * v, GfsSimulation * sim, guint n)
{
GfsDomain * domain = GFS_DOMAIN (sim);
DiffuseParms p;
- p.v = GFS_VARIABLE1 (event);
+ p.v = v;
p.tmp = gfs_temporary_variable (domain);
while (n--) {
@@ -445,6 +512,7 @@ static void variable_curvature_diffuse (GfsEvent * event, GfsSimulation * sim, g
static void variable_curvature_from_fraction (GfsEvent * event, GfsSimulation * sim)
{
GfsDomain * domain = GFS_DOMAIN (sim);
+ GfsVariable * kmax = GFS_VARIABLE_CURVATURE (event)->kmax;
gfs_domain_timer_start (domain, "variable_curvature");
@@ -453,13 +521,25 @@ 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, 1);
+ if (kmax) {
+ gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) kmax->fine_coarse, kmax);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, kmax);
+ variable_curvature_diffuse (kmax, sim, 1);
+ }
+ variable_curvature_diffuse (GFS_VARIABLE1 (event), sim, 1);
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
(FttCellTraverseFunc) fit_curvature, event);
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, 1);
+ if (kmax) {
+ gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) kmax->fine_coarse, kmax);
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, kmax);
+ variable_curvature_diffuse (kmax, sim, 1);
+ }
+ variable_curvature_diffuse (GFS_VARIABLE1 (event), sim, 1);
gfs_domain_timer_stop (domain, "variable_curvature");
}
@@ -549,7 +629,7 @@ static void variable_curvature_from_distance (GfsEvent * event, GfsSimulation *
for (c = 0; c < FTT_DIMENSION + 1; c++)
gts_object_destroy (GTS_OBJECT (n[c]));
- variable_curvature_diffuse (event, sim, 2);
+ variable_curvature_diffuse (GFS_VARIABLE1 (event), sim, 2);
gfs_domain_timer_stop (domain, "variable_curvature");
}
@@ -569,40 +649,12 @@ static gboolean variable_curvature_event (GfsEvent * event, GfsSimulation * sim)
static void variable_curvature_class_init (GtsObjectClass * klass)
{
+ klass->destroy = variable_curvature_destroy;
klass->read = variable_curvature_read;
klass->write = variable_curvature_write;
GFS_EVENT_CLASS (klass)->event = variable_curvature_event;
}
-static void curvature_coarse_fine (FttCell * parent, GfsVariable * v)
-{
- FttCellChildren child;
- guint n;
-
- ftt_cell_children (parent, &child);
- for (n = 0; n < FTT_CELLS; n++)
- if (child.c[n])
- GFS_VARIABLE (child.c[n], v->i) = GFS_VARIABLE (parent, v->i);
-}
-
-static void curvature_fine_coarse (FttCell * parent, GfsVariable * v)
-{
- FttCellChildren child;
- gdouble val = 0., sa = 0.;
- guint i;
-
- ftt_cell_children (parent, &child);
- for (i = 0; i < FTT_CELLS; i++)
- if (child.c[i] && GFS_VARIABLE (child.c[i], v->i) < G_MAXDOUBLE) {
- val += GFS_VARIABLE (child.c[i], v->i);
- sa += 1.;
- }
- if (sa > 0.)
- GFS_VARIABLE (parent, v->i) = val/sa;
- else
- GFS_VARIABLE (parent, v->i) = G_MAXDOUBLE;
-}
-
static void variable_curvature_init (GfsVariable * v)
{
v->coarse_fine = curvature_coarse_fine;
@@ -704,7 +756,7 @@ static gboolean variable_position_event (GfsEvent * event, GfsSimulation * sim)
(FttCellTraverseFunc) GFS_VARIABLE1 (event)->fine_coarse, event);
gfs_domain_bc (domain, FTT_TRAVERSE_ALL, -1, GFS_VARIABLE1 (event));
- variable_curvature_diffuse (event, sim, 2);
+ variable_curvature_diffuse (GFS_VARIABLE1 (event), sim, 2);
gfs_domain_timer_stop (domain, "variable_position");
return TRUE;
diff --git a/src/tension.h b/src/tension.h
index 6fcbcdd..c11e47c 100644
--- a/src/tension.h
+++ b/src/tension.h
@@ -100,7 +100,7 @@ struct _GfsVariableCurvature {
GfsVariable parent;
/*< public >*/
- GfsVariable * f;
+ GfsVariable * f, * kmax;
};
#define GFS_VARIABLE_CURVATURE(obj) GTS_OBJECT_CAST (obj,\
diff --git a/src/vof.c b/src/vof.c
index d169f41..1370b60 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -1535,6 +1535,7 @@ static gboolean curvature_along_direction (FttCell * cell,
GfsVariableTracerVOF * t,
FttComponent c,
gdouble * kappa,
+ gdouble * kmax,
GtsVector * interface,
guint * n)
{
@@ -1586,6 +1587,8 @@ static gboolean curvature_along_direction (FttCell * cell,
gdouble hx = (h[0] - h[1])/2.;
gdouble dnm = 1. + hx*hx;
*kappa = hxx/(size*sqrt (dnm*dnm*dnm));
+ if (kmax)
+ *kmax = fabs (*kappa);
}
#else /* 3D */
static FttComponent or[3][2] = { { FTT_Y, FTT_Z }, { FTT_X, FTT_Z }, { FTT_X, FTT_Y } };
@@ -1610,9 +1613,17 @@ static gboolean curvature_along_direction (FttCell * cell,
gdouble hyy = h[1][2] - 2.*H + h[1][0];
gdouble hx = (h[2][1] - h[0][1])/2.;
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 hxy = (h[2][2] + h[0][0] - h[2][0] - h[0][2])/4.;
gdouble dnm = 1. + hx*hx + hy*hy;
- *kappa = (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 - 2.*hxy*hx*hy)/(size*sqrt (dnm*dnm*dnm));
+ if (kmax) {
+ gdouble km = *kappa/2.;
+ /* Gaussian curvature */
+ gdouble kg = (hxx*hyy - hxy*hxy)/(size*size*dnm*dnm);
+ gdouble a = km*km - kg;
+ g_assert (a >= 0.);
+ *kmax = fabs (km) + sqrt (a);
+ }
}
#endif /* 3D */
return found_all_heights;
@@ -1746,10 +1757,11 @@ static void parabola_fit_solve (ParabolaFit * p)
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];
+ p->a[2] = M[2][0]*p->rhs[0] + M[2][1]*p->rhs[1] + M[2][2]*p->rhs[2];
gts_matrix_destroy (M);
}
else /* this may be a degenerate/isolated interface fragment */
- p->a[0] = p->a[1] = 0.;
+ p->a[0] = p->a[1] = p->a[2] = 0.;
# else
p->M[0][1] = p->M[2][2]; p->M[0][5] = p->M[3][3];
p->M[1][5] = p->M[4][4];
@@ -1768,33 +1780,46 @@ static void parabola_fit_solve (ParabolaFit * p)
}
else { /* this may be a degenerate/isolated interface fragment */
g_warning ("singular matrix");
- p->a[0] = p->a[1] = 0.;
+ p->a[0] = p->a[1] = p->a[2] = 0.;
}
# endif
#endif /* 3D */
}
-static gdouble parabola_fit_curvature (ParabolaFit * p, gdouble kappamax)
+static gdouble parabola_fit_curvature (ParabolaFit * p, gdouble kappamax,
+ gdouble * kmax)
{
gdouble kappa;
#if FTT_2D
gdouble dnm = 1. + p->a[1]*p->a[1];
kappa = 2.*p->a[0]/sqrt (dnm*dnm*dnm);
-#else
-# if PARABOLA_SIMPLER
- kappa = 2.*(p->a[0] + p->a[1]);
-# else
+ if (kmax)
+ *kmax = fabs (kappa);
+#else /* 3D */
gdouble hxx = 2.*p->a[0];
gdouble hyy = 2.*p->a[1];
- gdouble hx = p->a[3];
- gdouble hy = p->a[4];
gdouble hxy = p->a[2];
- gdouble dnm = 1. + hx*hx + hy*hy;
- kappa = (hxx + hyy + hxx*hy*hy + hyy*hx*hx - hxy*hx*hy)/sqrt (dnm*dnm*dnm);
+ gdouble hx, hy;
+# if PARABOLA_SIMPLER
+ hx = hy = 0.;
+# else
+ hx = p->a[3];
+ hy = p->a[4];
# endif
-#endif
- if (fabs (kappa) > kappamax)
+ gdouble dnm = 1. + hx*hx + hy*hy;
+ kappa = (hxx + hyy + hxx*hy*hy + hyy*hx*hx - 2.*hxy*hx*hy)/sqrt (dnm*dnm*dnm);
+ if (kmax) {
+ gdouble kg = (hxx*hyy - hxy*hxy)/(dnm*dnm);
+ gdouble a = kappa*kappa/4. - kg;
+ g_assert (a >= 0.);
+ *kmax = fabs (kappa/2.) + sqrt (a);
+ }
+#endif /* 3D */
+ if (fabs (kappa) > kappamax) {
+ if (kmax)
+ *kmax = kappamax;
return kappa > 0. ? kappamax : - kappamax;
+ }
return kappa;
}
@@ -1852,14 +1877,19 @@ static void fit_from_fractions (FttCell * cell, GfsVariable * v, ParabolaFit * f
* gfs_fit_curvature:
* @cell: a #FttCell containing an interface.
* @v: a #GfsVariableTracerVOF.
+ * @kmax: a pointer or %NULL.
*
* Computes an approximation of the curvature of the interface
- * contained in @cell using ellipsoid fitting of the centroids of the
+ * contained in @cell using paraboloid fitting of the centroids of the
* reconstructed interface segments.
*
- * Returns: the curvature of the interface contained in @cell.
+ * If @kmax is not %NULL, it is filled with the absolute value of the
+ * maximum surface curvature (note that in 2D this is just the absolute value of
+ * the mean curvature).
+ *
+ * Returns: (double in 3D) the mean curvature of the interface contained in @cell.
*/
-gdouble gfs_fit_curvature (FttCell * cell, GfsVariableTracerVOF * t)
+gdouble gfs_fit_curvature (FttCell * cell, GfsVariableTracerVOF * t, gdouble * kmax)
{
g_return_val_if_fail (cell != NULL, 0.);
g_return_val_if_fail (t != NULL, 0.);
@@ -1884,7 +1914,9 @@ gdouble gfs_fit_curvature (FttCell * cell, GfsVariableTracerVOF * t)
parabola_fit_add (&fit, &fc.x, area);
fit_from_fractions (cell, GFS_VARIABLE1 (t), &fit);
parabola_fit_solve (&fit);
- gdouble kappa = parabola_fit_curvature (&fit, 2.)/h;
+ gdouble kappa = parabola_fit_curvature (&fit, 2., kmax)/h;
+ if (kmax)
+ *kmax /= h;
parabola_fit_destroy (&fit);
return kappa;
}
@@ -1934,14 +1966,18 @@ static guint independent_positions (GtsVector * interface, guint n)
* gfs_height_curvature:
* @cell: a #FttCell containing an interface.
* @v: a #GfsVariableTracerVOF.
- * @k: a curvature vector.
+ * @kmax: a pointer or %NULL.
*
* An implementation of the Height-Function (HF) method generalised to
* adaptive meshes.
*
- * Returns: the curvature of the interface contained in @cell.
+ * If @kmax is not %NULL, it is filled with the absolute value of the
+ * maximum surface curvature (note that in 2D this is just the absolute value of
+ * the mean curvature).
+ *
+ * Returns: (double in 3D) the mean curvature of the interface contained in @cell.
*/
-gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
+gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t, gdouble * kmax)
{
g_return_val_if_fail (cell != NULL, 0.);
g_return_val_if_fail (t != NULL, 0.);
@@ -1962,7 +1998,7 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
GtsVector interface[FTT_DIMENSION*NI];
guint n = 0;
for (c = 0; c < FTT_DIMENSION; c++) /* try each direction */
- if (curvature_along_direction (cell, t, try[c], &kappa, interface, &n))
+ if (curvature_along_direction (cell, t, try[c], &kappa, kmax, interface, &n))
return kappa;
/* Could not compute curvature from the simple algorithm along any direction:
@@ -1990,7 +2026,9 @@ gdouble gfs_height_curvature (FttCell * cell, GfsVariableTracerVOF * t)
for (j = 0; j < n; j++)
parabola_fit_add (&fit, interface[j], 1.);
parabola_fit_solve (&fit);
- kappa = parabola_fit_curvature (&fit, 2.)/h;
+ kappa = parabola_fit_curvature (&fit, 2., kmax)/h;
+ if (kmax)
+ *kmax /= h;
parabola_fit_destroy (&fit);
return kappa;
}
diff --git a/src/vof.h b/src/vof.h
index 762bf0c..194f396 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -103,9 +103,11 @@ gdouble gfs_vof_interpolate (FttCell * cell,
guint level,
GfsVariableTracerVOF * t);
gdouble gfs_height_curvature (FttCell * cell,
- GfsVariableTracerVOF * t);
+ GfsVariableTracerVOF * t,
+ gdouble * kmax);
gdouble gfs_fit_curvature (FttCell * cell,
- GfsVariableTracerVOF * t);
+ GfsVariableTracerVOF * t,
+ gdouble * kmax);
#ifdef __cplusplus
}
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list