[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:38 UTC 2009
The following commit has been merged in the upstream branch:
commit 17218e8605aadc180433321c489278ddd6ce0e0a
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Feb 21 03:18:27 2006 +1100
Simplification of gfs_plane_volume()
darcs-hash:20060220161827-d4795-1009c06583db1ad422a838c26a9cc5d4da362039.gz
diff --git a/src/solid.c b/src/solid.c
index 6ba524a..1f98a0f 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -627,18 +627,25 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s, In
sym[c] = FALSE;
n += (&m.x)[c];
}
- g_assert (n > 0.);
- m.x /= n; m.y /= n; m.z /= n;
- alpha = m.x*ca.x + m.y*ca.y + m.z*ca.z;
- solid->a = gfs_plane_volume (&m, alpha, 1.);
- gfs_plane_center (&m, alpha, solid->a, &solid->cm);
- for (c = 0; c < FTT_DIMENSION; c++)
- (&solid->cm.x)[c] = (&o.x)[c] +
- (sym[c] ? 1. - (&solid->cm.x)[c] : (&solid->cm.x)[c])*(&h.x)[c];
- }
- else { /* this is a "thin" cell */
- p->thin++;
- deal_with_thin_cell (cell, p);
+ if (n > 0.) {
+ m.x /= n; m.y /= n; m.z /= n;
+ alpha = m.x*ca.x + m.y*ca.y + m.z*ca.z;
+ solid->a = gfs_plane_volume (&m, alpha, 1.);
+ gfs_plane_center (&m, alpha, solid->a, &solid->cm);
+ for (c = 0; c < FTT_DIMENSION; c++)
+ (&solid->cm.x)[c] = (&o.x)[c] +
+ (sym[c] ? 1. - (&solid->cm.x)[c] : (&solid->cm.x)[c])*(&h.x)[c];
+ }
+ else { /* degenerate intersections */
+ solid->a = 0.;
+ for (i = 0; i < FTT_NEIGHBORS; i++)
+ solid->a += solid->s[i];
+ solid->a /= FTT_NEIGHBORS;
+ if (solid->a == 0. || solid->a == 1.) {
+ g_free (solid);
+ GFS_STATE (cell)->solid = NULL;
+ }
+ }
}
if (solid->a == 0.)
GFS_VARIABLE (cell, p->status->i) = 1.;
diff --git a/src/vof.c b/src/vof.c
index c8c6113..d8d53b5 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -27,12 +27,11 @@
* gfs_line_area:
* @m: normal to the line.
* @alpha: line constant.
- * @c1: width of the cell.
*
- * Returns: the area of the fraction of a rectangular cell (c1,1)
- * lying under the line (@m, at alpha).
+ * Returns: the area of the fraction of a cell lying under the line
+ * (@m, at alpha).
*/
-gdouble gfs_line_area (FttVector * m, gdouble alpha, gdouble c1)
+gdouble gfs_line_area (FttVector * m, gdouble alpha)
{
FttVector n;
gdouble a, v;
@@ -43,14 +42,14 @@ gdouble gfs_line_area (FttVector * m, gdouble alpha, gdouble c1)
if (alpha <= 0.)
return 0.;
- if (alpha >= m->x*c1 + m->y || c1 == 0.)
- return c1;
+ if (alpha >= m->x + m->y)
+ return 1.;
n = *m; n.x += 1e-4; n.y += 1e-4;
v = alpha*alpha;
- a = alpha - n.x*c1;
+ a = alpha - n.x;
if (a > 0.)
v -= a*a;
@@ -117,12 +116,10 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
* gfs_plane_volume:
* @m: normal to the plane.
* @alpha: plane constant.
- * @c1: width of the cell.
*
- * Returns: the volume of a parallelepipedic cell (c1,1,1) lying under
- * the plane (@m, at alpha).
+ * Returns: the volume of a cell lying under the plane (@m, at alpha).
*/
-gdouble gfs_plane_volume (FttVector * m, gdouble alpha, gdouble c1)
+gdouble gfs_plane_volume (FttVector * m, gdouble alpha)
{
FttVector n;
gdouble a, amax, v;
@@ -135,29 +132,23 @@ gdouble gfs_plane_volume (FttVector * m, gdouble alpha, gdouble c1)
if (alpha <= 0.)
return 0.;
- if (alpha >= m->x*c1 + m->y + m->z || c1 == 0.)
- return c1;
+ if (alpha >= m->x + m->y + m->z)
+ return 1.;
n = *m; n.x += 1e-4; n.y += 1e-4; n.z += 1e-4;
- amax = n.x*c1 + n.y + n.z;
+ amax = n.x + n.y + n.z;
md = &n.x;
v = alpha*alpha*alpha;
- a = alpha - n.x*c1;
- if (a > 0.)
- v -= a*a*a;
- for (j = 1; j < 3; j++) {
+ for (j = 0; j < 3; j++) {
a = alpha - md[j];
if (a > 0.)
v -= a*a*a;
}
amax = alpha - amax;
- a = amax + n.x*c1;
- if (a > 0.)
- v += a*a*a;
- for (j = 1; j < 3; j++) {
+ for (j = 0; j < 3; j++) {
a = amax + md[j];
if (a > 0.)
v += a*a*a;
@@ -470,15 +461,19 @@ static void gfs_cell_vof_advected_face_values (FttCell * cell,
m.x /= 1. + u_right - u_left;
alpha += m.x*u_left;
- if (u_left < 0.)
- GFS_STATE (cell)->f[left].v =
- - gfs_plane_volume (&m, alpha - m.x*u_left, - u_left)/u_left;
+ if (u_left < 0.) {
+ m1 = m;
+ m1.x *= - u_left;
+ GFS_STATE (cell)->f[left].v = gfs_plane_volume (&m1, alpha + m1.x);
+ }
else
GFS_STATE (cell)->f[left].v = f;
- if (u_right > 0.)
- GFS_STATE (cell)->f[right].v =
- gfs_plane_volume (&m, alpha - m.x, u_right)/u_right;
+ if (u_right > 0.) {
+ m1 = m;
+ m1.x *= u_right;
+ GFS_STATE (cell)->f[right].v = gfs_plane_volume (&m1, alpha - m.x);
+ }
else
GFS_STATE (cell)->f[right].v = f;
}
@@ -667,7 +662,7 @@ void gfs_vof_coarse_fine (FttCell * parent, GfsVariable * v)
(&p.x)[c] = (&m.x)[c] < 0. ? (&p.x)[c] + 0.25 : 0.25 - (&p.x)[c];
alpha1 -= (&m1.x)[c]*(&p.x)[c];
}
- GFS_VARIABLE (child.c[i], v->i) = gfs_plane_volume (&m1, 2.*alpha1, 1.);
+ GFS_VARIABLE (child.c[i], v->i) = gfs_plane_volume (&m1, 2.*alpha1);
}
}
}
diff --git a/src/vof.h b/src/vof.h
index abe952e..9c8a317 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -29,8 +29,7 @@ extern "C" {
#define GFS_IS_FULL(f) ((f) < 1e-6 || (f) > 1. - 1.e-6)
gdouble gfs_line_area (FttVector * m,
- gdouble alpha,
- gdouble c1);
+ gdouble alpha);
void gfs_line_center (FttVector * m,
gdouble alpha,
gdouble a,
@@ -42,8 +41,7 @@ gdouble gfs_line_alpha (FttVector * m,
# define gfs_plane_alpha gfs_line_alpha
#else /* 3D */
gdouble gfs_plane_volume (FttVector * m,
- gdouble alpha,
- gdouble c1);
+ gdouble alpha);
gdouble gfs_plane_alpha (FttVector * m,
gdouble c);
void gfs_plane_center (FttVector * m,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list