[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
s.popinet at niwa.co.nz
Fri May 15 02:51:35 UTC 2009
The following commit has been merged in the upstream branch:
commit 4811b577b1bfa71b73a3a734de65bd028f82630a
Author: Stephane Popinet <s.popinet at niwa.co.nz>
Date: Fri Feb 4 14:39:22 2005 +1100
Replaced "marker" surface tension implementation with tensorial-CSF formulation
darcs-hash:20050204033922-fbd8f-853654b6317c0aedf842a2ddd67b95c774c86470.gz
diff --git a/src/tension.c b/src/tension.c
index d0587de..413da68 100644
--- a/src/tension.c
+++ b/src/tension.c
@@ -21,186 +21,36 @@
#include <stdlib.h>
#include "tension.h"
-
-static void surface_add_box (GtsSurface * s,
- gdouble x1, gdouble y1, gdouble z1,
- gdouble x2, gdouble y2, gdouble z2)
-{
- GtsVertex * v0 = gts_vertex_new (s->vertex_class, x1, y1, z1);
- GtsVertex * v1 = gts_vertex_new (s->vertex_class, x1, y1, z2);
- GtsVertex * v2 = gts_vertex_new (s->vertex_class, x1, y2, z2);
- GtsVertex * v3 = gts_vertex_new (s->vertex_class, x1, y2, z1);
- GtsVertex * v4 = gts_vertex_new (s->vertex_class, x2, y1, z1);
- GtsVertex * v5 = gts_vertex_new (s->vertex_class, x2, y1, z2);
- GtsVertex * v6 = gts_vertex_new (s->vertex_class, x2, y2, z2);
- GtsVertex * v7 = gts_vertex_new (s->vertex_class, x2, y2, z1);
-
- GtsEdge * e1 = gts_edge_new (s->edge_class, v0, v1);
- GtsEdge * e2 = gts_edge_new (s->edge_class, v1, v2);
- GtsEdge * e3 = gts_edge_new (s->edge_class, v2, v3);
- GtsEdge * e4 = gts_edge_new (s->edge_class, v3, v0);
- GtsEdge * e5 = gts_edge_new (s->edge_class, v0, v2);
-
- GtsEdge * e6 = gts_edge_new (s->edge_class, v4, v5);
- GtsEdge * e7 = gts_edge_new (s->edge_class, v5, v6);
- GtsEdge * e8 = gts_edge_new (s->edge_class, v6, v7);
- GtsEdge * e9 = gts_edge_new (s->edge_class, v7, v4);
- GtsEdge * e10 = gts_edge_new (s->edge_class, v4, v6);
-
- GtsEdge * e11 = gts_edge_new (s->edge_class, v3, v7);
- GtsEdge * e12 = gts_edge_new (s->edge_class, v2, v6);
- GtsEdge * e13 = gts_edge_new (s->edge_class, v1, v5);
- GtsEdge * e14 = gts_edge_new (s->edge_class, v0, v4);
-
- GtsEdge * e15 = gts_edge_new (s->edge_class, v1, v6);
- GtsEdge * e16 = gts_edge_new (s->edge_class, v2, v7);
- GtsEdge * e17 = gts_edge_new (s->edge_class, v3, v4);
- GtsEdge * e18 = gts_edge_new (s->edge_class, v0, v5);
-
- GtsFaceClass * klass = gts_face_class ();
-
- gts_surface_add_face (s, gts_face_new (klass, e1, e2, e5));
- gts_surface_add_face (s, gts_face_new (klass, e5, e3, e4));
- gts_surface_add_face (s, gts_face_new (klass, e6, e10, e7));
- gts_surface_add_face (s, gts_face_new (klass, e10, e9, e8));
- gts_surface_add_face (s, gts_face_new (klass, e2, e15, e12));
- gts_surface_add_face (s, gts_face_new (klass, e15, e13, e7));
- gts_surface_add_face (s, gts_face_new (klass, e3, e16, e11));
- gts_surface_add_face (s, gts_face_new (klass, e16, e12, e8));
- gts_surface_add_face (s, gts_face_new (klass, e17, e14, e4));
- gts_surface_add_face (s, gts_face_new (klass, e17, e11, e9));
- gts_surface_add_face (s, gts_face_new (klass, e18, e13, e1));
- gts_surface_add_face (s, gts_face_new (klass, e18, e14, e6));
-}
-
-static GtsBBox * bbox_cell (GtsBBoxClass * klass, FttCell * cell)
-{
- FttVector p;
- gdouble size;
-
- ftt_cell_pos (cell, &p);
- size = ftt_cell_size (cell)/2.;
- return gts_bbox_new (klass, cell,
- p.x - size, p.y - size, p.z - size,
- p.x + size, p.y + size, p.z + size);
-}
-
-static void Knds (GtsTriangle * t, FttVector * te)
-{
- GtsVertex * v1, * v2, * v3;
- FttComponent c;
- gdouble a = gts_triangle_area (t);
-
- gts_triangle_vertices (t, &v1, &v2, &v3);
- for (c = 0; c < FTT_DIMENSION; c++)
- (&te->x)[c] += a*(GTS_VERTEX_NORMAL (v1)->n[c] +
- GTS_VERTEX_NORMAL (v2)->n[c] +
- GTS_VERTEX_NORMAL (v3)->n[c]);
-}
-
-static void cell_tension (FttCell * cell,
- GtsSurface * s,
- GNode * stree,
- FttVector * t)
-{
- GtsBBox * bbox;
-
- t->x = t->y = t->z = 0.;
- bbox = bbox_cell (gts_bbox_class (), cell);
- if (gts_bb_tree_is_overlapping (stree, bbox)) {
- GtsSurface * cs;
- GNode * ctree;
- GtsSurfaceInter * si;
-
- if (GFS_IS_MIXED (cell))
- g_assert_not_implemented ();
-
- cs = gts_surface_new (gts_surface_class (),
- gts_face_class (),
- gts_edge_class (),
- gts_vertex_class ());
- surface_add_box (cs,
- bbox->x1, bbox->y1, bbox->z1,
- bbox->x2, bbox->y2, bbox->z2);
- ctree = gts_bb_tree_surface (cs);
- si = gts_surface_inter_new (gts_surface_inter_class (),
- s, cs, stree, ctree, FALSE, FALSE);
- if (si->edges != NULL) {
- GtsSurface * inter = gts_surface_new (gts_surface_class (),
- gts_face_class (),
- gts_edge_class (),
- gts_vertex_class ());
- FttComponent c;
- gdouble h = ftt_cell_size (cell);
- gdouble vol = 6.*h*h*h;
-
- gts_surface_inter_boolean (si, inter, GTS_1_IN_2);
- gts_surface_foreach_face (inter, (GtsFunc) Knds, t);
- for (c = 0; c < FTT_DIMENSION; c++)
- (&t->x)[c] /= vol;
-#if 0
- {
- GtsVector p;
-
- gts_surface_center_of_area (inter, p);
- printf ("VECT 1 2 0 2 0 %g %g %g %g %g %g\n",
- p[0], p[1], p[2],
- p[0] - 3.*t->x,
- p[1] - 3.*t->y,
- p[2] - 3.*t->z);
- }
-#endif
- gts_object_destroy (GTS_OBJECT (inter));
- }
- gts_object_destroy (GTS_OBJECT (si));
- gts_bb_tree_destroy (ctree, TRUE);
- gts_object_destroy (GTS_OBJECT (cs));
- }
- gts_object_destroy (GTS_OBJECT (bbox));
-}
-
-static void foreach_cell_tension (FttCell * cell, gpointer * data)
-{
- FttVector t;
- GfsSimulation * sim = data[0];
- GfsSourceTension * s = data[1];
- FttComponent c;
-
- cell_tension (cell, sim->interface, sim->itree, &t);
- for (c = 0; c < FTT_DIMENSION; c++)
- GFS_VARIABLE (cell, s->t[c]->i) = s->sigma*(&t.x)[c];
-}
-
-static void vertex_normal (GtsVertex * v, GtsSurface * s)
-{
- gts_vertex_mean_curvature_normal (v, s, GTS_VERTEX_NORMAL (v)->n);
-}
+#include "vof.h"
/* GfsSourceTension: Object */
static void gfs_source_tension_read (GtsObject ** o, GtsFile * fp)
{
GfsSourceTension * s = GFS_SOURCE_TENSION (*o);
+ GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
- if (GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->read)
- (* GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->read)
- (o, fp);
+ (* GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->read) (o, fp);
if (fp->type == GTS_ERROR)
return;
- s->t[0] = gfs_domain_add_variable (GFS_DOMAIN (gfs_object_simulation (*o)),
- "_gfs_source_tension_x");
- if (s->t[0] == NULL) {
- gts_file_error (fp, "only one GfsSourceTension is allowed");
+ if (fp->type != GTS_STRING) {
+ gts_file_error (fp, "expecting a variable (C)");
return;
}
- s->t[1] = gfs_domain_add_variable (GFS_DOMAIN (gfs_object_simulation (*o)),
- "_gfs_source_tension_y");
- g_assert (s->t[1]);
+ if ((s->c = gfs_variable_from_name (domain->variables, fp->token->str)) == NULL) {
+ gts_file_error (fp, "unknown variable `%s'", fp->token->str);
+ return;
+ }
+ gts_file_next_token (fp);
+
+ if ((s->t[0] = gfs_variable_from_name (domain->variables, "_gfs_source_tension_x")) == NULL)
+ s->t[0] = gfs_domain_add_variable (domain, "_gfs_source_tension_x");
+ if ((s->t[1] = gfs_variable_from_name (domain->variables, "_gfs_source_tension_y")) == NULL)
+ s->t[1] = gfs_domain_add_variable (domain, "_gfs_source_tension_y");
#if (!FTT_2D)
- s->t[2] = gfs_domain_add_variable (GFS_DOMAIN (gfs_object_simulation (*o)),
- "_gfs_source_tension_z");
- g_assert (s->t[2]);
+ if ((s->t[2] = gfs_variable_from_name (domain->variables, "_gfs_source_tension_z")) == NULL)
+ s->t[2] = gfs_domain_add_variable (domain, "_gfs_source_tension_z");
#endif /* 3D */
if (fp->type != GTS_INT && fp->type != GTS_FLOAT) {
@@ -214,28 +64,50 @@ static void gfs_source_tension_read (GtsObject ** o, GtsFile * fp)
static void gfs_source_tension_write (GtsObject * o, FILE * fp)
{
- if (GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->write)
- (* GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->write)
- (o, fp);
- fprintf (fp, " %g", GFS_SOURCE_TENSION (o)->sigma);
+ (* GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class->write) (o, fp);
+ fprintf (fp, " %s %g", GFS_SOURCE_TENSION (o)->c->name, GFS_SOURCE_TENSION (o)->sigma);
+}
+
+static void foreach_cell_normal (FttCell * cell, GfsVariable * v)
+{
+ FttVector n;
+ gdouble nn = 0.;
+ FttComponent c;
+
+ for (c = 0; c < FTT_DIMENSION; c++) {
+ (&n.x)[c] = gfs_youngs_gradient (cell, c, v);
+ nn += (&n.x)[c]*(&n.x)[c];
+ }
+ nn = sqrt (nn + 1e-50);
+ GFS_STATE (cell)->g[0] = n.x*n.x/nn;
+ GFS_STATE (cell)->g[1] = n.y*n.y/nn;
+ GFS_STATE (cell)->div = n.x*n.y/nn;
+}
+
+static void foreach_cell_tension (FttCell * cell, GfsSourceTension * s)
+{
+ gdouble h = ftt_cell_size (cell);
+ gdouble sigh = s->sigma/(h*h*h);
+
+ GFS_VARIABLE (cell, s->t[0]->i) = sigh*(gfs_youngs_gradient (cell, FTT_X, gfs_gy) -
+ gfs_youngs_gradient (cell, FTT_Y, gfs_div));
+ GFS_VARIABLE (cell, s->t[1]->i) = sigh*(gfs_youngs_gradient (cell, FTT_Y, gfs_gx) -
+ gfs_youngs_gradient (cell, FTT_X, gfs_div));
}
static gboolean gfs_source_tension_event (GfsEvent * event,
GfsSimulation * sim)
{
- if ((* GFS_EVENT_CLASS (GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class)->event) (event, sim)) {
- if (sim->interface) {
- gpointer data[2];
-
- gts_surface_foreach_vertex (sim->interface,
- (GtsFunc) vertex_normal, sim->interface);
- data[0] = sim;
- data[1] = event;
- gfs_domain_cell_traverse (GFS_DOMAIN (sim),
- FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) foreach_cell_tension,
- data);
- }
+ if ((* GFS_EVENT_CLASS (GTS_OBJECT_CLASS (gfs_source_tension_class ())->parent_class)->event)
+ (event, sim)) {
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim),
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) foreach_cell_normal,
+ GFS_SOURCE_TENSION (event)->c);
+ /* fixme: boundary conditions for normal */
+ gfs_domain_cell_traverse (GFS_DOMAIN (sim),
+ FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) foreach_cell_tension, event);
return TRUE;
}
return FALSE;
diff --git a/src/tension.h b/src/tension.h
index 306d90e..3fecd30 100644
--- a/src/tension.h
+++ b/src/tension.h
@@ -35,7 +35,7 @@ struct _GfsSourceTension {
GfsSourceVector parent;
/*< public >*/
- GfsVariable * t[FTT_DIMENSION];
+ GfsVariable * c, * t[FTT_DIMENSION];
gdouble sigma;
};
diff --git a/src/vof.c b/src/vof.c
index 8d21d94..8a83100 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -386,8 +386,20 @@ static gdouble avg (FttCell * cell, FttDirection d, GfsVariable * v)
return v1/(2*(FTT_DIMENSION - 1));
}
-static gdouble gfs_youngs_gradient (FttCell * cell, FttComponent c, GfsVariable * v)
+/**
+ * gfs_youngs_gradient:
+ * @cell: a #FttCell.
+ * @c: a component.
+ * @v: a #GfsVariable.
+ *
+ * Returns: the Youngs-averaged gradient of variable @v in direction
+ * @c normalised by the size of @cell.
+ */
+gdouble gfs_youngs_gradient (FttCell * cell, FttComponent c, GfsVariable * v)
{
+ g_return_val_if_fail (cell != NULL, 0.);
+ g_return_val_if_fail (v != NULL, 0.);
+
return avg (cell, 2*c, v) - avg (cell, 2*c + 1, v);
}
diff --git a/src/vof.h b/src/vof.h
index 40ee567..5d5d71e 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -49,6 +49,9 @@ void gfs_plane_center (FttVector * m,
gdouble a,
FttVector * p);
#endif /* 3D */
+gdouble gfs_youngs_gradient (FttCell * cell,
+ FttComponent c,
+ GfsVariable * v);
void gfs_cell_vof_advection (FttCell * cell,
FttComponent c,
GfsAdvectionParams * par);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list