[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sourceforge.net
Fri May 15 02:51:15 UTC 2009
The following commit has been merged in the upstream branch:
commit d033a637f114323d1279b99d65c94e07a99aa6d4
Author: Stephane Popinet <popinet at users.sourceforge.net>
Date: Thu Oct 14 13:31:51 2004 +1000
gfs_interpolate uses proper bi(tri)linear interpolation (gerris--mainline--0.5--patch-6)
gerris--mainline--0.5--patch-6
Keywords:
Corner values for the cell are computed using gfs_cell_corner_value
and then used to do the bi(tri)linear interpolation. This has been
tested succinctly and provides true continuous (C1)
interpolation. This was not the case before.
An important note is that the value at the center of the cell obtained
by (bi)trilinear interpolation is NOT equal to the variable value at
the center (it is equal to the mean of the corner values).
darcs-hash:20041014033151-aabb8-fcd5fa9c0ff1349fa32c0a10b163c69dcffcb0c7.gz
diff --git a/src/domain.c b/src/domain.c
index a17d1ce..f551ec7 100644
--- a/src/domain.c
+++ b/src/domain.c
@@ -1684,12 +1684,13 @@ FttCell * gfs_domain_locate (GfsDomain * domain,
* If @p is not contained within @domain, the coordinates are unchanged.
*/
void gfs_domain_advect_point (GfsDomain * domain,
- GtsPoint * p,
- gdouble dt)
+ GtsPoint * p,
+ gdouble dt)
{
FttCell * cell;
FttVector p0, p1;
FttComponent c;
+ GfsVariable * v1, * v;
g_return_if_fail (domain != NULL);
g_return_if_fail (p != NULL);
@@ -1700,13 +1701,15 @@ void gfs_domain_advect_point (GfsDomain * domain,
cell = gfs_domain_locate (domain, p0, -1);
if (cell == NULL)
return;
- for (c = 0; c < FTT_DIMENSION; c++)
- (&p1.x)[c] += dt*gfs_interpolate (cell, p0, GFS_VELOCITY_INDEX (c))/2.;
+ v1 = v = gfs_variable_from_name (domain->variables, "U");
+ for (c = 0; c < FTT_DIMENSION; c++, v = v->next)
+ (&p1.x)[c] += dt*gfs_interpolate (cell, p0, v)/2.;
cell = gfs_domain_locate (domain, p1, -1);
if (cell == NULL)
return;
- for (c = 0; c < FTT_DIMENSION; c++)
- (&p->x)[c] += dt*gfs_interpolate (cell, p1, GFS_VELOCITY_INDEX (c));
+ v = v1;
+ for (c = 0; c < FTT_DIMENSION; c++, v = v->next)
+ (&p->x)[c] += dt*gfs_interpolate (cell, p1, v);
}
static void count (FttCell * cell, guint * n)
@@ -2036,8 +2039,6 @@ void gfs_cell_read_binary (FttCell * cell, GtsFile * fp, GfsDomain * domain)
gfs_cell_init (cell, domain);
s = cell->data;
if (s0 >= 0.) {
- guint i;
-
s->solid = g_malloc0 (sizeof (GfsSolidVector));
s->solid->s[0] = s0;
diff --git a/src/fluid.c b/src/fluid.c
index b164be3..178de9b 100644
--- a/src/fluid.c
+++ b/src/fluid.c
@@ -19,6 +19,8 @@
#include <math.h>
+#include <stdlib.h>
+
#include "fluid.h"
#include "domain.h"
#include "solid.h"
@@ -731,6 +733,7 @@ static gboolean inverse (gdouble mi[N_CELLS - 1][N_CELLS - 1])
#endif /* 3D */
return TRUE;
}
+
#if (!FTT_2D)
static void draw_cell (FttCell * cell, gdouble r, gdouble g, gdouble b,
const gchar * name)
@@ -1796,64 +1799,84 @@ void gfs_cell_traverse_mixed (FttCell * root,
* gfs_interpolate:
* @cell: a #FttCell containing location @p.
* @p: the location at which to interpolate.
- * @v: a #GfsVariable index.
+ * @v: a #GfsVariable.
*
* Interpolates the @v variable of @cell, at location @p. Linear
* interpolation is used and the boundaries of the domain are treated
* as planes of symmetry for all variables.
*
- * Note that to work correctly this function needs a domain where the
- * values of variable @v are defined on all levels.
- *
* Returns: the interpolated value of variable @v at location @p.
*/
gdouble gfs_interpolate (FttCell * cell,
FttVector p,
- guint v)
+ GfsVariable * v)
{
- gboolean leveled;
FttVector o;
- gdouble f[FTT_DIMENSION], size, v0, val;
- FttComponent c;
+ FttDirection e[FTT_DIMENSION];
+ gdouble size;
g_return_val_if_fail (cell != NULL, 0.);
- do {
- guint level;
-
- ftt_cell_pos (cell, &o);
- level = ftt_cell_level (cell);
- leveled = TRUE;
-
- for (c = 0; c < FTT_DIMENSION && leveled; c++) {
- FttCell * neighbor;
-
- if (((gdouble *) &p)[c] >= ((gdouble *) &o)[c])
- neighbor = ftt_cell_neighbor (cell, 2*c);
- else
- neighbor = ftt_cell_neighbor (cell, 2*c + 1);
-
- if (neighbor) {
- if (ftt_cell_level (neighbor) != level)
- leveled = FALSE;
- else
- f[c] = GFS_VARIABLE (neighbor, v);
- }
- else /* symmetry conditions */
- f[c] = GFS_VARIABLE (cell, v);
- }
- if (!leveled) {
- cell = ftt_cell_parent (cell);
- g_assert (cell);
- }
- } while (!leveled);
-
- size = ftt_cell_size (cell);
- val = v0 = GFS_VARIABLE (cell, v);
- for (c = 0; c < FTT_DIMENSION; c++)
- val += fabs (((gdouble *) &p)[c] - ((gdouble *) &o)[c])*(f[c] - v0)/size;
-
- return val;
+ ftt_cell_pos (cell, &o);
+ size = ftt_cell_size (cell)/2.;
+ p.x = (p.x - o.x)/size;
+ p.y = (p.y - o.y)/size;
+#if FTT_2D
+ {
+ gdouble f[4], a, b, c, d;
+
+ e[0] = FTT_LEFT; e[1] = FTT_BOTTOM;
+ f[0] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_BOTTOM;
+ f[1] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_TOP;
+ f[2] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_LEFT; e[1] = FTT_TOP;
+ f[3] = gfs_cell_corner_value (cell, e, v, -1);
+ a = f[1] + f[2] - f[0] - f[3];
+ b = f[2] + f[3] - f[0] - f[1];
+ c = f[0] - f[1] + f[2] - f[3];
+ d = f[0] + f[1] + f[2] + f[3];
+ return (a*p.x + b*p.y + c*p.x*p.y + d)/4.;
+ }
+#else /* 3D */
+ {
+ gdouble f[8], c[8];
+
+ p.z = (p.z - o.z)/size;
+ e[0] = FTT_LEFT; e[1] = FTT_BOTTOM; e[2] = FTT_FRONT;
+ f[0] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_BOTTOM;
+ f[1] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_TOP;
+ f[2] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_LEFT; e[1] = FTT_TOP;
+ f[3] = gfs_cell_corner_value (cell, e, v, -1);
+
+ e[0] = FTT_LEFT; e[1] = FTT_BOTTOM; e[2] = FTT_BACK;
+ f[4] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_BOTTOM;
+ f[5] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_RIGHT; e[1] = FTT_TOP;
+ f[6] = gfs_cell_corner_value (cell, e, v, -1);
+ e[0] = FTT_LEFT; e[1] = FTT_TOP;
+ f[7] = gfs_cell_corner_value (cell, e, v, -1);
+
+ c[0] = - f[0] + f[1] + f[2] - f[3] - f[4] + f[5] + f[6] - f[7];
+ c[1] = - f[0] - f[1] + f[2] + f[3] - f[4] - f[5] + f[6] + f[7];
+ c[2] = f[0] + f[1] + f[2] + f[3] - f[4] - f[5] - f[6] - f[7];
+ c[3] = f[0] - f[1] + f[2] - f[3] + f[4] - f[5] + f[6] - f[7];
+ c[4] = - f[0] + f[1] + f[2] - f[3] + f[4] - f[5] - f[6] + f[7];
+ c[5] = - f[0] - f[1] + f[2] + f[3] + f[4] + f[5] - f[6] - f[7];
+ c[6] = f[0] - f[1] + f[2] - f[3] - f[4] + f[5] - f[6] + f[7];
+ c[7] = f[0] + f[1] + f[2] + f[3] + f[4] + f[5] + f[6] + f[7];
+
+ return (c[0]*p.x + c[1]*p.y + c[2]*p.z +
+ c[3]*p.x*p.y + c[4]*p.x*p.z + c[5]*p.y*p.z +
+ c[6]*p.x*p.y*p.z +
+ c[7])/8.;
+ }
+#endif /* 3D */
}
/* GfsVariable: Object */
diff --git a/src/fluid.h b/src/fluid.h
index 010de42..5859a2d 100644
--- a/src/fluid.h
+++ b/src/fluid.h
@@ -279,7 +279,7 @@ void gfs_cell_traverse_mixed (FttCell * root,
gpointer data);
gdouble gfs_interpolate (FttCell * cell,
FttVector p,
- guint v);
+ GfsVariable * v);
GfsVariable * gfs_variable_list_copy (GfsVariable * v,
GtsObject * parent);
void gfs_variable_list_destroy (GfsVariable * v);
diff --git a/src/graphic.c b/src/graphic.c
index 3520514..6b05275 100644
--- a/src/graphic.c
+++ b/src/graphic.c
@@ -531,7 +531,7 @@ static void iso_func (gdouble ** a, GtsCartesianGrid g, guint k,
if (cell == NULL)
a[i][j] = 0.;
else
- a[i][j] = gfs_interpolate (cell, p, v->i);
+ a[i][j] = gfs_interpolate (cell, p, v);
}
}
@@ -1290,7 +1290,7 @@ static GtsColor variable_color (GtsObject * o)
cell = gfs_domain_locate (domain, pos, -1);
if (cell) {
- val = gfs_interpolate (cell, pos, v->i);
+ val = gfs_interpolate (cell, pos, v);
c = colormap_color (colormap, (val - *min)/(*max - *min));
}
else
@@ -1827,6 +1827,7 @@ static GSList * grow_curve (GfsDomain * domain,
guint nstep = 0, nmax = 10000;
GtsPointClass * path_class = gfs_vertex_class ();
Colormap * colormap = NULL;
+ GfsVariable * U;
if (min < max)
colormap = colormap_jet ();
@@ -1848,11 +1849,13 @@ static GSList * grow_curve (GfsDomain * domain,
}
#endif /* 3D */
+ U = gfs_variable_from_name (domain->variables, "U");
p1 = p2 = p;
while ((cell = gfs_domain_locate (domain, p, -1)) != NULL && nmax--) {
gdouble h = delta*ftt_cell_size (cell);
FttVector u;
FttComponent c;
+ GfsVariable * v;
gdouble nu = 0.;
cost += curve_cost (p1, p2, p);
@@ -1861,7 +1864,7 @@ static GSList * grow_curve (GfsDomain * domain,
if (oldp == NULL || cost > maxcost) {
oldp = gts_point_new (path_class, p.x, p.y, p.z);
if (var)
- GFS_VERTEX (oldp)->v = gfs_interpolate (cell, p, var->i);
+ GFS_VERTEX (oldp)->v = gfs_interpolate (cell, p, var);
if (colormap)
GTS_COLORED_VERTEX (oldp)->c =
colormap_color (colormap, (GFS_VERTEX (oldp)->v - min)/(max - min));
@@ -1872,9 +1875,8 @@ static GSList * grow_curve (GfsDomain * domain,
nstep = 0;
}
- for (c = 0; c < FTT_DIMENSION; c++) {
- ((gdouble *) &u)[c] =
- direction*gfs_interpolate (cell, p, GFS_VELOCITY_INDEX (c));
+ for (c = 0, v = U; c < FTT_DIMENSION; c++, v = v->next) {
+ ((gdouble *) &u)[c] = direction*gfs_interpolate (cell, p, v);
nu += ((gdouble *) &u)[c]*((gdouble *) &u)[c];
}
if (nu > 0. && nstep++ < nmax) {
@@ -1890,8 +1892,8 @@ static GSList * grow_curve (GfsDomain * domain,
GtsVector dx;
dx[0] = p1.x - p.x; dx[1] = p1.y - p.y; dx[2] = p1.z - p.z;
- for (c = 0; c < FTT_DIMENSION; c++)
- rot[c] = gfs_interpolate (cell, p1, GFS_GRADIENT_INDEX (c));
+ for (c = 0, v = gfs_gx; c < FTT_DIMENSION; c++, v = v->next)
+ rot[c] = gfs_interpolate (cell, p1, v);
theta += gts_vector_scalar (rot, dx)/nu;
}
#endif /* 3D */
@@ -1904,7 +1906,7 @@ static GSList * grow_curve (GfsDomain * domain,
if (cell) {
oldp = gts_point_new (path_class, p2.x, p2.y, p2.z);
if (var)
- GFS_VERTEX (oldp)->v = gfs_interpolate (cell, p2, var->i);
+ GFS_VERTEX (oldp)->v = gfs_interpolate (cell, p2, var);
if (twist)
GFS_TWISTED_VERTEX (oldp)->theta = theta;
path = g_slist_prepend (path, oldp);
diff --git a/src/output.c b/src/output.c
index d076cfe..9ca5773 100644
--- a/src/output.c
+++ b/src/output.c
@@ -1102,7 +1102,7 @@ static gboolean gfs_output_location_event (GfsEvent * event,
location->p.x, location->p.y, location->p.z);
while (v) {
if (v->name)
- fprintf (fp, " %g", gfs_interpolate (cell, location->p, v->i));
+ fprintf (fp, " %g", gfs_interpolate (cell, location->p, v));
v = v->next;
}
fputc ('\n', fp);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list