[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:55 UTC 2009
The following commit has been merged in the upstream branch:
commit 08f7f2137e0b79ae95de1b96dcc5fb173567e87e
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Dec 11 12:54:57 2006 +1100
VariableCurvature uses either height-function method or levelset
According to the arguments (i.e. the second argument specifies a
VariableTracer or a VariableDistance).
darcs-hash:20061211015457-d4795-b44adfe37699bd26fd0304edce008953046d7ddf.gz
diff --git a/src/levelset.c b/src/levelset.c
index 674194a..b0a8a0c 100644
--- a/src/levelset.c
+++ b/src/levelset.c
@@ -72,7 +72,8 @@ static gdouble vof_distance2 (FttCell * cell, GtsPoint * t, gpointer v)
g_assert_not_implemented ();
#endif
- g_assert (n == 2);
+ if (n != 2)
+ return G_MAXDOUBLE;
GtsPoint p1, p2;
p1.x = p[0].x; p1.y = p[0].y; p1.z = 0.;
diff --git a/src/poisson.c b/src/poisson.c
index 3a5c862..b13076d 100644
--- a/src/poisson.c
+++ b/src/poisson.c
@@ -441,7 +441,7 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
gdouble w2 = c2*(1. - c2);
gdouble k1 = GFS_VARIABLE (face->cell, kappa->i);
gdouble k2 = GFS_VARIABLE (face->neighbor, kappa->i);
-
+
if (w1 + w2 > 0.)
v *= (w1*k1 + w2*k2)/(w1 + w2);
else {
@@ -456,6 +456,7 @@ static void tension_coeff (FttCellFace * face, gpointer * data)
else
v = 1e6;
}
+ g_assert (v <= 1e6);
if (alpha)
v *= gfs_function_face_value (alpha, face);
diff --git a/src/tension.c b/src/tension.c
index cb6fa9e..635ecd6 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -22,6 +22,7 @@
#include "tension.h"
#include "vof.h"
+#include "levelset.h"
/* GfsSourceTensionGeneric: Object */
@@ -327,7 +328,7 @@ static void variable_curvature_read (GtsObject ** o, GtsFile * fp)
return;
if (fp->type != GTS_STRING) {
- gts_file_error (fp, "expecting a string (f)");
+ gts_file_error (fp, "expecting a string (fraction or distance)");
return;
}
domain = GFS_DOMAIN (gfs_object_simulation (*o));
@@ -335,12 +336,22 @@ static void variable_curvature_read (GtsObject ** o, GtsFile * fp)
gts_file_error (fp, "unknown variable `%s'", fp->token->str);
return;
}
- gts_file_next_token (fp);
-
if (GFS_VARIABLE1 (v)->description)
g_free (GFS_VARIABLE1 (v)->description);
- GFS_VARIABLE1 (v)->description = g_strjoin (" ", "Curvature of the interface defined by tracer",
- v->f->name, NULL);
+ if (GFS_IS_VARIABLE_TRACER (v->f))
+ GFS_VARIABLE1 (v)->description = g_strjoin (" ",
+ "Curvature of the interface defined by tracer",
+ v->f->name, NULL);
+ 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);
+ 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)
@@ -389,18 +400,9 @@ static void diffuse (FttCell * cell, DiffuseParms * p)
}
}
-static void variable_curvature_event_half (GfsEvent * event, GfsSimulation * sim)
+static void variable_curvature_diffuse (GfsEvent * event, GfsSimulation * sim)
{
GfsDomain * domain = GFS_DOMAIN (sim);
-
- gfs_domain_timer_start (domain, "variable_curvature");
-
- gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) 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));
-
DiffuseParms p;
p.v = GFS_VARIABLE1 (event);
p.tmp = gfs_temporary_variable (domain);
@@ -415,11 +417,123 @@ static void variable_curvature_event_half (GfsEvent * event, GfsSimulation * sim
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, p.v);
}
- gts_object_destroy (GTS_OBJECT (p.tmp));
+ gts_object_destroy (GTS_OBJECT (p.tmp));
+}
+
+static void variable_curvature_from_fraction (GfsEvent * event, GfsSimulation * sim)
+{
+ GfsDomain * domain = GFS_DOMAIN (sim);
+
+ gfs_domain_timer_start (domain, "variable_curvature");
+
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) 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);
+
+ gfs_domain_timer_stop (domain, "variable_curvature");
+}
+
+
+static void normal (FttCell * cell, gpointer * data)
+{
+ GfsVariable ** nv = data[0];
+ GfsVariable * d = GFS_VARIABLE_CURVATURE (data[1])->f;
+ GtsVector n = { 0., 0., 0. };
+ FttComponent c;
+
+ for (c = 0; c < FTT_DIMENSION; c++)
+ n[c] = gfs_center_gradient (cell, c, d->i);
+ gts_vector_normalize (n);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ GFS_VARIABLE (cell, nv[c]->i) = n[c];
+}
+
+static void distance_curvature (FttCell * cell, gpointer * data)
+{
+ GfsVariable ** nv = data[0];
+ gdouble kappa = 0.;
+ FttComponent c;
+
+ for (c = 0; c < FTT_DIMENSION; c++)
+ kappa += gfs_center_gradient (cell, c, nv[c]->i);
+ GFS_VARIABLE (cell, nv[FTT_DIMENSION]->i) = kappa/ftt_cell_size (cell);
+}
+
+static void interface_curvature (FttCell * cell, gpointer * data)
+{
+ GfsVariable * v = data[1];
+ GfsVariableCurvature * k = GFS_VARIABLE_CURVATURE (v);
+ gdouble f = GFS_VARIABLE (cell, GFS_VARIABLE_DISTANCE (k->f)->v->i);
+
+ if (GFS_IS_FULL (f))
+ GFS_VARIABLE (cell, v->i) = G_MAXDOUBLE;
+ else {
+ GfsVariable ** nv = data[0];
+ gdouble h = ftt_cell_size (cell)/2.;
+ FttCell * target = cell;
+ FttComponent c;
+ FttVector p;
+
+ ftt_cell_pos (cell, &p);
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ gdouble delta = GFS_VARIABLE (cell, k->f->i)*GFS_VARIABLE (cell, nv[c]->i);
+ (&p.x)[c] -= delta;
+ if (fabs (delta) > h)
+ target = NULL;
+ }
+ if (!target)
+ target = gfs_domain_locate (v->domain, p, -1);
+ GFS_VARIABLE (cell, v->i) = gfs_interpolate (target, p, nv[FTT_DIMENSION]);
+ }
+}
+
+static void variable_curvature_from_distance (GfsEvent * event, GfsSimulation * sim)
+{
+ GfsVariable * n[FTT_DIMENSION + 1];
+ GfsDomain * domain = GFS_DOMAIN (sim);
+ gpointer data[2];
+ FttComponent c;
+
+ gfs_domain_timer_start (domain, "variable_curvature");
+
+ for (c = 0; c < FTT_DIMENSION + 1; c++) {
+ n[c] = gfs_temporary_variable (domain);
+ gfs_variable_set_vector (n[c], c);
+ }
+ data[0] = n;
+ data[1] = event;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) normal, data);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, n[c]);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) distance_curvature, data);
+ gfs_domain_copy_bc (domain, FTT_TRAVERSE_LEAFS, -1,
+ GFS_VARIABLE1 (event), n[FTT_DIMENSION]);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) interface_curvature, data);
+ 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));
+ for (c = 0; c < FTT_DIMENSION + 1; c++)
+ gts_object_destroy (GTS_OBJECT (n[c]));
+
+ variable_curvature_diffuse (event, sim);
gfs_domain_timer_stop (domain, "variable_curvature");
}
+static void variable_curvature_event_half (GfsEvent * event, GfsSimulation * sim)
+{
+ if (GFS_IS_VARIABLE_TRACER (GFS_VARIABLE_CURVATURE (event)->f))
+ variable_curvature_from_fraction (event, sim);
+ else /* distance */
+ variable_curvature_from_distance (event, sim);
+}
+
static gboolean variable_curvature_event (GfsEvent * event, GfsSimulation * sim)
{
if ((* GFS_EVENT_CLASS (GTS_OBJECT_CLASS (gfs_variable_curvature_class ())->parent_class)->event)
@@ -566,6 +680,8 @@ static void variable_position_event_half (GfsEvent * event, GfsSimulation * sim)
(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_position");
}
diff --git a/src/vof.c b/src/vof.c
index fb63306..a6096d5 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -870,7 +870,7 @@ guint gfs_vof_facet (FttCell * cell, GfsVariable * v, FttVector * p, FttVector *
p[n].x = q.x + h*(0.5 - x); p[n].y = q.y - h/2.; p[n++].z = 0.;
}
}
- g_assert (n == 2);
+ g_assert (n <= 2);
#else /* 3D */
gdouble max = fabs (m->x);
FttComponent c = FTT_X;
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list