[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:49 UTC 2009
The following commit has been merged in the upstream branch:
commit 3104081504ced15d50f85faac7562b35e1216c11
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Dec 6 11:20:18 2006 +1100
Functions to compute cell/plane intersections have been moved from GfsView
They are used to compute 3D VOF facets.
Note also that the interface of gfs_vof_facet() has changed.
darcs-hash:20061206002018-d4795-c315a6758384b3a9019496e7d4c4ba449a199193.gz
diff --git a/src/Makefile.am b/src/Makefile.am
index 0977427..5ded918 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -52,6 +52,7 @@ GFS_HDS = \
utils.h \
ocean.h \
levelset.h \
+ isocube.h \
$(MPI_HDS)
pkginclude_HEADERS = \
@@ -90,7 +91,7 @@ SRC = \
libgfs3D_la_LDFLAGS = $(NO_UNDEFINED)\
-version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE)\
-release $(LT_RELEASE) -export-dynamic
-libgfs3D_la_SOURCES = $(SRC) isocube.h
+libgfs3D_la_SOURCES = $(SRC)
libgfs2D_la_LDFLAGS = $(NO_UNDEFINED)\
-version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE)\
@@ -101,7 +102,7 @@ libgfs2D_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D=1
libgfs2D3_la_LDFLAGS = $(NO_UNDEFINED)\
-version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE)\
-release $(LT_RELEASE) -export-dynamic
-libgfs2D3_la_SOURCES = $(SRC) isocube.h
+libgfs2D3_la_SOURCES = $(SRC)
libgfs2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
CLEANFILES = $(BUILT_SOURCES)
@@ -109,8 +110,7 @@ CLEANFILES = $(BUILT_SOURCES)
EXTRA_DIST = \
mpi_boundary.c \
mpi_boundary.h \
- ftt_internal.c \
- isocube.h
+ ftt_internal.c
bin_PROGRAMS = gerris2D gerris3D gerris2D3
diff --git a/src/graphic.c b/src/graphic.c
index b6b5bc6..3b03cf3 100644
--- a/src/graphic.c
+++ b/src/graphic.c
@@ -26,6 +26,10 @@
#include "solid.h"
#include "variable.h"
+#if !FTT_2D
+# include "isocube.h"
+#endif /* 3D */
+
typedef struct {
GPtrArray * colors;
gboolean reversed;
@@ -1979,3 +1983,158 @@ void gfs_draw_streamline (GfsDomain * domain,
gfs_streamline_draw (path, fp);
gfs_streamline_destroy (path);
}
+
+#if !FTT_2D
+
+static gdouble point_orientation (FttVector p[3], FttVector * c)
+{
+ gdouble adx, bdx, cdx;
+ gdouble ady, bdy, cdy;
+ gdouble adz, bdz, cdz;
+
+ adx = p[0].x - c->x;
+ bdx = p[1].x - c->x;
+ cdx = p[2].x - c->x;
+ ady = p[0].y - c->y;
+ bdy = p[1].y - c->y;
+ cdy = p[2].y - c->y;
+ adz = p[0].z - c->z;
+ bdz = p[1].z - c->z;
+ cdz = p[2].z - c->z;
+
+ return (adx * (bdy * cdz - bdz * cdy) +
+ bdx * (cdy * adz - cdz * ady) +
+ cdx * (ady * bdz - adz * bdy));
+}
+
+/**
+ * gfs_plane_cuts_cell:
+ * @plane: three points belonging to the plane.
+ * @cell: a #FttCell.
+ *
+ * Returns: %TRUE if @plane cuts @cell, %FALSE otherwise.
+ */
+gboolean gfs_plane_cuts_cell (FttVector plane[3], FttCell * cell)
+{
+ FttVector o;
+ gdouble h = ftt_cell_size (cell)*SLIGHTLY_LARGER;
+ guint i;
+
+ g_return_val_if_fail (cell != NULL, FALSE);
+
+ ftt_cell_pos (cell, &o);
+ o.x -= h/2.; o.y -= h/2.; o.z -= h/2.;
+ for (i = 0; i < 12; i++) {
+ FttVector e, d;
+ gdouble a, b;
+ d.x = o.x + h*edge[i][0].x; d.y = o.y + h*edge[i][0].y; d.z = o.z + h*edge[i][0].z;
+ e.x = o.x + h*edge[i][1].x; e.y = o.y + h*edge[i][1].y; e.z = o.z + h*edge[i][1].z;
+ a = point_orientation (plane, &e);
+ b = point_orientation (plane, &d);
+ if ((a <= 0. && b > 0.) || (a >= 0. && b < 0.))
+ return TRUE;
+ }
+ return FALSE;
+}
+
+static void cube_plane_intersection (FttCell * cell,
+ FttVector * O,
+ FttVector * n,
+ FttVector p[12],
+ gint orient[12],
+ GfsVariable * var,
+ gdouble v[12],
+ gint max_level)
+{
+ FttVector o;
+ gdouble h = ftt_cell_size (cell)*SLIGHTLY_LARGER, vc[8];
+ guint i;
+
+ if (var)
+ for (i = 0; i < 8; i++)
+ vc[i] = G_MAXDOUBLE;
+
+ ftt_cell_pos (cell, &o);
+ o.x -= h/2.; o.y -= h/2.; o.z -= h/2.;
+ for (i = 0; i < 12; i++) {
+ FttVector e, d;
+ d.x = o.x + h*edge[i][0].x; d.y = o.y + h*edge[i][0].y; d.z = o.z + h*edge[i][0].z;
+ e.x = o.x + h*edge[i][1].x; e.y = o.y + h*edge[i][1].y; e.z = o.z + h*edge[i][1].z;
+ gdouble den = n->x*(e.x - d.x) + n->y*(e.y - d.y) + n->z*(e.z - d.z);
+ orient[i] = -1;
+ if (fabs (den) > 1e-10) {
+ gdouble t = (n->x*(O->x - d.x) + n->y*(O->y - d.y) + n->z*(O->z - d.z))/den;
+ if (t >= 0. && t < 1.) {
+ p[i].x = d.x + t*(e.x - d.x); p[i].y = d.y + t*(e.y - d.y); p[i].z = d.z + t*(e.z - d.z);
+ orient[i] = (n->x*(e.x - O->x) + n->y*(e.y - O->y) + n->z*(e.z - O->z) > 0.);
+ if (var) {
+ guint j = edge1[i][0];
+ if (vc[j] == G_MAXDOUBLE)
+ vc[j] = gfs_cell_corner_value (cell, corner[j], var, max_level);
+ d.z = vc[j];
+ j = edge1[i][1];
+ if (vc[j] == G_MAXDOUBLE)
+ vc[j] = gfs_cell_corner_value (cell, corner[j], var, max_level);
+ e.z = vc[j];
+ v[i] = d.z + t*(e.z - d.z);
+ }
+ }
+ }
+ }
+}
+
+/**
+ * gfs_cut_cube_vertices:
+ * @cell: a #FttCell.
+ * @maxlevel: the maximum level to consider (or -1).
+ * @p: a point on the plane.
+ * @n: the normal to the plane.
+ * @v: where to return the vertices coordinates.
+ * @d: where to return the direction.
+ * @var: a #GfsVariable or %NULL.
+ * @val: where to return the values of @var or %NULL.
+ *
+ * Fills @v, @d and @val with the coordinates/values of the vertices,
+ * intersections of @cell with the plane defined by @p and @n.
+ *
+ * The vertices are ordered consistently to define a consistent,
+ * oriented polygon.
+ *
+ * Returns: the number of vertices (0 if the plane does not cut the cell).
+ */
+guint gfs_cut_cube_vertices (FttCell * cell, gint maxlevel,
+ FttVector * p, FttVector * n,
+ FttVector v[12], FttDirection d[12],
+ GfsVariable * var,
+ gdouble val[12])
+{
+ FttVector a[12];
+ gdouble vv[12];
+ gint orient[12];
+ guint i;
+
+ g_return_val_if_fail (cell != NULL, 0);
+ g_return_val_if_fail (p != NULL, 0);
+ g_return_val_if_fail (n != NULL, 0);
+ g_return_val_if_fail ((var == NULL && val == NULL) || (var != NULL && val != NULL), 0);
+
+ cube_plane_intersection (cell, p, n, a, orient, var, vv, maxlevel);
+ for (i = 0; i < 12; i++) {
+ guint nv = 0, e = i;
+ while (orient[e] >= 0) {
+ guint m = 0, * ne = connect[e][orient[e]];
+ d[nv] = ne[3];
+ if (var)
+ val[nv] = vv[e];
+ v[nv++] = a[e];
+ orient[e] = -1;
+ while (m < 3 && orient[e] < 0)
+ e = ne[m++];
+ }
+ if (nv > 2)
+ return nv;
+ }
+ return 0;
+}
+
+#endif /* 3D */
diff --git a/src/graphic.h b/src/graphic.h
index 1380bae..c7b7cca 100644
--- a/src/graphic.h
+++ b/src/graphic.h
@@ -123,6 +123,14 @@ void gfs_draw_stream_ribbon (GfsDomain * domain,
void gfs_draw_streamline (GfsDomain * domain,
FttVector p,
FILE * fp);
+gboolean gfs_plane_cuts_cell (FttVector plane[3],
+ FttCell * cell);
+guint gfs_cut_cube_vertices (FttCell * cell,
+ gint maxlevel,
+ FttVector * p, FttVector * n,
+ FttVector v[12], FttDirection d[12],
+ GfsVariable * var,
+ gdouble val[12]);
#ifdef __cplusplus
}
diff --git a/src/isocube.h b/src/isocube.h
index 2c3819f..3abe95c 100644
--- a/src/isocube.h
+++ b/src/isocube.h
@@ -44,17 +44,55 @@ static guint face_v[6][4] = {
{1,5,7,3},/* front */
{0,4,6,2} /* back */
};
-static guint connect[12][2][3] = {
- {{9, 1, 8}, {4, 3, 7}}, /* 0 */
- {{6, 2, 5}, {8, 0, 9}}, /* 1 */
- {{10, 3, 11}, {5, 1, 6}}, /* 2 */
- {{7, 0, 4}, {11, 2, 10}}, /* 3 */
- {{3, 7, 0}, {8, 5, 11}}, /* 4 */
- {{11, 4, 8}, {1, 6, 2}}, /* 5 */
- {{2, 5, 1}, {9, 7, 10}}, /* 6 */
- {{10, 6, 9}, {0, 4, 3}}, /* 7 */
- {{5, 11, 4}, {0, 9, 1}}, /* 8 */
- {{1, 8, 0}, {7, 10, 6}}, /* 9 */
- {{6, 9, 7}, {3, 11, 2}}, /* 10 */
- {{2, 10, 3}, {4, 8, 5}} /* 11 */
+/* first index is the edge number, second index is the edge orientation
+ (0 or 1), third index are the edges which this edge may connect to
+ in order and the corresponding face direction */
+static guint connect[12][2][4] = {
+ {{9, 1, 8, FTT_BOTTOM}, {4, 3, 7, FTT_BACK}}, /* 0 */
+ {{6, 2, 5, FTT_FRONT}, {8, 0, 9, FTT_BOTTOM}}, /* 1 */
+ {{10, 3, 11, FTT_TOP}, {5, 1, 6, FTT_FRONT}}, /* 2 */
+ {{7, 0, 4, FTT_BACK}, {11, 2, 10, FTT_TOP}}, /* 3 */
+ {{3, 7, 0, FTT_BACK}, {8, 5, 11, FTT_LEFT}}, /* 4 */
+ {{11, 4, 8, FTT_LEFT}, {1, 6, 2, FTT_FRONT}}, /* 5 */
+ {{2, 5, 1, FTT_FRONT}, {9, 7, 10, FTT_RIGHT}}, /* 6 */
+ {{10, 6, 9, FTT_RIGHT}, {0, 4, 3, FTT_BACK}}, /* 7 */
+ {{5, 11, 4, FTT_LEFT}, {0, 9, 1, FTT_BOTTOM}}, /* 8 */
+ {{1, 8, 0, FTT_BOTTOM}, {7, 10, 6, FTT_RIGHT}}, /* 9 */
+ {{6, 9, 7, FTT_RIGHT}, {3, 11, 2, FTT_TOP}}, /* 10 */
+ {{2, 10, 3, FTT_TOP}, {4, 8, 5, FTT_LEFT}} /* 11 */
};
+static FttVector edge[12][2] = {
+ {{0.,0.,0.},{1.,0.,0.}},{{0.,0.,1.},{1.,0.,1.}},{{0.,1.,1.},{1.,1.,1.}},{{0.,1.,0.},{1.,1.,0.}},
+ {{0.,0.,0.},{0.,1.,0.}},{{0.,0.,1.},{0.,1.,1.}},{{1.,0.,1.},{1.,1.,1.}},{{1.,0.,0.},{1.,1.,0.}},
+ {{0.,0.,0.},{0.,0.,1.}},{{1.,0.,0.},{1.,0.,1.}},{{1.,1.,0.},{1.,1.,1.}},{{0.,1.,0.},{0.,1.,1.}}
+};
+static FttDirection corner[8][3] = {
+ {FTT_LEFT, FTT_BOTTOM, FTT_BACK},
+ {FTT_LEFT, FTT_BOTTOM, FTT_FRONT},
+ {FTT_LEFT, FTT_TOP, FTT_BACK},
+ {FTT_LEFT, FTT_TOP, FTT_FRONT},
+ {FTT_RIGHT, FTT_BOTTOM, FTT_BACK},
+ {FTT_RIGHT, FTT_BOTTOM, FTT_FRONT},
+ {FTT_RIGHT, FTT_TOP, FTT_BACK},
+ {FTT_RIGHT, FTT_TOP, FTT_FRONT}
+};
+static guint connectv[12][2][4] = {
+ {{4, 5, 1, 0}, {0, 2, 6, 4}}, /* 0 */
+ {{5, 7, 3, 1}, {1, 0, 4, 5}}, /* 1 */
+ {{7, 6, 2, 3}, {3, 1, 5, 7}}, /* 2 */
+ {{6, 4, 0, 2}, {2, 3, 7, 6}}, /* 3 */
+ {{2, 6, 4, 0}, {0, 1, 3, 2}}, /* 4 */
+ {{3, 2, 0, 1}, {1, 5, 7, 3}}, /* 5 */
+ {{7, 3, 1, 5}, {5, 4, 6, 7}}, /* 6 */
+ {{6, 7, 5, 4}, {4, 0, 2, 6}}, /* 7 */
+ {{1, 3, 2, 0}, {0, 4, 5, 1}}, /* 8 */
+ {{5, 1, 0, 4}, {4, 6, 7, 5}}, /* 9 */
+ {{7, 5, 4, 6}, {6, 2, 3, 7}}, /* 10 */
+ {{3, 7, 6, 2}, {2, 0, 1, 3}} /* 11 */
+};
+static FttVector cvertex[8] = {
+ {0., 0., 0.}, {0., 0., 1.}, {0., 1., 0.}, {0., 1., 1.},
+ {1., 0., 0.}, {1., 0., 1.}, {1., 1., 0.}, {1., 1., 1.}
+};
+
+#define SLIGHTLY_LARGER 1.001
diff --git a/src/levelset.c b/src/levelset.c
index 641fb9a..98097fe 100644
--- a/src/levelset.c
+++ b/src/levelset.c
@@ -60,25 +60,21 @@ static gdouble vof_distance2 (FttCell * cell, GtsPoint * t, gpointer v)
if (!FTT_CELL_IS_LEAF (cell))
return ftt_cell_point_distance2_min (cell, t);
else {
- gdouble h = ftt_cell_size (cell), d2;
- GSList * l = gfs_vof_facet (cell, v);
- GtsPoint * p1 = l->data, * p2 = l->next->data;
- GtsSegment s;
- FttVector p;
+ FttVector p[2], m;
+ guint n = gfs_vof_facet (cell, v, p, &m);
#if FTT_3D
- g_assert_not_implemented ();
+ g_assert_not_implemented ();
#endif
- ftt_cell_pos (cell, &p);
- p1->x = p.x + h*p1->x; p1->y = p.y + h*p1->y;
- p2->x = p.x + h*p2->x; p2->y = p.y + h*p2->y;
- s.v1 = (GtsVertex *) p1; s.v2 = (GtsVertex *) p2;
- d2 = gts_point_segment_distance2 (t, &s);
- gts_object_destroy (GTS_OBJECT (p1));
- gts_object_destroy (GTS_OBJECT (p2));
- g_slist_free (l);
- return d2;
+ g_assert (n == 2);
+
+ 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);
}
}
diff --git a/src/vof.c b/src/vof.c
index 1b7ec81..e82e7bd 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -20,6 +20,7 @@
#include <math.h>
#include "vof.h"
#include "variable.h"
+#include "graphic.h"
#define THRESHOLD(c) {if ((c) < 0.) c = 0.; else if ((c) > 1.) c = 1.;}
@@ -816,51 +817,78 @@ gboolean gfs_vof_plane (FttCell * cell, GfsVariable * v,
* gfs_vof_facet:
* @cell: a #FttCell.
* @v: a #GfsVariable.
+ * @p: a #FttVector array (of size 2 in 2D and 6 in 3D)
+ * @m: a #FttVector.
+ *
+ * Fills @p with the coordinates of points defining the
+ * VOF-reconstructed interface facet defined by @v.
*
- * Returns: a list of newly allocated #GtsPoint defining the
- * VOF-reconstructed interface facet defined by @v, or %NULL if @cell
- * is not cut by the interface.
+ * Fills @m with the normal to the interface.
+ *
+ * Returns: the number of points defining the facet.
*/
-GSList * gfs_vof_facet (FttCell * cell, GfsVariable * v)
+guint gfs_vof_facet (FttCell * cell, GfsVariable * v, FttVector * p, FttVector * m)
{
- FttVector m;
gdouble alpha;
- g_return_val_if_fail (cell != NULL, NULL);
- g_return_val_if_fail (v != NULL, NULL);
+ g_return_val_if_fail (cell != NULL, 0);
+ g_return_val_if_fail (v != NULL, 0);
+ g_return_val_if_fail (p != NULL, 0);
+ g_return_val_if_fail (m != NULL, 0);
-#if FTT_3D
- g_assert_not_implemented ();
-#endif
-
- if (!gfs_vof_plane (cell, v, &m, &alpha))
- return NULL;
+ if (!gfs_vof_plane (cell, v, m, &alpha))
+ return 0;
else {
- GSList * l = NULL;
+ guint n = 0;
+ FttVector q;
+ ftt_cell_pos (cell, &q);
+ gdouble h = ftt_cell_size (cell);
+#if FTT_2D
gdouble x, y;
- if (m.y != 0.) {
- y = alpha/m.y;
- if (y >= 0. && y <= 1.)
- l = g_slist_prepend (l, gts_point_new (gts_point_class (), 0.5, 0.5 - y, 0.));
+ if (m->y != 0.) {
+ y = alpha/m->y;
+ if (y >= 0. && y <= 1.) {
+ p[n].x = q.x + h/2.; p[n].y = q.y + h*(0.5 - y); p[n++].z = 0.;
+ }
}
- if (m.x != 0.) {
- x = alpha/m.x;
- if (x >= 0. && x <= 1.)
- l = g_slist_prepend (l, gts_point_new (gts_point_class (), 0.5 - x, 0.5, 0.));
+ if (m->x != 0.) {
+ x = alpha/m->x;
+ if (x >= 0. && x <= 1.) {
+ p[n].x = q.x + h*(0.5 - x); p[n].y = q.y + h/2.; p[n++].z = 0.;
+ }
}
- if (m.y != 0.) {
- y = (alpha - m.x)/m.y;
- if (y >= 0. && y <= 1.)
- l = g_slist_prepend (l, gts_point_new (gts_point_class (), -0.5, 0.5 - y, 0.));
+ if (m->y != 0.) {
+ y = (alpha - m->x)/m->y;
+ if (y >= 0. && y <= 1.) {
+ p[n].x = q.x - h/2.; p[n].y = q.y + h*(0.5 - y); p[n++].z = 0.;
+ }
}
- if (m.x != 0.) {
- x = (alpha - m.y)/m.x;
- if (x >= 0. && x <= 1.)
- l = g_slist_prepend (l, gts_point_new (gts_point_class (), 0.5 - x, -0.5, 0.));
+ if (m->x != 0.) {
+ x = (alpha - m->y)/m->x;
+ if (x >= 0. && x <= 1.) {
+ p[n].x = q.x + h*(0.5 - x); p[n].y = q.y - h/2.; p[n++].z = 0.;
+ }
+ }
+ g_assert (n == 2);
+#else /* 3D */
+ gdouble max = fabs (m->x);
+ FttComponent c = FTT_X;
+ if (fabs (m->y) > max) {
+ max = fabs (m->y);
+ c = FTT_Y;
}
- g_assert (g_slist_length (l) == 2);
- return l;
+ if (fabs (m->z) > max)
+ c = FTT_Z;
+ q.x -= h/2.; q.y -= h/2.; q.z -= h/2.;
+ (&q.x)[c] -= h*(alpha - m->x - m->y - m->z)/(&m->x)[c];
+ gts_vector_normalize (&m->x);
+
+ FttDirection d[12];
+ n = gfs_cut_cube_vertices (cell, -1, &q, m, p, d, NULL, NULL);
+ g_assert (n <= 6);
+#endif /* 3D */
+ return n;
}
}
diff --git a/src/vof.h b/src/vof.h
index b359d97..68f7143 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -66,8 +66,10 @@ gboolean gfs_vof_plane (FttCell * cell,
GfsVariable * v,
FttVector * m,
gdouble * alpha);
-GSList * gfs_vof_facet (FttCell * cell,
- GfsVariable * v);
+guint gfs_vof_facet (FttCell * cell,
+ GfsVariable * v,
+ FttVector * p,
+ FttVector * m);
gdouble gfs_vof_interpolate (FttCell * cell,
FttVector * p,
guint level,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list