[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:53:47 UTC 2009
The following commit has been merged in the upstream branch:
commit ac932930bb14cbd10c248c4e21d6b2d81c7a1215
Author: Stephane Popinet <popinet at users.sf.net>
Date: Tue Nov 7 16:54:44 2006 +1100
More robust algorithm for computation of local interface height
darcs-hash:20061107055444-d4795-2d93a9161727ab151cac1adf0f416a0157cd2492.gz
diff --git a/src/vof.c b/src/vof.c
index 0ef76a2..f5adc61 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -897,42 +897,41 @@ static guint local_height (FttVector * p,
gdouble * H)
{
gdouble h = ftt_level_size (level);
- gdouble f1 = fraction (p, level, v);
- guint n = !GFS_IS_FULL (f1);
-
- *H = f1;
- FttVector i = *p;
- (&i.x)[c] -= h;
- f1 = fraction (&i, level, v);
- while (!GFS_IS_FULL (f1) && n < NMAX) {
- // fprintf (stderr, " f1: %g (%g,%g)\n", f1, i.x, i.y);
- *H += f1; n++;
- (&i.x)[c] -= h;
- f1 = fraction (&i, level, v);
- }
- if (f1 > 0.5)
- *H += ((&i.x)[c] - (&origin->x)[c])/h + 0.5;
-
- i = *p;
- (&i.x)[c] += h;
- gdouble f2 = fraction (&i, level, v);
- while (!GFS_IS_FULL (f2) && n < NMAX) {
- // fprintf (stderr, " f2: %g (%g,%g)\n", f2, i.x, i.y);
- *H += f2; n++;
- (&i.x)[c] += h;
- f2 = fraction (&i, level, v);
- }
- if (f2 > 0.5) {
- if (f1 > 0.5)
- return 0;
- *H = 0.5 + *H - ((&i.x)[c] - (&origin->x)[c])/h;
- }
- else if (f1 < 0.5)
+ gdouble right = fraction (p, level, v), left = right;
+ FttVector pright = *p, pleft = pright;
+ guint n = 1;
+
+ *H = right;
+ while (fabs (right - left) < 1. && n < NMAX)
+ if (GFS_IS_FULL (left)) {
+ (&pright.x)[c] += h;
+ right = fraction (&pright, level, v); n++;
+ *H += right;
+ if (fabs (right - left) < 1.) {
+ (&pleft.x)[c] -= h;
+ left = fraction (&pleft, level, v); n++;
+ *H += left;
+ }
+ }
+ else {
+ (&pleft.x)[c] -= h;
+ left = fraction (&pleft, level, v); n++;
+ *H += left;
+ if (fabs (right - left) < 1.) {
+ (&pright.x)[c] += h;
+ right = fraction (&pright, level, v); n++;
+ *H += right;
+ }
+ }
+ if (fabs (right - left) < 1.)
return 0;
-
- // fprintf (stderr, " %d %g ", n, *H);
-
- return n >= NMAX ? 0 : n;
+ if (left > 0.5)
+ *H += ((&pleft.x)[c] - (&origin->x)[c])/h - 0.5;
+ else {
+ g_assert (right > 0.5);
+ *H -= ((&pright.x)[c] - (&origin->x)[c])/h + 0.5;
+ }
+ return n;
}
static FttComponent orientation (FttVector * m)
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list