[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203

[Stephane Popinet]

The following commit has been merged in the upstream branch:
commit 50df3f5ed636c2447f40689e3a0f456b7d72bb91
Author: Stephane Popinet <popinet at users.sf.net>
Date:   Mon Mar 19 16:07:37 2007 +1100

    gfs_plane_center() works with any normal
    
    Not just with normals in the positive quadrant.
    
    darcs-hash:20070319050737-d4795-174b3fa96c31179c55cb3c9cbfa48186f7564433.gz

diff --git a/src/vof.c b/src/vof.c
index c02f076..0d8cd27 100644
--- a/src/vof.c
+++ b/src/vof.c
@@ -128,21 +128,31 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
 
   g_return_if_fail (m != NULL);
   g_return_if_fail (p != NULL);
-  g_return_if_fail (m->x >= 0. && m->y >= 0.);
 
+  n = *m;
+  if (n.x < 0.) {
+    alpha -= n.x;
+    n.x = - n.x;
+  }
+  if (n.y < 0.) {
+    alpha -= n.y;
+    n.y = - n.y;
+  }
+
+  p->z = 0.;
   if (alpha <= 0.) {
     p->x = p->y = 0.;
     return;
   }
 
-  if (alpha >= m->x + m->y) {
+  if (alpha >= n.x + n.y) {
     p->x = p->y = 0.5;
     return;
   }
 
   g_return_if_fail (a > 0. && a < 1.);
 
-  n = *m; n.x += 1e-4; n.y += 1e-4;
+  n.x += 1e-4; n.y += 1e-4;
 
   p->x = p->y = alpha*alpha*alpha;
 
@@ -160,6 +170,11 @@ void gfs_line_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
   
   p->x /= 6.*n.x*n.x*n.y*a;
   p->y /= 6.*n.x*n.y*n.y*a;
+
+  if (m->x < 0.)
+    p->x = 1. - p->x;
+  if (m->y < 0.)
+    p->y = 1. - p->y;
 }
 
 #if (!FTT_2D)
@@ -315,21 +330,34 @@ void gfs_plane_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
 
   g_return_if_fail (m != NULL);
   g_return_if_fail (p != NULL);
-  g_return_if_fail (m->x >= 0. && m->y >= 0. && m->z >= 0.);
+
+  n = *m;
+  if (n.x < 0.) {
+    alpha -= n.x;
+    n.x = - n.x;
+  }
+  if (n.y < 0.) {
+    alpha -= n.y;
+    n.y = - n.y;
+  }  
+  if (n.z < 0.) {
+    alpha -= n.z;
+    n.z = - n.z;
+  }  
 
   if (alpha <= 0.) {
     p->x = p->y = p->z = 0.;
     return;
   }
 
-  if (alpha >= m->x + m->y + m->z) {
+  if (alpha >= n.x + n.y + n.z) {
     p->x = p->y = p->z = 0.5;
     return;
   }
 
   g_return_if_fail (a > 0. && a < 1.);
 
-  n = *m; n.x += 1e-4; n.y += 1e-4; n.z += 1e-4;
+  n.x += 1e-4; n.y += 1e-4; n.z += 1e-4;
 
   amax = n.x + n.y + n.z;
   p->x = p->y = p->z = alpha*alpha*alpha*alpha;
@@ -375,6 +403,10 @@ void gfs_plane_center (FttVector * m, gdouble alpha, gdouble a, FttVector * p)
 
   b = 24.*n.x*n.y*n.z*a;
   p->x /= b*n.x; p->y /= b*n.y; p->z /= b*n.z;
+
+  if (m->x < 0.) p->x = 1. - p->x;
+  if (m->y < 0.) p->y = 1. - p->y;
+  if (m->z < 0.) p->z = 1. - p->z;
 }
 #endif /* 3D */
 
diff --git a/src/vof.h b/src/vof.h
index bf06d15..68ed68e 100644
--- a/src/vof.h
+++ b/src/vof.h
@@ -40,6 +40,7 @@ gdouble gfs_line_alpha             (FttVector * m,
 #if FTT_2D
 #  define gfs_plane_volume         gfs_line_area
 #  define gfs_plane_alpha          gfs_line_alpha
+#  define gfs_plane_center         gfs_line_center
 #else /* 3D */
 gdouble gfs_plane_volume           (FttVector * m, 
 				    gdouble alpha);

-- 
Gerris Flow Solver