[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:47 UTC 2009
The following commit has been merged in the upstream branch:
commit 3a6fadae90048c2c5f7f6485a471100d861f51ac
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Oct 16 12:46:17 2008 +1100
Consistent z-scaling of terrain
darcs-hash:20081016014617-d4795-346d037c2a603d6d2a8e7213208ccb3d080ca13f.gz
diff --git a/modules/map.mod b/modules/map.mod
index a159bb7..858ff82 100644
--- a/modules/map.mod
+++ b/modules/map.mod
@@ -32,7 +32,7 @@ struct _GfsMapProjection {
gdouble cosa, sina;
/*< public >*/
- gdouble lon, lat, angle, zscale;
+ gdouble lon, lat, angle;
};
#define GFS_MAP_PROJECTION(obj) GTS_OBJECT_CAST (obj,\
@@ -55,14 +55,12 @@ static void gfs_map_projection_read (GtsObject ** o, GtsFile * fp)
{GTS_DOUBLE, "lon", TRUE},
{GTS_DOUBLE, "lat", TRUE},
{GTS_DOUBLE, "angle", TRUE},
- {GTS_DOUBLE, "zscale", TRUE},
{GTS_NONE}
};
GfsMapProjection * map = GFS_MAP_PROJECTION (*o);
var[0].data = &map->lon;
var[1].data = &map->lat;
var[2].data = &map->angle;
- var[3].data = &map->zscale;
gts_file_assign_variables (fp, var);
if (fp->type == GTS_ERROR)
@@ -70,8 +68,6 @@ static void gfs_map_projection_read (GtsObject ** o, GtsFile * fp)
map->cosa = cos (map->angle*DEG_TO_RAD);
map->sina = sin (map->angle*DEG_TO_RAD);
- if (!var[3].set)
- map->zscale = gfs_object_simulation (map)->physical_params.L;
char * parms[] = {
"proj=lcc", /* Lambert Conformal Conic */
@@ -94,8 +90,8 @@ static void gfs_map_projection_write (GtsObject * o, FILE * fp)
{
(* GTS_OBJECT_CLASS (gfs_map_projection_class ())->parent_class->write) (o, fp);
GfsMapProjection * map = GFS_MAP_PROJECTION (o);
- fprintf (fp, " { lon = %.8g lat = %.8g angle = %g zscale = %g }",
- map->lon, map->lat, map->angle, map->zscale);
+ fprintf (fp, " { lon = %.8g lat = %.8g angle = %g }",
+ map->lon, map->lat, map->angle);
}
static void gfs_map_projection_destroy (GtsObject * object)
@@ -115,7 +111,7 @@ static void projection_transform (GfsMap * map, const FttVector * src, FttVector
odata = pj_fwd (idata, m->pj);
dest->x = odata.u*m->cosa - odata.v*m->sina;
dest->y = odata.v*m->cosa + odata.u*m->sina;
- dest->z = src->z/m->zscale*gfs_object_simulation (map)->physical_params.L;
+ dest->z = src->z;
}
static void projection_inverse (GfsMap * map, const FttVector * src, FttVector * dest)
@@ -128,7 +124,7 @@ static void projection_inverse (GfsMap * map, const FttVector * src, FttVector *
odata = pj_inv (idata, GFS_MAP_PROJECTION (map)->pj);
dest->x = odata.u*RAD_TO_DEG;
dest->y = odata.v*RAD_TO_DEG;
- dest->z = src->z*m->zscale/gfs_object_simulation (map)->physical_params.L;
+ dest->z = src->z;
}
static void gfs_map_projection_class_init (GfsMapClass * klass)
@@ -145,7 +141,6 @@ static void gfs_map_projection_init (GfsMapProjection * object)
/* Wellington */
object->lon = 174.777222;
object->lat = -41.288889;
- object->zscale = 1.;
object->angle = 0.; object->cosa = 1.; object->sina = 0.;
object->pj = NULL;
}
diff --git a/modules/terrain.mod b/modules/terrain.mod
index 1bce143..2652c87 100644
--- a/modules/terrain.mod
+++ b/modules/terrain.mod
@@ -43,7 +43,7 @@ struct _GfsRefineTerrain {
#if !FTT_2D
GfsVariable * min, * max;
- gdouble front;
+ gdouble front, scale;
#endif
RSurface ** rs;
@@ -696,25 +696,31 @@ static void reset_terrain (FttCell * cell, GfsRefineTerrain * t)
# define traverse_boundary(domain,order,flags,depth,func,data) \
gfs_domain_cell_traverse_boundary(domain,FTT_FRONT,order,flags,depth,func,data)
-static void terrain_min_max (FttCell * cell, GfsVariable * h[NM], gdouble minmax[2])
+static void terrain_min_max (gdouble H[NM], gdouble minmax[2], gdouble scale)
{
gdouble dx, dy;
- gdouble H0 = GFS_VALUE (cell, h[0]), H1 = GFS_VALUE (cell, h[1]);
- gdouble H2 = GFS_VALUE (cell, h[2]), H3 = GFS_VALUE (cell, h[3]);
minmax[0] = G_MAXDOUBLE; minmax[1] = - G_MAXDOUBLE;
for (dx = -1.; dx <= 1.; dx += 2.)
for (dy = -1.; dy <= 1.; dy += 2.) {
- gdouble v = H0 + dx*H1 + dy*H2 + dx*dy*H3;
+ gdouble v = H[0] + dx*H[1] + dy*H[2] + dx*dy*H[3];
if (v < minmax[0]) minmax[0] = v;
if (v > minmax[1]) minmax[1] = v;
}
+ minmax[0] *= scale;
+ minmax[1] *= scale;
}
static void min_max (FttCell * cell, GfsRefineTerrain * t)
{
gdouble minmax[2] = { G_MAXDOUBLE, - G_MAXDOUBLE };
if (FTT_CELL_IS_LEAF (cell)) {
- terrain_min_max (cell, t->h, minmax);
+ gdouble h[4];
+ h[0] = GFS_VALUE (cell, t->h[0]);
+ h[1] = GFS_VALUE (cell, t->h[1]);
+ h[2] = GFS_VALUE (cell, t->h[2]);
+ h[3] = GFS_VALUE (cell, t->h[3]);
+ terrain_min_max (h, minmax, t->scale);
+
FttVector p;
ftt_cell_pos (cell, &p);
if (p.z > t->front)
@@ -883,15 +889,10 @@ static void terrain_coarse_fine (FttCell * parent, GfsVariable * v)
hc[3] = h[3]/4.;
#if !FTT_2D
ftt_cell_pos (child.c[n], &p);
- gdouble zmin = p.z - size, zmax = p.z + size;
- gdouble dx, dy;
- gdouble minmax[2] = { G_MAXDOUBLE, - G_MAXDOUBLE };
- for (dx = -1.; dx <= 1.; dx += 2.)
- for (dy = -1.; dy <= 1.; dy += 2.) {
- gdouble v = hc[0] + dx*hc[1] + dy*hc[2] + dx*dy*hc[3];
- if (v < minmax[0]) minmax[0] = v;
- if (v > minmax[1]) minmax[1] = v;
- }
+ gdouble zmin = p.z - size, zmax = p.z + size, minmax[2];
+ p.z = 1.;
+ gfs_simulation_map (GFS_SIMULATION (v->domain), &p);
+ terrain_min_max (hc, minmax, p.z);
if (minmax[0] > zmax || minmax[1] < zmin)
GFS_VALUE (child.c[n], v) = G_MAXDOUBLE;
else
@@ -938,6 +939,9 @@ static void terrain_refine (GfsRefine * refine, GfsSimulation * sim)
t->min = gfs_temporary_variable (domain);
t->max = gfs_temporary_variable (domain);
t->front = - G_MAXDOUBLE;
+ FttVector p = {0.,0.,1.};
+ gfs_simulation_map (sim, &p);
+ t->scale = p.z;
gfs_domain_cell_traverse_boundary (domain, FTT_FRONT, FTT_POST_ORDER, FTT_TRAVERSE_ALL, -1,
(FttCellTraverseFunc) min_max, t);
gts_container_foreach (GTS_CONTAINER (domain), (GtsFunc) refine_box, t);
@@ -960,13 +964,15 @@ static void refine_terrain_destroy (GtsObject * object)
GfsRefineTerrain * t = GFS_REFINE_TERRAIN (object);
GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (object));
- gchar * dname = g_strconcat (t->name, "min", NULL);
- gfs_domain_remove_derived_variable (domain, dname);
- g_free (dname);
-
- dname = g_strconcat (t->name, "max", NULL);
- gfs_domain_remove_derived_variable (domain, dname);
- g_free (dname);
+ if (t->name) {
+ gchar * dname = g_strconcat (t->name, "min", NULL);
+ gfs_domain_remove_derived_variable (domain, dname);
+ g_free (dname);
+
+ dname = g_strconcat (t->name, "max", NULL);
+ gfs_domain_remove_derived_variable (domain, dname);
+ g_free (dname);
+ }
g_free (t->name);
g_free (t->basename);
@@ -1222,6 +1228,7 @@ struct _GfsSurfaceTerrain {
/*< private >*/
GfsGenericSurface parent;
GfsVariable * h[NM];
+ gdouble scale;
/*< public >*/
gchar * name;
@@ -1282,6 +1289,16 @@ static GfsGenericSurface * cell_is_cut (FttCell * cell, GfsGenericSurface * s1,
return GFS_VALUE (cell, GFS_SURFACE_TERRAIN (s1)->h[0]) != G_MAXDOUBLE ? s1 : NULL;
}
+static gdouble zscale (GfsSurfaceTerrain * t)
+{
+ if (t->scale == 0.) {
+ FttVector p = {0.,0.,1.};
+ gfs_simulation_map (gfs_object_simulation (t), &p);
+ t->scale = p.z;
+ }
+ return t->scale;
+}
+
static guint surface_segment_intersection (GfsGenericSurface * s1,
FttCell * cell,
GfsSegment * I)
@@ -1293,8 +1310,9 @@ static guint surface_segment_intersection (GfsGenericSurface * s1,
FttVector pE, pD;
pE.x = I->E->x; pE.y = I->E->y;
pD.x = I->D->x; pD.y = I->D->y;
- gdouble vE = I->E->z - cell_value (cell, GFS_SURFACE_TERRAIN (s1)->h, pE);
- gdouble vD = I->D->z - cell_value (cell, GFS_SURFACE_TERRAIN (s1)->h, pD);
+ GfsSurfaceTerrain * t = GFS_SURFACE_TERRAIN (s1);
+ gdouble vE = I->E->z - cell_value (cell, t->h, pE)*zscale (t);
+ gdouble vD = I->D->z - cell_value (cell, t->h, pD)*zscale (t);
if ((vE > 0. && vD <= 0.) || (vE <= 0. && vD > 0.)) {
I->n = 1;
@@ -1341,7 +1359,7 @@ static void surface_segment_normal (GfsGenericSurface * s1,
p.y = (p.y - q.y)/size;
n[0] = - (GFS_VALUE (cell, h[1]) + GFS_VALUE (cell, h[3])*p.y)/size;
n[1] = - (GFS_VALUE (cell, h[2]) + GFS_VALUE (cell, h[3])*p.x)/size;
- n[2] = 1.;
+ n[2] = 1./zscale (GFS_SURFACE_TERRAIN (s1));
}
static void gfs_surface_terrain_class_init (GfsGenericSurfaceClass * klass)
diff --git a/src/simulation.c b/src/simulation.c
index e61c192..5571616 100644
--- a/src/simulation.c
+++ b/src/simulation.c
@@ -1489,7 +1489,7 @@ void gfs_simulation_map (GfsSimulation * sim, FttVector * p)
i = i->next;
}
FttComponent c;
- for (c = 0; c < FTT_DIMENSION; c++)
+ for (c = 0; c < 3; c++)
(&p->x)[c] *= (&GFS_DOMAIN (sim)->lambda.x)[c]/sim->physical_params.L;
}
@@ -1507,7 +1507,7 @@ void gfs_simulation_map_inverse (GfsSimulation * sim, FttVector * p)
g_return_if_fail (p != NULL);
FttComponent c;
- for (c = 0; c < FTT_DIMENSION; c++)
+ for (c = 0; c < 3; c++)
(&p->x)[c] *= sim->physical_params.L/(&GFS_DOMAIN (sim)->lambda.x)[c];
GSList * i = sim->maps->items;
while (i) {
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list