[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:57 UTC 2009
The following commit has been merged in the upstream branch:
commit 7685c8d76555444fb06a64abe232ff32d037206e
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Feb 12 12:46:56 2007 +1100
Mixed Youngs-Centered VOF normal calculation
darcs-hash:20070212014656-d4795-d705d4950a4c896fcd5f6ad2df615c617084d618.gz
diff --git a/AUTHORS b/AUTHORS
index 0e3f87f..678bf29 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -5,4 +5,5 @@ St
Contributors
------------
Marcelo E. Magallon: Debian packages.
-Ruben Scardovelli: author of the Fortran version of gfs_plane_alpha()
+Ruben Scardovelli: - author of the Fortran version of gfs_plane_alpha()
+ - Mixed Youngs-Centered VOF normal calculation
diff --git a/src/Makefile.am b/src/Makefile.am
index 5ded918..224a48b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -86,6 +86,7 @@ SRC = \
utils.c \
ocean.c \
levelset.c \
+ myc.h \
$(GFS_HDS)
libgfs3D_la_LDFLAGS = $(NO_UNDEFINED)\
diff --git a/src/myc.h b/src/myc.h
new file mode 100644
index 0000000..1cf061c
--- /dev/null
+++ b/src/myc.h
@@ -0,0 +1,145 @@
+#define NOT_ZERO 1.e-30
+
+/*-----------------------------------------------------*
+ *MYC - Mixed Youngs and Central Scheme *
+ *-----------------------------------------------------*/
+/*
+
+Known problems: the index [1][1][1], i.e. the central cell
+in the block, never occurs: neither in the central scheme
+nor in Youngs' method. Therefore an isolated droplet will have
+a normal with all components to zero. I took care of the
+division-by-zero issue, but not of this one.
+
+Ruben
+
+*/
+static void mycs(double c[3][3][3],double mxyz[3])
+{
+ double m1,m2,m[4][3],t0,t1,t2;
+ int cn;
+
+ /* write the plane as: sgn(mx) X = my Y + mz Z + alpha
+ m00 X = m01 Y + m02 Z + alpha */
+ m1 = c[0][1][0] + c[0][1][2] + c[0][0][1] + c[0][2][1] +
+ c[0][1][1];
+ m2 = c[2][1][0] + c[2][1][2] + c[2][0][1] + c[2][2][1] +
+ c[2][1][1];
+ t0 = fabs(m1-m2) + NOT_ZERO;
+ m[0][0] = (m1-m2)/t0;
+
+ m1 = c[0][0][1]+ c[2][0][1]+ c[1][0][1];
+ m2 = c[0][2][1]+ c[2][2][1]+ c[1][2][1];
+ m[0][1] = 0.5*(m1-m2);
+
+ m1 = c[0][1][0]+ c[2][1][0]+ c[1][1][0];
+ m2 = c[0][1][2]+ c[2][1][2]+ c[1][1][2];
+ m[0][2] = 0.5*(m1-m2);
+
+ /* write the plane as: sgn(my) Y = mx X + mz Z + alpha
+ m11 Y = m10 X + m12 Z + alpha */
+ m1 = c[0][0][1] + c[0][2][1] + c[0][1][1];
+ m2 = c[2][0][1] + c[2][2][1] + c[2][1][1];
+ m[1][0] = 0.5*(m1-m2);
+
+ m1 = c[1][0][0] + c[1][0][2] + c[2][0][1] + c[0][0][1] +
+ c[1][0][1];
+ m2 = c[1][2][0] + c[1][2][2] + c[2][2][1] + c[0][2][1] +
+ c[1][2][1];
+ t0 = fabs(m1-m2) + NOT_ZERO;
+ m[1][1] = (m1-m2)/t0;
+
+ m1 = c[1][0][0]+ c[1][1][0]+ c[1][2][0];
+ m2 = c[1][0][2]+ c[1][1][2]+ c[1][2][2];
+ m[1][2] = 0.5*(m1-m2);
+
+ /* write the plane as: sgn(mz) Z = mx X + my Y + alpha
+ m22 Z = m20 X + m21 Y + alpha */
+
+ m1 = c[0][1][0]+ c[0][1][2]+ c[0][1][1];
+ m2 = c[2][1][0]+ c[2][1][2]+ c[2][1][1];
+ m[2][0] = 0.5*(m1-m2);
+
+ m1 = c[1][0][0]+ c[1][0][2]+ c[1][0][1];
+ m2 = c[1][2][0]+ c[1][2][2]+ c[1][2][1];
+ m[2][1] = 0.5*(m1-m2);
+
+ m1 = c[0][1][0] + c[2][1][0] + c[1][0][0] + c[1][2][0] +
+ c[1][1][0];
+ m2 = c[0][1][2] + c[2][1][2] + c[1][0][2] + c[1][2][2] +
+ c[1][1][2];
+ t0 = fabs(m1-m2) + NOT_ZERO;
+ m[2][2] = (m1-m2)/t0;
+
+ /* normalize each set (mx,my,mz): |mx|+|my|+|mz| = 1 */
+ t0 = fabs(m[0][0]) + fabs(m[0][1]) + fabs(m[0][2]) + NOT_ZERO;
+ m[0][0] /= t0;
+ m[0][1] /= t0;
+ m[0][2] /= t0;
+
+ t0 = fabs(m[1][0]) + fabs(m[1][1]) + fabs(m[1][2]) + NOT_ZERO;
+ m[1][0] /= t0;
+ m[1][1] /= t0;
+ m[1][2] /= t0;
+
+ t0 = fabs(m[2][0]) + fabs(m[2][1]) + fabs(m[2][2]) + NOT_ZERO;
+ m[2][0] /= t0;
+ m[2][1] /= t0;
+ m[2][2] /= t0;
+
+ /* choose among the three central scheme */
+ t0 = fabs(m[0][0]);
+ t1 = fabs(m[1][1]);
+ t2 = fabs(m[2][2]);
+
+ cn = 0;
+ if (t1 > t0 && t1 > t2)
+ cn = 1;
+ if (t2 > t0 && t2 > t1)
+ cn = 2;
+
+ /* Youngs-CIAM scheme */
+ m1 = c[0][0][0] + c[0][2][0] + c[0][0][2] + c[0][2][2] +
+ 2.*(c[0][0][1] + c[0][2][1] + c[0][1][0] + c[0][1][2]) +
+ 4.*c[0][1][1];
+ m2 = c[2][0][0] + c[2][2][0] + c[2][0][2] + c[2][2][2] +
+ 2.*(c[2][0][1] + c[2][2][1] + c[2][1][0] + c[2][1][2]) +
+ 4.*c[2][1][1];
+ m[3][0] = m1-m2;
+
+ m1 = c[0][0][0] + c[0][0][2] + c[2][0][0] + c[2][0][2] +
+ 2.*( c[0][0][1] + c[2][0][1] + c[1][0][0] + c[1][0][2]) +
+ 4.*c[1][0][1];
+ m2 = c[0][2][0] + c[0][2][2] + c[2][2][0] + c[2][2][2] +
+ 2.*(c[0][2][1] + c[2][2][1] + c[1][2][0] + c[1][2][2]) +
+ 4.*c[1][2][1];
+ m[3][1] = m1-m2;
+
+ m1 = c[0][0][0] + c[0][2][0] + c[2][0][0] + c[2][2][0] +
+ 2.*(c[0][1][0] + c[2][1][0] + c[1][0][0] + c[1][2][0]) +
+ 4.*c[1][1][0];
+ m2 = c[0][0][2] + c[0][2][2] + c[2][0][2] + c[2][2][2] +
+ 2.*(c[0][1][2] + c[2][1][2] + c[1][0][2] + c[1][2][2]) +
+ 4.*c[1][1][2];
+ m[3][2] = m1-m2;
+
+ /* normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1 */
+ t0 = fabs(m[3][0]) + fabs(m[3][1]) + fabs(m[3][2]) + NOT_ZERO;
+ m[3][0] /= t0;
+ m[3][1] /= t0;
+ m[3][2] /= t0;
+
+ /* choose between the previous choice and Youngs-CIAM */
+ t0 = MAX(fabs(m[3][0]),fabs(m[3][1]));
+ t0 = MAX(t0,fabs(m[3][2]));
+
+ if (fabs(m[cn][cn]) > t0)
+ cn = 3;
+
+ /* components of the normal vector */
+ mxyz[0] = m[cn][0];
+ mxyz[1] = m[cn][1];
+ mxyz[2] = m[cn][2];
+
+ return;
+}
diff --git a/src/vof.c b/src/vof.c
index efc0ae6..29a21a2 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -562,6 +562,23 @@ static void youngs_normal (FttCell * cell, GfsVariable * v, FttVector * n)
#endif /* 3D */
}
+#include "myc.h"
+
+static void myc_normal (FttCell * cell, GfsVariable * v, FttVector * n)
+{
+ gdouble f[3][3][3];
+
+ stencil (cell, v, f);
+#if FTT_2D
+ /* fixme: this is still Youngs */
+ n->x = (f[0][2][1] + 2.*f[0][1][1] + f[0][0][1] - 2.*f[2][1][1] - f[2][2][1] - f[2][0][1])/8.;
+ n->y = (f[2][0][1] + 2.*f[1][0][1] + f[0][0][1] - 2.*f[1][2][1] - f[2][2][1] - f[0][2][1])/8.;
+ n->z = 0.;
+#else /* 3D */
+ mycs (f, &n->x);
+#endif /* 3D */
+}
+
static void vof_plane (FttCell * cell, GfsVariable * v)
{
if (FTT_CELL_IS_LEAF (cell)) {
@@ -580,7 +597,7 @@ static void vof_plane (FttCell * cell, GfsVariable * v)
FttVector m;
gdouble n = 0.;
- youngs_normal (cell, v, &m);
+ myc_normal (cell, v, &m);
for (c = 0; c < FTT_DIMENSION; c++)
n += fabs ((&m.x)[c]);
if (n > 0.)
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list