[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:40 UTC 2009
The following commit has been merged in the upstream branch:
commit 035c7a68f6e0f311676d25b8d64d9b731d07eefb
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Sep 3 14:56:36 2008 +1000
VariableDistance now works also in 3D
darcs-hash:20080903045636-d4795-c117a44854b4ae04d569e2b5720896ea27404aa6.gz
diff --git a/src/levelset.c b/src/levelset.c
index b0a8a0c..0142072 100644
--- a/src/levelset.c
+++ b/src/levelset.c
@@ -47,6 +47,15 @@ static void variable_distance_read (GtsObject ** o, GtsFile * fp)
"Distance to the interface defined by tracer",
fp->token->str, NULL);
gts_file_next_token (fp);
+
+ if (fp->type == '{') {
+ GtsFileVariable var[] = {
+ {GTS_INT, "stencil", TRUE},
+ {GTS_NONE}
+ };
+ var[0].data = &GFS_VARIABLE_DISTANCE (*o)->stencil;
+ gts_file_assign_variables (fp, var);
+ }
}
static void variable_distance_write (GtsObject * o, FILE * fp)
@@ -54,58 +63,55 @@ static void variable_distance_write (GtsObject * o, FILE * fp)
(* GTS_OBJECT_CLASS (gfs_variable_distance_class ())->parent_class->write) (o, fp);
fprintf (fp, " %s", GFS_VARIABLE_DISTANCE (o)->v->name);
+ if (GFS_VARIABLE_DISTANCE (o)->stencil)
+ fputs (" { stencil = 1 }", fp);
}
static gdouble vof_distance2 (FttCell * cell, GtsPoint * t, gpointer v)
{
- gdouble f = GFS_VARIABLE (cell, GFS_VARIABLE1 (v)->i);
+ gdouble f = GFS_VALUE (cell, GFS_VARIABLE1 (v));
if (GFS_IS_FULL (f))
return G_MAXDOUBLE;
if (!FTT_CELL_IS_LEAF (cell))
return ftt_cell_point_distance2_min (cell, t);
- else {
- FttVector p[2], m;
- guint n = gfs_vof_facet (cell, v, p, &m);
-
-#if FTT_3D
- g_assert_not_implemented ();
-#endif
-
- if (n != 2)
- return G_MAXDOUBLE;
-
- GtsPoint p1, p2;
- p1.x = p[0].x; p1.y = p[0].y; p1.z = 0.;
- p2.x = p[1].x; p2.y = p[1].y; p2.z = 0.;
- GtsSegment s;
- s.v1 = (GtsVertex *) &p1; s.v2 = (GtsVertex *) &p2;
- return gts_point_segment_distance2 (t, &s);
- }
+ else
+ return gfs_vof_facet_distance2 (cell, v, t);
}
-static void distance (FttCell * cell, gpointer * data)
+static void distance_for_stencil (FttCell * cell, gpointer * data)
{
GfsVariable * v = data[0];
GfsVariable * s2 = data[2];
- if (GFS_VARIABLE (cell, s2->i)) {
+ if (GFS_VALUE (cell, s2)) {
GfsVariableDistance * l = GFS_VARIABLE_DISTANCE (v);
GtsPoint p;
gdouble d2;
ftt_cell_pos (cell, (FttVector *) &p.x);
d2 = gfs_domain_cell_point_distance2 (v->domain, &p, vof_distance2, l->v, NULL);
- GFS_VARIABLE (cell, v->i) = GFS_VARIABLE (cell, l->v->i) > 0.5 ? sqrt (d2) : -sqrt (d2);
+ GFS_VALUE (cell, v) = GFS_VALUE (cell, l->v) > 0.5 ? sqrt (d2) : -sqrt (d2);
}
else
- GFS_VARIABLE (cell, v->i) = G_MAXDOUBLE;
+ GFS_VALUE (cell, v) = G_MAXDOUBLE;
+}
+
+static void distance (FttCell * cell, GfsVariable * v)
+{
+ GfsVariableDistance * l = GFS_VARIABLE_DISTANCE (v);
+ GtsPoint p;
+ gdouble d2;
+
+ ftt_cell_pos (cell, (FttVector *) &p.x);
+ d2 = gfs_domain_cell_point_distance2 (v->domain, &p, vof_distance2, l->v, NULL);
+ GFS_VALUE (cell, v) = GFS_VALUE (cell, l->v) > 0.5 ? sqrt (d2) : -sqrt (d2);
}
static void stencil_interpolate (FttCell * cell, gpointer * data)
{
GfsVariableDistance * v = data[0];
- gdouble f = GFS_VARIABLE (cell, v->v->i);
+ gdouble f = GFS_VALUE (cell, v->v);
if (!GFS_IS_FULL (f))
gfs_interpolate_stencil (cell, data[1]);
@@ -117,7 +123,7 @@ static void stencil_gradient (FttCell * cell, gpointer * data)
GfsVariable * s2 = data[2];
FttComponent c;
- if (GFS_VARIABLE (cell, s1->i))
+ if (GFS_VALUE (cell, s1))
for (c = 0; c < FTT_DIMENSION; c++)
gfs_center_gradient_stencil (cell, c, s2->i);
}
@@ -126,34 +132,41 @@ static void variable_distance_event_half (GfsEvent * event, GfsSimulation * sim)
{
GfsDomain * domain = GFS_DOMAIN (sim);
GfsVariableDistance * v = GFS_VARIABLE_DISTANCE (event);
- gpointer data[3], tmp;
gfs_domain_timer_start (domain, "distance");
- data[0] = v;
- data[1] = gfs_temporary_variable (domain);
- data[2] = gfs_temporary_variable (domain);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_reset, data[1]);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) stencil_interpolate, data);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) gfs_cell_reset, data[2]);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) stencil_gradient, data);
- tmp = data[1]; data[1] = data[2]; data[2] = tmp;
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) stencil_gradient, data);
-
- gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
- (FttCellTraverseFunc) v->v->fine_coarse, v->v);
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) distance, data);
+ if (v->stencil) { /* fixme: this "acceleration technique"
+ i.e. computing distance only in a band around
+ the interface seems to be slower than computing
+ the distance function everywhere! */
+ gpointer data[3], tmp;
+ data[0] = v;
+ data[1] = gfs_temporary_variable (domain);
+ data[2] = gfs_temporary_variable (domain);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, data[1]);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) stencil_interpolate, data);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) gfs_cell_reset, data[2]);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) stencil_gradient, data);
+ tmp = data[1]; data[1] = data[2]; data[2] = tmp;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) stencil_gradient, data);
+
+ gfs_domain_cell_traverse (domain, FTT_POST_ORDER, FTT_TRAVERSE_NON_LEAFS, -1,
+ (FttCellTraverseFunc) v->v->fine_coarse, v->v);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) distance_for_stencil, data);
+ gts_object_destroy (data[1]);
+ gts_object_destroy (data[2]);
+ }
+ else
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) distance, v);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, GFS_VARIABLE1 (event));
- gts_object_destroy (data[1]);
- gts_object_destroy (data[2]);
-
gfs_domain_timer_stop (domain, "distance");
}
diff --git a/src/levelset.h b/src/levelset.h
index 2e8fede..bb84d16 100644
--- a/src/levelset.h
+++ b/src/levelset.h
@@ -37,6 +37,7 @@ struct _GfsVariableDistance {
/*< public >*/
GfsVariable * v;
+ gboolean stencil;
};
#define GFS_VARIABLE_DISTANCE(obj) GTS_OBJECT_CAST (obj,\
diff --git a/src/vof.c b/src/vof.c
index 61cc6f6..d2570db 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -1656,6 +1656,71 @@ guint gfs_vof_facet (FttCell * cell,
}
/**
+ * gfs_vof_facet_distance2:
+ * @cell: a #FttCell.
+ * @t: a #GfsVariableTracerVOF.
+ * @p: a #GtsPoint.
+ *
+ * Returns: the square of the distance between point @p and the
+ * VOF-reconstructed interface facet defined by @t or %G_MAXDOUBLE if
+ * @cell does not contain an interface.
+ */
+gdouble gfs_vof_facet_distance2 (FttCell * cell,
+ GfsVariableTracerVOF * t,
+ GtsPoint * p)
+{
+ g_return_val_if_fail (cell != NULL, G_MAXDOUBLE);
+ g_return_val_if_fail (t != NULL, G_MAXDOUBLE);
+ g_return_val_if_fail (p != NULL, G_MAXDOUBLE);
+
+ if (GFS_IS_FULL (GFS_VALUE (cell, GFS_VARIABLE1 (t))))
+ return G_MAXDOUBLE;
+
+ FttVector q, m;
+ ftt_cell_pos (cell, &q);
+ gdouble h = ftt_cell_size (cell), lambda = 0., norm2 = 0.;
+ FttComponent c;
+ q.x -= h/2.; q.y -= h/2.; q.z -= h/2.;
+ /* compute position of closest point on VOF plane m*x + m*y + m*z = alpha */
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ (&m.x)[c] = GFS_VALUE (cell, t->m[c]);
+ lambda += (&m.x)[c]*((&p->x)[c] - (&q.x)[c])/h;
+ norm2 += (&m.x)[c]*(&m.x)[c];
+ }
+ gdouble alpha = GFS_VALUE (cell, t->alpha);
+ g_assert (norm2 > 0.);
+ lambda = (lambda - alpha)/norm2;
+
+ FttVector o;
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ (&o.x)[c] = ((&p->x)[c] - (&q.x)[c])/h - lambda*(&m.x)[c];
+ if ((&o.x)[c] <= 0. || (&o.x)[c] >= 1.) {
+ /* closest point on VOF plane is not within cell
+ return minimum distance from facet edges */
+ FttVector q[FTT_DIMENSION*(FTT_DIMENSION - 1) + 1];
+ gdouble dmin = G_MAXDOUBLE;
+ guint i, n = gfs_vof_facet (cell, t, q, &m);
+#if !FTT_2D
+ if (n > 2)
+ q[n++] = q[0];
+#endif
+ for (i = 0; i < n - 1; i++) {
+ GtsPoint p1, p2;
+ p1.x = q[i].x; p1.y = q[i].y; p1.z = q[i].z;
+ p2.x = q[i + 1].x; p2.y = q[i + 1].y; p2.z = q[i + 1].z;
+ GtsSegment s;
+ s.v1 = (GtsVertex *) &p1; s.v2 = (GtsVertex *) &p2;
+ gdouble d = gts_point_segment_distance2 (p, &s);
+ if (d < dmin)
+ dmin = d;
+ }
+ return dmin;
+ }
+ }
+ return h*h*lambda*lambda*norm2;
+}
+
+/**
* gfs_vof_center:
* @cell: a #FttCell.
* @t: a #GfsVariableTracerVOF.
diff --git a/src/vof.h b/src/vof.h
index aa33102..5141596 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -90,6 +90,9 @@ guint gfs_vof_facet (FttCell * cell,
GfsVariableTracerVOF * t,
FttVector * p,
FttVector * m);
+gdouble gfs_vof_facet_distance2 (FttCell * cell,
+ GfsVariableTracerVOF * t,
+ GtsPoint * p);
gdouble gfs_vof_center (FttCell * cell,
GfsVariableTracerVOF * t,
FttVector * p);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list