[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sourceforge.net
Fri May 15 02:51:20 UTC 2009
The following commit has been merged in the upstream branch:
commit 7ea1168f312cc7722aa8486881f85dcb765db116
Author: Stephane Popinet <popinet at users.sourceforge.net>
Date: Tue Nov 2 10:32:17 2004 +1100
New calculation of 3D center of mass of solid fractions (gerris--mainline--0.7--patch-19)
gerris--mainline--0.7--patch-19
Keywords:
Uses the new 3D VOF gfs_plane_center() function.
darcs-hash:20041101233217-aabb8-057abdd4f66f84f172fc3dd7ee0a1a2cb599ccbb.gz
diff --git a/src/solid.c b/src/solid.c
index 4327fc7..4ace14f 100644
--- a/src/solid.c
+++ b/src/solid.c
@@ -379,11 +379,16 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
}
/* now compute cell fraction, center of area, center of mass */
- ca.x /= n1; ca.y /= n1; ca.z /= n1;
+
+ /* fixme: the calculation of ca below is not strictly equivalent to
+ the true center of area of the piece of surface contained in the
+ cell if there are more than 4 intersection points (n1 > 4) */
+ ca.x /= n1; ca.y /= n1; ca.z /= n1;
solid->ca = ca;
{
FttVector m;
gdouble alpha, n = 0.;
+ gboolean sym[FTT_DIMENSION];
FttComponent c;
for (c = 0; c < FTT_DIMENSION; c++) {
@@ -392,13 +397,19 @@ static void set_solid_fractions_from_surface (FttCell * cell, GtsSurface * s)
if ((&m.x)[c] < 0.) {
(&m.x)[c] = - (&m.x)[c];
(&ca.x)[c] = 1. - (&ca.x)[c];
+ sym[c] = TRUE;
}
+ else
+ sym[c] = FALSE;
n += (&m.x)[c];
}
if (n > 0.) {
- m.x /= n; m.y /= n; m.z /= n;
+ m.x = m.x/n + 1e-6; m.y = m.y/n + 1e-6; m.z = m.z/n + 1e-6;
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;
}
else { /* degenerate intersections */
solid->a = 0.;
diff --git a/src/vof.c b/src/vof.c
index 57a6f60..91a0649 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -62,19 +62,19 @@ gdouble gfs_line_area (FttVector * m, gdouble alpha, gdouble c1)
* gfs_line_center:
* @m: normal to the line.
* @alpha: line constant.
- * @c: area of cell fraction.
+ * @a: area of cell fraction.
* @p: a #FttVector.
*
* Fills @p with the position of the center of mass of the fraction of
* a square cell lying under the line (@m, at alpha).
*/
-void gfs_line_center (FttVector * m, gdouble alpha, gdouble c, FttVector * p)
+void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
{
- gdouble a;
+ gdouble b;
g_return_if_fail (m != NULL);
g_return_if_fail (p != NULL);
- g_return_if_fail (c > 0. && c < 1.);
+ g_return_if_fail (a > 0. && a < 1.);
if (alpha <= 0.) {
p->x = p->y = 0.;
@@ -90,20 +90,20 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble c, FttVector * p)
p->x = p->y = alpha*alpha*alpha;
- a = alpha - m->x;
- if (a > 0.) {
- p->x -= a*a*(alpha + 2.*m->x);
- p->y -= a*a*a;
+ b = alpha - m->x;
+ if (b > 0.) {
+ p->x -= b*b*(alpha + 2.*m->x);
+ p->y -= b*b*b;
}
- a = alpha - m->y;
- if (a > 0.) {
- p->y -= a*a*(alpha + 2.*m->y);
- p->x -= a*a*a;
+ b = alpha - m->y;
+ if (b > 0.) {
+ p->y -= b*b*(alpha + 2.*m->y);
+ p->x -= b*b*b;
}
- p->x /= 6.*m->x*m->x*m->y*c;
- p->y /= 6.*m->x*m->y*m->y*c;
+ p->x /= 6.*m->x*m->x*m->y*a;
+ p->y /= 6.*m->x*m->y*m->y*a;
}
#if (!FTT_2D)
@@ -157,6 +157,85 @@ gdouble gfs_plane_volume (FttVector * m, gdouble alpha, gdouble c1)
return v/(6.*m->x*m->y*m->z);
}
+
+/**
+ * gfs_plane_center:
+ * @m: normal to the plane.
+ * @alpha: plane constant.
+ * @a: volume of cell fraction.
+ * @p: a #FttVector.
+ *
+ * Fills @p with the position of the center of mass of the fraction of
+ * a cubic cell lying under the plane (@m, at alpha).
+ */
+void gfs_plane_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
+{
+ gdouble b, amax;
+
+ g_return_if_fail (m != NULL);
+ g_return_if_fail (a > 0. && a < 1.);
+ g_return_if_fail (p != NULL);
+
+ if (alpha <= 0.) {
+ p->x = p->y = p->z = 0.;
+ return;
+ }
+
+ amax = m->x + m->y + m->z;
+ if (alpha >= amax) {
+ p->x = p->y = p->z = 0.5;
+ return;
+ }
+
+ g_assert (m->x >= 1e-9 && m->y >= 1e-9 && m->z >= 1e-9);
+
+ p->x = p->y = p->z = alpha*alpha*alpha;
+ p->x *= alpha/m->x;
+ p->y *= alpha/m->y;
+ p->z *= alpha/m->z;
+
+ b = alpha - m->x;
+ if (b > 0.) {
+ p->x -= b*b*b*(3. + alpha/m->x);
+ p->y -= b*b*b*b/m->y;
+ p->z -= b*b*b*b/m->z;
+ }
+ b = alpha - m->y;
+ if (b > 0.) {
+ p->y -= b*b*b*(3. + alpha/m->y);
+ p->x -= b*b*b*b/m->x;
+ p->z -= b*b*b*b/m->z;
+ }
+ b = alpha - m->z;
+ if (b > 0.) {
+ p->z -= b*b*b*(3. + alpha/m->z);
+ p->x -= b*b*b*b/m->x;
+ p->y -= b*b*b*b/m->y;
+ }
+
+ amax = alpha - amax;
+ b = amax + m->x;
+ if (b > 0.) {
+ p->y += b*b*b*(3. + (alpha - m->z)/m->y);
+ p->z += b*b*b*(3. + (alpha - m->y)/m->z);
+ p->x += b*b*b*b/m->x;
+ }
+ b = amax + m->y;
+ if (b > 0.) {
+ p->x += b*b*b*(3. + (alpha - m->z)/m->x);
+ p->z += b*b*b*(3. + (alpha - m->x)/m->z);
+ p->y += b*b*b*b/m->y;
+ }
+ b = amax + m->z;
+ if (b > 0.) {
+ p->x += b*b*b*(3. + (alpha - m->y)/m->x);
+ p->y += b*b*b*(3. + (alpha - m->x)/m->y);
+ p->z += b*b*b*b/m->z;
+ }
+
+ b = 24.*m->x*m->y*m->z*a;
+ p->x /= b; p->y /= b; p->z /= b;
+}
#endif /* 3D */
static gdouble line_area_derivative_ratio (FttVector * m,
diff --git a/src/vof.h b/src/vof.h
index 82b6495..40ee567 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -31,7 +31,7 @@ gdouble gfs_line_area (FttVector * m,
gdouble c1);
void gfs_line_center (FttVector * m,
gdouble alpha,
- gdouble c,
+ gdouble a,
FttVector * p);
gdouble gfs_line_alpha (FttVector * m,
gdouble c);
@@ -44,6 +44,10 @@ gdouble gfs_plane_volume (FttVector * m,
gdouble c1);
gdouble gfs_plane_alpha (FttVector * m,
gdouble c);
+void gfs_plane_center (FttVector * m,
+ gdouble alpha,
+ gdouble a,
+ FttVector * p);
#endif /* 3D */
void gfs_cell_vof_advection (FttCell * cell,
FttComponent c,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list