[SCM] Gerris Flow Solver branch, upstream, updated. e8f73a07832050124d2b8bf6c6f35b33180e65a8
Stephane Popinet
popinet at users.sf.net
Tue Nov 24 12:24:19 UTC 2009
The following commit has been merged in the upstream branch:
commit aa8f66bb24574cde7cd1a15b14645f6c8a7217c9
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Jul 1 14:29:53 2009 +1000
GfsInitFraction should now work also for parallel simulations
It used to fail for domains containing disconnected GfsBoxes (which can
happen for some domain decompositions in parallel).
darcs-hash:20090701042953-d4795-defe3fdda69ceba41dcdd5c630d3ef7cb948bda7.gz
diff --git a/src/solid.c b/src/solid.c
index e215375..225ae3b 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -814,14 +814,16 @@ static void solid_fractions_from_children (FttCell * cell, InitSolidParams * p)
GFS_VALUE (cell, p->status) = GFS_STATUS_SOLID;
else {
gfs_cell_init_solid_fractions_from_children (cell);
- if (p->destroy_solid)
- GFS_VALUE (cell, p->status) = GFS_STATUS_UNDEFINED;
- else if (!GFS_IS_MIXED (cell)) {
+ GFS_VALUE (cell, p->status) = GFS_STATUS_UNDEFINED;
+ if (!p->destroy_solid && !GFS_IS_MIXED (cell)) {
ftt_cell_children (cell, &child);
- GFS_VALUE (cell, p->status) = GFS_STATUS_SOLID;
for (i = 0; i < FTT_CELLS; i++)
- if (child.c[i] && GFS_VALUE (child.c[i], p->status) == GFS_STATUS_FLUID)
- GFS_VALUE (cell, p->status) = GFS_STATUS_FLUID;
+ if (child.c[i]) {
+ if (GFS_VALUE (cell, p->status) == GFS_STATUS_UNDEFINED)
+ GFS_VALUE (cell, p->status) = GFS_VALUE (child.c[i], p->status);
+ else
+ g_assert (GFS_VALUE (cell, p->status) == GFS_VALUE (child.c[i], p->status));
+ }
}
}
}
@@ -1169,21 +1171,14 @@ static void save_solid (FttCell * cell, GfsVariable * c)
static void restore_solid (FttCell * cell, gpointer * data)
{
- GfsVariable * c = data[0];
- gboolean * not_cut = data[1];
- GfsVariable * status = data[2];
+ GfsVariable * status = data[0];
+ GfsVariable * c = data[1];
GfsSolidVector * solid = GFS_STATE (cell)->solid;
GFS_STATE (cell)->solid = GFS_DOUBLE_TO_POINTER (GFS_VALUE (cell, c));
if (solid) {
GFS_VALUE (cell, c) = solid->a;
g_free (solid);
- *not_cut = FALSE;
- }
- else if (GFS_VALUE (cell, status) == GFS_STATUS_UNDEFINED) {
- /* fixme: this can fail for non-contiguous domains (e.g. non-connected GfsBoxes) */
- g_assert (*not_cut);
- GFS_VALUE (cell, c) = 0.;
}
else {
g_assert (GFS_VALUE (cell, status) == GFS_STATUS_SOLID ||
@@ -1192,6 +1187,28 @@ static void restore_solid (FttCell * cell, gpointer * data)
}
}
+static void set_status (FttCell * cell, gpointer * data)
+{
+ GfsVariable * status = data[0];
+ gdouble * val = data[2];
+ GFS_VALUE (cell, status) = *val;
+}
+
+static void check_status (GfsBox * box, gpointer * data)
+{
+ GfsVariable * status = data[0];
+ if (!GFS_IS_MIXED (box->root) && GFS_VALUE (box->root, status) == GFS_STATUS_UNDEFINED) {
+ GfsGenericSurface * s = data[1];
+ FttVector pos;
+ ftt_cell_pos (box->root, &pos);
+ gdouble val = gfs_surface_point_is_inside (s, &pos) > 0 ?
+ GFS_STATUS_FLUID : GFS_STATUS_SOLID;
+ data[2] = &val;
+ ftt_cell_traverse (box->root, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) set_status, data);
+ }
+}
+
/**
* gfs_domain_init_fraction:
* @domain: a #GfsDomain.
@@ -1205,7 +1222,6 @@ void gfs_domain_init_fraction (GfsDomain * domain,
GfsGenericSurface * s,
GfsVariable * c)
{
- gboolean not_cut = TRUE;
gpointer data[3];
GfsVariable * status;
@@ -1222,9 +1238,10 @@ void gfs_domain_init_fraction (GfsDomain * domain,
GSList * l = g_slist_prepend (NULL, &tmp);
gfs_domain_init_solid_fractions (domain, l, FALSE, NULL, NULL, status);
g_slist_free (l);
- data[0] = c;
- data[1] = ¬_cut;
- data[2] = status;
+ data[0] = status;
+ data[1] = s;
+ gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) check_status, data);
+ data[1] = c;
gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) restore_solid, data);
gfs_domain_bc (domain, FTT_TRAVERSE_LEAFS, -1, c);
diff --git a/src/surface.c b/src/surface.c
index b67f635..5c9dbfa 100644
--- a/src/surface.c
+++ b/src/surface.c
@@ -240,6 +240,24 @@ void gfs_surface_segment_normal (GfsGenericSurface * s,
(* GFS_GENERIC_SURFACE_CLASS (GTS_OBJECT (s)->klass)->segment_normal) (s, cell, I, n);
}
+/**
+ * gfs_surface_point_is_inside:
+ * @s: a #GfsGenericSurface.
+ * @p: a #FttVector.
+ *
+ * Returns: 1 if @p is inside @s, 0 if @p lies on the boundary of @s,
+ * -1 otherwise.
+ */
+gint gfs_surface_point_is_inside (GfsGenericSurface * s,
+ FttVector * p)
+{
+ g_return_val_if_fail (s != NULL, 0);
+ g_return_val_if_fail (p != NULL, 0);
+
+ g_assert (GFS_GENERIC_SURFACE_CLASS (GTS_OBJECT (s)->klass)->point_is_inside);
+ return (* GFS_GENERIC_SURFACE_CLASS (GTS_OBJECT (s)->klass)->point_is_inside) (s, p);
+}
+
/* GfsSurface: Object */
static void check_solid_surface (GtsSurface * s,
@@ -744,6 +762,37 @@ static void surface_segment_normal (GfsGenericSurface * s1,
}
}
+static void add_tetrahedron_volume (GtsTriangle * t, gpointer * data)
+{
+ GtsPoint * p = data[0];
+ GtsVertex * v1, * v2, * v3;
+ gdouble * vol = data[1];
+ gts_triangle_vertices (t, &v1, &v2, &v3);
+ *vol += gts_point_orientation_3d (GTS_POINT (v1), GTS_POINT (v2), GTS_POINT (v3), p);
+}
+
+static gint point_is_inside_surface (GfsGenericSurface * s1, FttVector * p)
+{
+ GfsSurface * s = GFS_SURFACE (s1);
+ GtsPoint q;
+ q.x = p->x; q.y = p->y; q.z = p->z;
+
+ if (s->s) {
+ gdouble vol = 0.;
+ gpointer data[2];
+ data[0] = &q;
+ data[1] = &vol;
+ gts_surface_foreach_face (s->s, (GtsFunc) add_tetrahedron_volume, data);
+ fprintf (stderr, "vol: %g\n", vol);
+ return vol > 0. ? -1 : 1;
+ }
+ else {
+ gdouble v = gfs_surface_implicit_value (s, q);
+ return v == 0. ? 0 : v < 0. ? -1 : 1;
+ }
+ return 0;
+}
+
static void gfs_surface_class_init (GtsObjectClass * klass)
{
klass->read = surface_read;
@@ -753,6 +802,7 @@ static void gfs_surface_class_init (GtsObjectClass * klass)
GFS_GENERIC_SURFACE_CLASS (klass)->cell_is_cut = cell_is_cut;
GFS_GENERIC_SURFACE_CLASS (klass)->segment_intersection = surface_segment_intersection;
GFS_GENERIC_SURFACE_CLASS (klass)->segment_normal = surface_segment_normal;
+ GFS_GENERIC_SURFACE_CLASS (klass)->point_is_inside = point_is_inside_surface;
}
static void gfs_surface_init (GfsSurface * s)
diff --git a/src/surface.h b/src/surface.h
index ec7573a..430ae11 100644
--- a/src/surface.h
+++ b/src/surface.h
@@ -52,6 +52,8 @@ struct _GfsGenericSurfaceClass {
FttCell * cell,
GfsSegment * I,
GtsVector n);
+ gint (* point_is_inside) (GfsGenericSurface * s,
+ FttVector * p);
};
#define GFS_GENERIC_SURFACE(obj) GTS_OBJECT_CAST (obj,\
@@ -71,6 +73,8 @@ void gfs_surface_segment_normal (GfsGenericSurface * s,
FttCell * cell,
GfsSegment * I,
GtsVector n);
+gint gfs_surface_point_is_inside (GfsGenericSurface * s,
+ FttVector * p);
GfsGenericSurface * gfs_cell_is_cut (FttCell * cell,
GfsGenericSurface * s,
gboolean flatten,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list