[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:15 UTC 2009
The following commit has been merged in the upstream branch:
commit d1e9e6aa1f16e1a55685ef848498052257496218
Author: Stephane Popinet <popinet at users.sf.net>
Date: Wed Apr 23 07:31:04 2008 +1000
Terrain module uses R*-tree
darcs-hash:20080422213104-d4795-863e43aed6016bc0c6cbde658a58ef470b5d950c.gz
diff --git a/modules/Makefile.am b/modules/Makefile.am
index 597e99e..df9bef8 100644
--- a/modules/Makefile.am
+++ b/modules/Makefile.am
@@ -9,15 +9,19 @@ if HAS_LIBPROJ
MAP = libmap2D.la libmap3D.la libmap2D3.la
endif
-TERRAIN = libterrain2D.la libterrain3D.la libterrain2D3.la
-
pkglib_LTLIBRARIES = \
$(MAP) \
- $(TERRAIN) \
+ libterrain2D.la \
+ libterrain3D.la \
+ libterrain2D3.la \
libtide2D.la \
libtide3D.la \
libtide2D3.la
+bin_PROGRAMS = \
+ xyz2rsurface \
+ rsurfacequery
+
EXTRA_DIST = \
map.mod \
tide.mod \
@@ -41,14 +45,16 @@ libmap2D3_la_SOURCES = map.c
libmap2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
libmap2D3_la_LIBADD = -lproj
-libterrain3D_la_SOURCES = terrain.c
-libterrain3D_la_LIBADD = -lpq
-libterrain2D_la_SOURCES = terrain.c
+libterrain3D_la_SOURCES = terrain.c rsurface.c
+libterrain3D_la_LIBADD = -LRStarTree -lcSmRST
+
libterrain2D_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D=1
-libterrain2D_la_LIBADD = -lpq
-libterrain2D3_la_SOURCES = terrain.c
+libterrain2D_la_SOURCES = terrain.c rsurface.c
+libterrain2D_la_LIBADD = -LRStarTree -lcSmRST
+
+libterrain2D3_la_SOURCES = terrain.c rsurface.c
libterrain2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
-libterrain2D3_la_LIBADD = -lpq
+libterrain2D3_la_LIBADD = -LRStarTree -lcSmRST
libtide3D_la_SOURCES = tide.c
libtide3D_la_LIBADD = -L/home/popinet/local/lib -lfes -lnetcdf -lgsl -lgslcblas -lm -lpthread
@@ -59,6 +65,12 @@ libtide2D3_la_SOURCES = tide.c
libtide2D3_la_CFLAGS = $(AM_CFLAGS) -DFTT_2D3=1
libtide2D3_la_LIBADD = -L/home/popinet/local/lib -lfes -lnetcdf -lgsl -lgslcblas -lm -lpthread
+xyz2rsurface_SOURCES = xyz2rsurface.c rsurface.c
+xyz2rsurface_LDADD = -LRStarTree -lcSmRST
+
+rsurfacequery_SOURCES = rsurfacequery.c rsurface.c
+rsurfacequery_LDADD = -LRStarTree -lcSmRST
+
if HAVE_MODULES
%.c : %.mod
@echo "/* $@" > $@
diff --git a/modules/rsurface.c b/modules/rsurface.c
new file mode 100644
index 0000000..82bf601
--- /dev/null
+++ b/modules/rsurface.c
@@ -0,0 +1,113 @@
+#include <stdio.h>
+#include <string.h>
+#include <malloc.h>
+#include "RStarTree/RStarTree.h"
+#include "rsurface.h"
+
+/* RSurface */
+
+struct _RSurface {
+ RSTREE t;
+ char * name;
+};
+
+RSurface * r_surface_new (const char * fname, int size, FILE * fp)
+{
+ return NULL;
+}
+
+RSurface * r_surface_open (const char * fname, const char * mode, int size)
+{
+ RSurface * rt = malloc (sizeof (RSurface));
+ if (!strcmp (mode, "w")) {
+ RemoveRST (fname);
+ if (!CreateRST (fname, size, FALSE)) {
+ free (rt);
+ return NULL;
+ }
+ }
+ rt->t = NULL;
+ if (!OpenRST (&rt->t, fname)) {
+ free (rt);
+ return NULL;
+ }
+ rt->name = malloc (sizeof (char)*(strlen (fname) + 1));
+ strcpy (rt->name, fname);
+ return rt;
+}
+
+void r_surface_close (RSurface * rt)
+{
+ CloseRST (&rt->t);
+ free (rt->name);
+ free (rt);
+}
+
+int r_surface_insert (RSurface * rt, double p[3], int id)
+{
+ typrect rect;
+ typinfo info;
+ boolean inserted;
+ rect[0].l = p[0]; rect[0].h = p[0];
+ rect[1].l = p[1]; rect[1].h = p[1];
+ info.height = p[2];
+ return (InsertRecord (rt->t, rect, &info, &inserted) && inserted);
+}
+
+void r_surface_query (RSurface * rt,
+ double a[2], double b[2], double c[2], double d[2],
+ RSurfaceQuery q, void * user_data)
+{
+}
+
+static boolean Intersects (RSTREE R,
+ typrect RSTrect,
+ typrect queryrect,
+ typrect unused)
+{
+ int maxdim= NumbOfDim -1;
+ boolean inter;
+ int d;
+
+ d= -1;
+ do {
+ d++;
+ inter= RSTrect[d].l <= queryrect[d].h &&
+ RSTrect[d].h >= queryrect[d].l;
+ } while (inter && d != maxdim);
+ return inter;
+}
+
+static void ManageQuery (RSTREE R,
+ typrect rectangle,
+ refinfo infoptr,
+ void ** data,
+ boolean *modify,
+ boolean *finish)
+{
+ RSurfaceQuery q = data[0];
+ double p[3];
+ p[0] = rectangle[0].l; p[1] = rectangle[1].l;
+ p[2] = infoptr->height;
+ (*q) (p, data[1]);
+ *modify = FALSE;
+ *finish = FALSE;
+}
+
+void r_surface_query_region (RSurface * rt,
+ double min[2], double max[2],
+ RSurfaceQuery q, void * user_data)
+{
+ typrect rect, unused;
+ rect[0].l = min[0]; rect[0].h = max[0];
+ rect[1].l = min[1]; rect[1].h = max[1];
+ void * data[2];
+ data[0] = q;
+ data[1] = user_data;
+ RegionQuery (rt->t, rect, unused, Intersects, Intersects, ManageQuery, data);
+}
+
+const char * r_surface_name (RSurface * rt)
+{
+ return rt->name;
+}
diff --git a/modules/rsurface.h b/modules/rsurface.h
new file mode 100644
index 0000000..0112937
--- /dev/null
+++ b/modules/rsurface.h
@@ -0,0 +1,15 @@
+typedef struct _RSurface RSurface;
+
+RSurface * r_surface_new (const char * fname, int size, FILE * fp);
+RSurface * r_surface_open (const char * fname, const char * mode, int size);
+void r_surface_close (RSurface * rt);
+int r_surface_insert (RSurface * rt, double p[3], int id);
+
+typedef int (* RSurfaceQuery) (double p[3], void * user_data);
+void r_surface_query (RSurface * rt,
+ double a[2], double b[2], double c[2], double d[2],
+ RSurfaceQuery q, void * user_data);
+void r_surface_query_region (RSurface * rt,
+ double min[2], double max[2],
+ RSurfaceQuery q, void * user_data);
+const char * r_surface_name (RSurface * rt);
diff --git a/modules/rsurfacequery.c b/modules/rsurfacequery.c
new file mode 100644
index 0000000..e11a23a
--- /dev/null
+++ b/modules/rsurfacequery.c
@@ -0,0 +1,31 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "rsurface.h"
+
+static int query (double p[3], void * data)
+{
+ printf ("%.8f %.8f %g\n", p[0], p[1], p[2]);
+ return 0;
+}
+
+int main (int argc, char** argv)
+{
+ if (argc != 2) {
+ fprintf (stderr, "Usage: %s basename\n", argv[0]);
+ return -1;
+ }
+
+ RSurface * rs = r_surface_open (argv[1], "r", 0);
+ double min[2], max[2];
+ int count = 0;
+ while (scanf ("%lf %lf %lf %lf", &min[0], &min[1], &max[0], &max[1]) == 4) {
+ r_surface_query_region (rs, min, max, query, NULL);
+ if (count > 0 && count % 1000 == 0)
+ fprintf (stderr, "\r%d", count);
+ count++;
+ }
+ fputc ('\n', stderr);
+ r_surface_close (rs);
+
+ return 0.;
+}
diff --git a/modules/terrain.mod b/modules/terrain.mod
index d0ee45a..d2b49be 100644
--- a/modules/terrain.mod
+++ b/modules/terrain.mod
@@ -17,54 +17,31 @@
* 02111-1307, USA.
*/
#include <stdlib.h>
-#include <postgresql/libpq-fe.h>
-#include "refine.h"
-
-#define DBMAX "1000"
-
-#ifdef G_LITTLE_ENDIAN
-# define DBYTE "NDR"
-#else
-# define DBYTE "XDR"
+#include <glob.h>
+#if GSL
+# include <gsl/gsl_linalg.h>
#endif
-#define DOUBLE_OID 701
-#define INT32_OID 23
-
-static gboolean inverse (gdouble m[3][3], gdouble mi[3][3])
-{
- gdouble det;
-
- det = (m[0][0]*(m[1][1]*m[2][2] - m[2][1]*m[1][2]) -
- m[0][1]*(m[1][0]*m[2][2] - m[2][0]*m[1][2]) +
- m[0][2]*(m[1][0]*m[2][1] - m[2][0]*m[1][1]));
- if (det == 0.)
- return FALSE;
-
- mi[0][0] = (m[1][1]*m[2][2] - m[1][2]*m[2][1])/det;
- mi[0][1] = (m[2][1]*m[0][2] - m[0][1]*m[2][2])/det;
- mi[0][2] = (m[0][1]*m[1][2] - m[1][1]*m[0][2])/det;
- mi[1][0] = (m[1][2]*m[2][0] - m[1][0]*m[2][2])/det;
- mi[1][1] = (m[0][0]*m[2][2] - m[2][0]*m[0][2])/det;
- mi[1][2] = (m[1][0]*m[0][2] - m[0][0]*m[1][2])/det;
- mi[2][0] = (m[1][0]*m[2][1] - m[2][0]*m[1][1])/det;
- mi[2][1] = (m[2][0]*m[0][1] - m[0][0]*m[2][1])/det;
- mi[2][2] = (m[0][0]*m[1][1] - m[0][1]*m[1][0])/det;
+#include "refine.h"
+#include "rsurface.h"
- return TRUE;
-}
+static gchar * default_dir = "/home/popinet/Projects/GIS/rsurfaces";
/* GfsRefineTerrain: Header */
typedef struct _GfsRefineTerrain GfsRefineTerrain;
+#define NM 4
+
struct _GfsRefineTerrain {
GfsRefine parent;
- PGconn * conn;
- int wkb_oid;
+ guint level;
+ gboolean refined;
+
+ gchar * name, * basename;
+ RSurface ** rs;
+ guint nrs;
- gchar * name;
- gchar * host, * dbname, * user, * passwd, * tables;
- GfsVariable * h, * hx, * hy, * he, * hn;
+ GfsVariable * h[NM], * he, * hn, * cond;
GfsFunction * criterion;
};
@@ -78,272 +55,520 @@ GfsRefineClass * gfs_refine_terrain_class (void);
/* GfsRefineTerrain: Object */
-static void update_terrain (FttCell * cell, GfsRefineTerrain * t)
+typedef struct {
+ FttVector c;
+ FttVector p[4];
+ gdouble min[2], max[2], h;
+ GfsRefineTerrain * t;
+ FttCell * cell;
+} Polygon;
+
+static void map_inverse (GfsSimulation * sim, FttVector * p, gdouble min[2], gdouble max[2])
{
- if (GFS_VALUE (cell, t->h) == G_MAXDOUBLE) {
- FttVector p, q[4];
- ftt_cell_pos (cell, &p);
- gdouble h = ftt_cell_size (cell)/2.;
- GfsSimulation * sim = gfs_object_simulation (t);
-
- /* build the cell boundaries in the non-projected coordinates system */
- q[0].x = p.x + h; q[0].y = p.y + h;
- gfs_simulation_map_inverse (sim, &q[0]);
- q[1].x = p.x - h; q[1].y = p.y + h;
- gfs_simulation_map_inverse (sim, &q[1]);
- q[2].x = p.x - h; q[2].y = p.y - h;
- gfs_simulation_map_inverse (sim, &q[2]);
- q[3].x = p.x + h; q[3].y = p.y - h;
- gfs_simulation_map_inverse (sim, &q[3]);
-
- /* get all the terrain data within this polygon */
- gchar * polygon =
- g_strdup_printf ("SRID=4326;POLYGON((%.8g %.8g,%.8g %.8g,%.8g %.8g,%.8g %.8g,%.8g %.8g))",
- q[0].x, q[0].y,
- q[1].x, q[1].y,
- q[2].x, q[2].y,
- q[3].x, q[3].y,
- q[0].x, q[0].y);
- // fprintf (stderr, "%s\n", polygon);
- PGresult * res = PQexecPrepared (t->conn, "get_points_in_polygon", 1,
- (const char *const*) &polygon, NULL, NULL, 1);
- g_free (polygon);
- if (!res || PQresultStatus (res) != PGRES_TUPLES_OK) {
- g_log (G_LOG_DOMAIN, G_LOG_LEVEL_CRITICAL,
- "PostgreSQL command failed\n%s", PQerrorMessage (t->conn));
- PQclear (res);
- return;
- }
+ gfs_simulation_map_inverse (sim, p);
+ if (p->x < min[0]) min[0] = p->x;
+ if (p->x > max[0]) max[0] = p->x;
+ if (p->y < min[1]) min[1] = p->y;
+ if (p->y > max[1]) max[1] = p->y;
+}
- /* least-squares fit */
- guint i, n = PQntuples (res);
- gdouble H[4] = {0.,0.,0.,0.}, m[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}, mi[3][3];
- m[0][0] = n;
- for (i = 0; i < n; i++) {
- gint32 height = g_ntohl (*((gint32 *) PQgetvalue (res, i, 0)));
- gchar * wkb = PQgetvalue (res, i, 1);
- FttVector r;
- memcpy (&r.x, &wkb[5], 8);
- memcpy (&r.y, &wkb[5+8], 8);
- gfs_simulation_map (sim, &r);
- r.x -= p.x; r.y -= p.y;
- m[0][1] += r.x; m[0][2] += r.y;
- m[1][1] += r.x*r.x; m[1][2] += r.x*r.y;
- m[2][2] += r.y*r.y;
- H[0] += height;
- H[1] += r.x*height;
- H[2] += r.y*height;
- H[3] += height*height;
- }
- PQclear (res);
-
- m[1][0] = m[0][1]; m[2][0] = m[0][2]; m[2][1] = m[1][2];
- if (n > 2 && inverse (m, mi)) {
- gdouble h0, hx, hy;
- h0 = GFS_VALUE (cell, t->h) = mi[0][0]*H[0] + mi[0][1]*H[1] + mi[0][2]*H[2];
- hx = GFS_VALUE (cell, t->hx) = mi[1][0]*H[0] + mi[1][1]*H[1] + mi[1][2]*H[2];
- hy = GFS_VALUE (cell, t->hy) = mi[2][0]*H[0] + mi[2][1]*H[1] + mi[2][2]*H[2];
- GFS_VALUE (cell, t->he) =
- sqrt (fabs (h0*(h0*m[0][0] + 2.*(hx*m[0][1] + hy*m[0][2] - H[0])) +
- hx*(hx*m[1][1] + 2.*(hy*m[1][2] - H[1])) +
- hy*(hy*m[2][2] - 2.*H[2]) +
- H[3]));
- }
- else { /* not enough points for fit, use parent info */
- FttCell * parent = ftt_cell_parent (cell);
- if (parent) {
- gdouble h0 = GFS_VALUE (parent, t->h);
- gdouble hx = GFS_VALUE (parent, t->hx);
- gdouble hy = GFS_VALUE (parent, t->hy);
- FttVector pp;
- ftt_cell_pos (parent, &pp);
- GFS_VALUE (cell, t->h) = h0 + (p.x - pp.x)*hx + (p.y - pp.y)*hy;
- GFS_VALUE (cell, t->hx) = hx;
- GFS_VALUE (cell, t->hy) = hy;
- GFS_VALUE (cell, t->he) = GFS_VALUE (parent, t->he);
- }
- else {
- g_warning ("cannot initialise terrain for cell at level %d", ftt_cell_level (cell));
- GFS_VALUE (cell, t->h) = G_MAXDOUBLE;
- GFS_VALUE (cell, t->hx) = G_MAXDOUBLE;
- GFS_VALUE (cell, t->hy) = G_MAXDOUBLE;
- GFS_VALUE (cell, t->he) = G_MAXDOUBLE;
- }
- }
- GFS_VALUE (cell, t->hn) = n;
- }
+static void polygon_init (Polygon * p, FttCell * cell, GfsRefineTerrain * t)
+{
+ GfsSimulation * sim = gfs_object_simulation (t);
+ FttVector q;
+ ftt_cell_pos (cell, &q);
+ p->cell = cell;
+ p->t = t;
+ p->c = q;
+ p->h = ftt_cell_size (cell)/2.;
+ p->min[0] = p->min[1] = G_MAXDOUBLE;
+ p->max[0] = p->max[1] = - G_MAXDOUBLE;
+ p->p[0].x = q.x + p->h; p->p[0].y = q.y + p->h;
+ map_inverse (sim, &p->p[0], p->min, p->max);
+ p->p[1].x = q.x - p->h; p->p[1].y = q.y + p->h;
+ map_inverse (sim, &p->p[1], p->min, p->max);
+ p->p[2].x = q.x - p->h; p->p[2].y = q.y - p->h;
+ map_inverse (sim, &p->p[2], p->min, p->max);
+ p->p[3].x = q.x + p->h; p->p[3].y = q.y - p->h;
+ map_inverse (sim, &p->p[3], p->min, p->max);
+}
+
+static gboolean right (const double a[2], const double b[2], const double c[2])
+{
+ return (b[0] - a[0])*(c[1] - a[1]) - (b[1] - a[1])*(c[0] - a[0]) < 0.;
}
-static gboolean refine_terrain (FttCell * cell, GfsRefine * refine)
+static gboolean polygon_contains (Polygon * p, gdouble q[2])
{
- update_terrain (cell, GFS_REFINE_TERRAIN (refine));
- if (ftt_cell_level (cell) >= gfs_function_value (refine->maxlevel, cell))
+ if (right (&p->p[0].x, &p->p[1].x, q))
+ return FALSE;
+ if (right (&p->p[1].x, &p->p[2].x, q))
+ return FALSE;
+ if (right (&p->p[2].x, &p->p[3].x, q))
+ return FALSE;
+ if (right (&p->p[3].x, &p->p[0].x, q))
return FALSE;
- return gfs_function_value (GFS_REFINE_TERRAIN (refine)->criterion, cell);
+ return TRUE;
}
-static void refine_box (GfsBox * box, GfsRefine * refine)
+typedef struct {
+ gdouble H[NM+1], m[NM][NM];
+ gdouble h[NM], he, cond, min, max;
+ Polygon * p;
+ gboolean relative;
+} RMS;
+
+static void rms_init (RMS * rms, Polygon * p, gboolean relative)
{
- ftt_cell_refine (box->root,
- (FttCellRefineFunc) refine_terrain, refine,
- (FttCellInitFunc) gfs_cell_fine_init, gfs_box_domain (box));
+ guint i, j;
+ for (i = 0; i < NM + 1; i++)
+ rms->H[i] = 0.;
+ for (i = 0; i < NM; i++)
+ for (j = 0; j < NM; j++)
+ rms->m[i][j] = 0.;
+ rms->p = p;
+ rms->relative = relative;
+ rms->min = G_MAXDOUBLE;
+ rms->max = - G_MAXDOUBLE;
}
-static void set_undefined (FttCell * cell, GfsVariable * h)
+static gdouble rms_minimum (RMS * rms)
{
- GFS_VALUE (cell, h) = G_MAXDOUBLE;
+ if (rms->m[0][0] == 0.)
+ return 0.;
+ return sqrt (fabs (rms->h[0]*(rms->h[0]*rms->m[0][0] +
+ 2.*(rms->h[1]*rms->m[0][1] +
+ rms->h[2]*rms->m[0][2] +
+ rms->h[3]*rms->m[0][3] - rms->H[0])) +
+ rms->h[1]*(rms->h[1]*rms->m[1][1] +
+ 2.*(rms->h[2]*rms->m[1][2] +
+ rms->h[3]*rms->m[1][3] - rms->H[1])) +
+ rms->h[2]*(rms->h[2]*rms->m[2][2] +
+ 2.*(rms->h[3]*rms->m[2][3] - rms->H[2])) +
+ rms->h[3]*(rms->h[3]*rms->m[3][3] +
+ - 2.*rms->H[3]) +
+ rms->H[4])/rms->m[0][0]);
}
-static void set_undefined_coarse_fine (FttCell * parent, GfsVariable * h)
+static void function_from_corners (gdouble h[4], gdouble H[4])
{
- FttCellChildren child;
- guint n;
+ h[0] = (H[0] + H[1] + H[2] + H[3])/4.;
+ h[1] = (H[0] - H[1] - H[2] + H[3])/4.;
+ h[2] = (H[0] + H[1] - H[2] - H[3])/4.;
+ h[3] = (H[0] - H[1] + H[2] - H[3])/4.;
+}
- ftt_cell_children (parent, &child);
- for (n = 0; n < FTT_CELLS; n++)
- if (child.c[n])
- GFS_VALUE (child.c[n], h) = G_MAXDOUBLE;
+static gdouble cell_value (FttCell * cell, GfsRefineTerrain * t, FttVector p)
+{
+ gdouble h = ftt_cell_size (cell)/2.;
+ FttVector q;
+ ftt_cell_pos (cell, &q);
+ p.x = (p.x - q.x)/h;
+ p.y = (p.y - q.y)/h;
+ return GFS_VALUE (cell, t->h[0]) +
+ GFS_VALUE (cell, t->h[1])*p.x +
+ GFS_VALUE (cell, t->h[2])*p.y +
+ GFS_VALUE (cell, t->h[3])*p.x*p.y;
}
-static void terrain_refine (GfsRefine * refine, GfsSimulation * sim)
+static void corners_from_parent (FttCell * cell, GfsRefineTerrain * t, gdouble H[4])
{
- GfsRefineTerrain * t = GFS_REFINE_TERRAIN (refine);
- gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) set_undefined, t->h);
- t->h->coarse_fine = set_undefined_coarse_fine;
- gts_container_foreach (GTS_CONTAINER (sim), (GtsFunc) refine_box, refine);
- t->h->coarse_fine = gfs_cell_coarse_fine;
- gfs_domain_cell_traverse (GFS_DOMAIN (sim), FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
- (FttCellTraverseFunc) update_terrain, refine);
+ gdouble size = ftt_cell_size (cell);
+ FttCell * parent = ftt_cell_parent (cell);
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ p.x += size/2.; p.y += size/2.;
+ H[0] = cell_value (parent, t, p);
+ p.x -= size;
+ H[1] = cell_value (parent, t, p);
+ p.y -= size;
+ H[2] = cell_value (parent, t, p);
+ p.x += size;
+ H[3] = cell_value (parent, t, p);
}
-static void refine_terrain_destroy (GtsObject * object)
+static gdouble clamp (gdouble x, gdouble min, gdouble max)
{
- GfsRefineTerrain * t = GFS_REFINE_TERRAIN (object);
- g_free (t->name);
- if (t->conn)
- PQfinish (t->conn);
- gts_object_destroy (GTS_OBJECT (t->criterion));
- g_free (t->host);
- g_free (t->dbname);
- g_free (t->user);
- g_free (t->passwd);
- g_free (t->tables);
- (* GTS_OBJECT_CLASS (gfs_refine_terrain_class ())->parent_class->destroy) (object);
+ if (x > max) return max;
+ if (x < min) return min;
+ return x;
}
-static PGresult * query (PGconn * conn, const gchar * s, int status, GtsFile * fp)
+static void variance_check (RMS * rms)
{
- PGresult * res = PQexec (conn, s);
- if (!res || PQresultStatus (res) != status) {
- gts_file_error (fp, "PostgreSQL command failed\n%s", PQerrorMessage (conn));
- PQclear (res);
- PQfinish (conn);
- return NULL;
- }
- if (status != PGRES_TUPLES_OK)
- PQclear (res);
- return res;
-}
-
-static gboolean postgres_connection (GfsRefineTerrain * t, GtsFile * fp)
-{
- gchar * q = g_strjoin (" ",
- "host =", t->host,
- "dbname =", t->dbname,
- "user =", t->user,
- "password =", t->passwd,
- NULL);
- PGconn * conn = PQconnectdb (q);
- g_free (q);
- if (PQstatus (conn) == CONNECTION_BAD) {
- gts_file_error (fp, "could not connect to database\n%s", PQerrorMessage (conn));
- PQfinish (conn);
- return FALSE;
- }
- PGresult * res = query (conn, "SELECT OID FROM pg_type WHERE typname = 'bytea'",
- PGRES_TUPLES_OK, fp);
- if (!res)
- return FALSE;
- if (PQntuples (res) != 1) {
- gts_file_error (fp, "query to find OID of geometry did not return 1 row");
- PQclear (res);
- PQfinish (conn);
- return FALSE;
- }
- t->wkb_oid = atoi (PQgetvalue (res, 0, 0));
- PQclear (res);
-
- q = g_strjoin (NULL,
- "BEGIN; DECLARE mycursor BINARY CURSOR FOR SELECT height, ASBINARY(geom,'"
- DBYTE "') FROM ",
- t->tables,
- " LIMIT 1",
- NULL);
- if (!query (conn, q, PGRES_COMMAND_OK, fp)) {
- g_free (q);
- return FALSE;
- }
- g_free (q);
- if (!(res = query (conn, "FETCH ALL IN mycursor", PGRES_TUPLES_OK, fp)))
- return FALSE;
- if (PQftype (res, 0) != INT32_OID) {
- gts_file_error (fp, "invalid type (%d) for field 'height'", PQftype (res, 0));
- PQclear (res);
- PQfinish (conn);
- return FALSE;
+ g_assert (rms->m[0][0] >= NM);
+ gdouble H[4];
+ if (rms->relative) {
+ gdouble H0[4];
+ corners_from_parent (rms->p->cell, rms->p->t, H0);
+ H[0] = clamp (rms->h[0] + rms->h[1] + rms->h[2] + rms->h[3],
+ rms->min - H0[0], rms->max - H0[0]);
+ H[1] = clamp (rms->h[0] - rms->h[1] + rms->h[2] - rms->h[3],
+ rms->min - H0[1], rms->max - H0[1]);
+ H[2] = clamp (rms->h[0] - rms->h[1] - rms->h[2] + rms->h[3],
+ rms->min - H0[2], rms->max - H0[2]);
+ H[3] = clamp (rms->h[0] + rms->h[1] - rms->h[2] - rms->h[3],
+ rms->min - H0[3], rms->max - H0[3]);
}
- if (PQftype (res, 1) != t->wkb_oid) {
- gts_file_error (fp, "invalid type (%d) for field 'geom'", PQftype (res, 1));
- PQclear (res);
- PQfinish (conn);
- return FALSE;
+ else {
+ H[0] = clamp (rms->h[0] + rms->h[1] + rms->h[2] + rms->h[3], rms->min, rms->max);
+ H[1] = clamp (rms->h[0] - rms->h[1] + rms->h[2] - rms->h[3], rms->min, rms->max);
+ H[2] = clamp (rms->h[0] - rms->h[1] - rms->h[2] + rms->h[3], rms->min, rms->max);
+ H[3] = clamp (rms->h[0] + rms->h[1] - rms->h[2] - rms->h[3], rms->min, rms->max);
}
- if (PQntuples (res) == 0) {
- gts_file_error (fp, "database seems to be empty");
- PQclear (res);
- PQfinish (conn);
- return FALSE;
+ function_from_corners (rms->h, H);
+}
+
+static void rms_update (RMS * rms)
+{
+ guint i;
+ if (rms->m[0][0] == 0.) {
+ for (i = 1; i < NM; i++)
+ rms->h[i] = 0.;
+ rms->h[0] = G_MAXDOUBLE;
+ rms->he = 0.;
+ rms->cond = G_MAXDOUBLE;
+ return;
}
- gchar * wkb = PQgetvalue (res, 0, 1);
- guint32 type;
- memcpy (&type, wkb + 1, 4);
- PQclear (res);
-#ifdef G_LITTLE_ENDIAN
- if (wkb[0] != 1) {
+ else if (rms->m[0][0] >= NM) {
+ guint j;
+ for (i = 1; i < NM; i++)
+ for (j = 0; j < i; j++)
+ rms->m[i][j] = rms->m[j][i];
+#if GSL
+ double m[NM*NM], v[NM*NM], s[NM];
+ gsl_matrix_view gv = gsl_matrix_view_array (v, NM, NM);
+ gsl_vector_view gs = gsl_vector_view_array (s, NM);
+ for (i = 0; i < NM; i++)
+ for (j = 0; j < NM; j++)
+ m[i+NM*j] = rms->m[i][j];
+ gsl_matrix_view gm = gsl_matrix_view_array (m, NM, NM);
+ gsl_linalg_SV_decomp_jacobi (&gm.matrix, &gv.matrix, &gs.vector);
+ rms->cond = s[NM - 1] > 0. ? s[0]/s[NM - 1] : G_MAXDOUBLE;
+ if (rms->cond < 10000.) {
+ gsl_vector_view gH = gsl_vector_view_array (rms->H, NM);
+ gsl_vector_view gh = gsl_vector_view_array (rms->h, NM);
+ gsl_linalg_SV_solve (&gm.matrix, &gv.matrix, &gs.vector, &gH.vector, &gh.vector);
+ variance_check (rms);
+ rms->he = rms_minimum (rms);
+ return;
+ }
#else
- if (wkb[0] != 0) {
+ gdouble ** m = gfs_matrix_new (NM, sizeof (gdouble));
+ for (i = 0; i < NM; i++)
+ for (j = 0; j < NM; j++)
+ m[i][j] = rms->m[i][j];
+ if (gfs_matrix_inverse (m, NM, 1e-5)) {
+ for (i = 0; i < NM; i++) {
+ rms->h[i] = 0.;
+ for (j = 0; j < NM; j++)
+ rms->h[i] += m[i][j]*rms->H[j];
+ }
+ gfs_matrix_free (m);
+ variance_check (rms);
+ rms->he = rms_minimum (rms);
+ return;
+ }
+ gfs_matrix_free (m);
#endif
- gts_file_error (fp, "endianness do not match");
- PQfinish (conn);
- return FALSE;
}
- if (type != 1) {
- gts_file_error (fp, "geometry type (%d) is not POINT", type);
- PQfinish (conn);
- return FALSE;
+ rms->h[0] = rms->H[0]/rms->m[0][0];
+ for (i = 1; i < NM; i++)
+ rms->h[i] = 0.;
+ rms->he = rms_minimum (rms);
+}
+
+#if DEBUG
+static gdouble rms_value (RMS * rms, FttVector * p)
+{
+ return rms->h[0] + rms->h[1]*p->x + rms->h[2]*p->y + rms->h[3]*p->x*p->y;
+}
+
+static void rms_write (RMS * rms, Polygon * p)
+{
+ FttVector q, r;
+ q.x = p->c.x + p->h; q.y = p->c.y + p->h;
+ r.x = 1.; r.y = 1.;
+ fprintf (stderr, "%g %g %g\n", q.x, q.y, rms_value (rms, &r)/20000.);
+ q.x = p->c.x + p->h; q.y = p->c.y - p->h;
+ r.x = 1.; r.y = - 1.;
+ fprintf (stderr, "%g %g %g\n", q.x, q.y, rms_value (rms, &r)/20000.);
+ q.x = p->c.x - p->h; q.y = p->c.y - p->h;
+ r.x = - 1.; r.y = - 1.;
+ fprintf (stderr, "%g %g %g\n", q.x, q.y, rms_value (rms, &r)/20000.);
+ q.x = p->c.x - p->h; q.y = p->c.y + p->h;
+ r.x = - 1.; r.y = 1.;
+ fprintf (stderr, "%g %g %g\n", q.x, q.y, rms_value (rms, &r)/20000.);
+}
+
+static int write_points (double p[3], RMS * rms)
+{
+ if (polygon_contains (rms->p, p))
+ fprintf (stderr, "aa %g %g %g\n", p[0], p[1], p[2]);
+}
+#endif
+
+#define RAW 0. /* fitted but not continuous (C0) */
+#define FAIR 1. /* fitted and C0 */
+#define REFINED 2. /* non-fitted refined cell */
+#define NEW_CHILD 3. /* non-fitted child extrapolated from its parent */
+
+static int rms_add (double p[3], RMS * rms)
+{
+ if (polygon_contains (rms->p, p)) {
+ gfs_simulation_map (gfs_object_simulation (rms->p->t), (FttVector *) p);
+ if (p[2] > rms->max) rms->max = p[2];
+ if (p[2] < rms->min) rms->min = p[2];
+ if (rms->relative)
+ p[2] -= cell_value (ftt_cell_parent (rms->p->cell), rms->p->t, *((FttVector *) p));
+ p[0] = (p[0] - rms->p->c.x)/rms->p->h; p[1] = (p[1] - rms->p->c.y)/rms->p->h;
+ rms->m[0][1] += p[0]; rms->m[0][2] += p[1]; rms->m[0][3] += p[0]*p[1];
+ rms->m[1][1] += p[0]*p[0]; rms->m[1][2] += p[0]*p[1]; rms->m[1][3] += p[0]*p[0]*p[1];
+ rms->m[2][2] += p[1]*p[1]; rms->m[2][3] += p[0]*p[1]*p[1];
+ rms->m[3][3] += p[0]*p[0]*p[1]*p[1];
+ rms->H[0] += p[2];
+ rms->H[1] += p[0]*p[2];
+ rms->H[2] += p[1]*p[2];
+ rms->H[3] += p[0]*p[1]*p[2];
+ rms->H[4] += p[2]*p[2];
+ rms->m[0][0] += 1.;
}
- if ((type & 0x80000000) != 0) {
- gts_file_error (fp, "geometry type should be 2D");
- PQfinish (conn);
- return FALSE;
+ return 0;
+}
+
+static void update_terrain (FttCell * cell, GfsRefineTerrain * t)
+{
+ g_assert (GFS_VALUE (cell, t->cond) == REFINED);
+ Polygon p;
+ polygon_init (&p, cell, t);
+ RMS rms;
+ rms_init (&rms, &p, TRUE);
+ guint i;
+ for (i = 0; i < t->nrs; i++)
+ r_surface_query_region (t->rs[i], p.min, p.max, (RSurfaceQuery) rms_add, &rms);
+ rms_update (&rms);
+ for (i = 0; i < NM; i++)
+ GFS_VALUE (cell, t->h[i]) = rms.h[i];
+ GFS_VALUE (cell, t->he) = rms.he;
+ GFS_VALUE (cell, t->hn) = rms.m[0][0];
+ GFS_VALUE (cell, t->cond) = RAW;
+}
+
+static void function_from_parent (FttCell * cell, GfsRefineTerrain * t, gdouble h[4])
+{
+ gdouble H[4];
+ corners_from_parent (cell, t, H);
+ function_from_corners (h, H);
+}
+
+static void cell_fine_init (FttCell * parent, GfsRefineTerrain * t)
+{
+ gfs_cell_fine_init (parent, GFS_DOMAIN (gfs_object_simulation (t)));
+ FttCellChildren child;
+ ftt_cell_children (parent, &child);
+ guint i;
+ for (i = 0; i < FTT_CELLS; i++)
+ if (child.c[i]) {
+ gdouble h[NM];
+ function_from_parent (child.c[i], t, h);
+ guint j;
+ for (j = 0; j < NM; j++)
+ GFS_VALUE (child.c[i], t->h[j]) = h[j];
+ GFS_VALUE (child.c[i], t->he) = GFS_VALUE (parent, t->he);
+ GFS_VALUE (child.c[i], t->hn) = GFS_VALUE (parent, t->hn)/FTT_CELLS;
+ GFS_VALUE (child.c[i], t->cond) = NEW_CHILD;
+ }
+}
+
+static gdouble corner_value (GfsRefineTerrain * t, FttVector * p, gdouble eps, guint level)
+{
+ GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (t));
+ gdouble v = 0., w = 0.;
+ gint i, j;
+ for (i = -1; i <= 1; i += 2)
+ for (j = -1; j <= 1; j += 2) {
+ FttVector q;
+ q.x = p->x + eps*i; q.y = p->y + eps*j;
+ FttCell * cell = gfs_domain_locate (domain, q, level);
+ if (cell) {
+ if (ftt_cell_level (cell) < level)
+ return 0.;
+ else if (GFS_VALUE (cell, t->cond) == FAIR)
+ return cell_value (cell, t, q);
+ gdouble n = GFS_VALUE (cell, t->hn);
+ if (n > 0.) {
+ g_assert (GFS_VALUE (cell, t->cond) == RAW);
+ v += cell_value (cell, t, q);
+ w += 1.;
+ }
+ }
+ }
+ return w > 0 ? v/w : 0.;
+}
+
+static void update_error_estimate (FttCell * cell, GfsRefineTerrain * t, gboolean relative)
+{
+ if (GFS_VALUE (cell, t->hn) > 0.) {
+ Polygon poly;
+ polygon_init (&poly, cell, t);
+ RMS rms;
+ rms_init (&rms, &poly, relative);
+ guint i;
+ for (i = 0; i < t->nrs; i++)
+ r_surface_query_region (t->rs[i], poly.min, poly.max, (RSurfaceQuery) rms_add, &rms);
+ for (i = 0; i < NM; i++)
+ rms.h[i] = GFS_VALUE (cell, t->h[i]);
+ GFS_VALUE (cell, t->he) = rms_minimum (&rms);
}
- if (!query (conn, "CLOSE mycursor; COMMIT", PGRES_COMMAND_OK, fp))
- return FALSE;
- q = g_strjoin (NULL,
- "PREPARE get_points_in_polygon (text) AS"
- " SELECT height, ASBINARY(geom,'" DBYTE "') FROM ", t->tables,
- " WHERE geom && $1 AND ST_Contains($1,geom)"
- // " ORDER BY randid LIMIT " DBMAX,
- " LIMIT " DBMAX,
- NULL);
- if (!query (conn, q, PGRES_COMMAND_OK, fp)) {
- g_free (q);
- return FALSE;
+ else
+ GFS_VALUE (cell, t->he) = 0.;
+}
+
+static void remove_knots (FttCell * cell, GfsRefineTerrain * t)
+{
+ gdouble size = ftt_cell_size (cell), eps = size/1000.;
+ guint level = ftt_cell_level (cell);
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ gdouble h[4], H[4];
+ p.x += size/2.; p.y += size/2.;
+ H[0] = corner_value (t, &p, eps, level);
+ p.x -= size;
+ H[1] = corner_value (t, &p, eps, level);
+ p.y -= size;
+ H[2] = corner_value (t, &p, eps, level);
+ p.x += size;
+ H[3] = corner_value (t, &p, eps, level);
+ function_from_corners (h, H);
+ guint i;
+ for (i = 0; i < NM; i++)
+ GFS_VALUE (cell, t->h[i]) = h[i];
+ GFS_VALUE (cell, t->cond) = FAIR;
+
+ update_error_estimate (cell, t, TRUE);
+}
+
+static void update_height_and_check_for_refinement (FttCell * cell, GfsRefineTerrain * t)
+{
+ if (GFS_VALUE (cell, t->cond) == FAIR) {
+ gdouble h[4];
+ function_from_parent (cell, t, h);
+ guint i;
+ for (i = 0; i < NM; i++)
+ GFS_VALUE (cell, t->h[i]) += h[i];
+
+ if (ftt_cell_level (cell) < gfs_function_value (GFS_REFINE (t)->maxlevel, cell) &&
+ gfs_function_value (t->criterion, cell)) {
+ g_assert (FTT_CELL_IS_LEAF (cell));
+ ftt_cell_refine_single (cell, (FttCellInitFunc) cell_fine_init, t);
+ FttCellChildren child;
+ guint i;
+ ftt_cell_children (cell, &child);
+ for (i = 0; i < FTT_CELLS; i++)
+ GFS_VALUE (child.c[i], t->cond) = REFINED;
+ t->refined = TRUE;
+ }
}
- g_free (q);
- t->conn = conn;
- return TRUE;
+ else
+ g_assert (GFS_VALUE (cell, t->cond) == NEW_CHILD);
+}
+
+#if DEBUG
+static void draw_terrain (FttCell * cell, gpointer * data)
+{
+ GfsRefineTerrain * t = data[0];
+ FILE * fp = data[1];
+ gdouble h = ftt_cell_size (cell);
+ FttVector p;
+ ftt_cell_pos (cell, &p);
+ p.x += h/2.; p.y += h/2.;
+ fprintf (fp, "%g %g %g\n", p.x, p.y, cell_value (cell, t, p)/20000.);
+ p.x -= h;
+ fprintf (fp, "%g %g %g\n", p.x, p.y, cell_value (cell, t, p)/20000.);
+ p.y -= h;
+ fprintf (fp, "%g %g %g\n", p.x, p.y, cell_value (cell, t, p)/20000.);
+ p.x += h;
+ fprintf (fp, "%g %g %g\n", p.x, p.y, cell_value (cell, t, p)/20000.);
+}
+
+static void draw_level (GfsDomain * domain, GfsRefine * refine, guint level, const gchar * name)
+{
+ gpointer data[2];
+ data[0] = refine;
+ data[1] = fopen (name, "w");
+ fprintf (data[1], "QUAD\n");
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, level,
+ (FttCellTraverseFunc) draw_terrain, data);
+ fclose (data[1]);
+}
+
+static void draw_all (GfsDomain * domain, GfsRefine * refine, const gchar * name)
+{
+ gpointer data[2];
+ data[0] = refine;
+ data[1] = fopen (name, "w");
+ fprintf (data[1], "QUAD\n");
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,
+ (FttCellTraverseFunc) draw_terrain, data);
+ fclose (data[1]);
+}
+#endif
+
+static void reset_terrain (FttCell * cell, GfsRefineTerrain * t)
+{
+ guint i;
+ for (i = 0; i < NM; i++)
+ GFS_VALUE (cell, t->h[i]) = 0.;
+ GFS_VALUE (cell, t->cond) = REFINED;
+ if (FTT_CELL_IS_LEAF (cell) && ftt_cell_level (cell) < t->level)
+ t->level = ftt_cell_level (cell);
+}
+
+static void terrain_refine (GfsRefine * refine, GfsSimulation * sim)
+{
+ GfsDomain * domain = GFS_DOMAIN (sim);
+ GfsRefineTerrain * t = GFS_REFINE_TERRAIN (refine);
+ t->level = G_MAXINT/2;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_ALL, -1,
+ (FttCellTraverseFunc) reset_terrain, refine);
+ do {
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, t->level,
+ (FttCellTraverseFunc) update_terrain, refine);
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, t->level,
+ (FttCellTraverseFunc) remove_knots, refine);
+ t->refined = FALSE;
+ gfs_domain_cell_traverse (domain, FTT_PRE_ORDER, FTT_TRAVERSE_LEVEL, t->level,
+ (FttCellTraverseFunc) update_height_and_check_for_refinement,
+ refine);
+#if DEBUG
+ GfsNorm norm = gfs_domain_norm_variable (domain, t->he, NULL, FTT_TRAVERSE_LEAFS, -1);
+ fprintf (stderr, "level: %d bias: %g 1: %g 2: %g inf: %g\n",
+ t->level, norm.bias, norm.first, norm.second, norm.infty);
+ fprintf (stderr, "level: %d depth: %d\n", t->level, gfs_domain_depth (domain));
+ gchar name[] = "/tmp/level-x";
+ name[11] = 48 + t->level;
+ draw_level (domain, refine, t->level, name);
+#endif
+ t->level++;
+ } while (t->refined);
+#if DEBUG
+ draw_all (domain, refine, "/tmp/all");
+#endif
+}
+
+static void refine_terrain_destroy (GtsObject * object)
+{
+ GfsRefineTerrain * t = GFS_REFINE_TERRAIN (object);
+ g_free (t->name);
+ g_free (t->basename);
+ if (t->rs) {
+ guint i;
+ for (i = 0; i < t->nrs; i++)
+ r_surface_close (t->rs[i]);
+ g_free (t->rs);
+ }
+ gts_object_destroy (GTS_OBJECT (t->criterion));
+ (* GTS_OBJECT_CLASS (gfs_refine_terrain_class ())->parent_class->destroy) (object);
}
static void refine_terrain_read (GtsObject ** o, GtsFile * fp)
@@ -360,56 +585,113 @@ static void refine_terrain_read (GtsObject ** o, GtsFile * fp)
t->name = g_strdup (fp->token->str);
gts_file_next_token (fp);
+ gchar * dir = NULL;
if (fp->type == '{') {
GtsFileVariable var[] = {
- {GTS_STRING, "host", TRUE},
- {GTS_STRING, "dbname", TRUE},
- {GTS_STRING, "user", TRUE},
- {GTS_STRING, "password", TRUE},
- {GTS_STRING, "tables", TRUE},
+ {GTS_STRING, "basename", TRUE},
+ {GTS_STRING, "dir", TRUE},
{GTS_NONE}
};
- gchar * host = NULL, * dbname = NULL, * user = NULL, * passwd = NULL, * tables = NULL;
- var[0].data = &host;
- var[1].data = &dbname;
- var[2].data = &user;
- var[3].data = &passwd;
- var[4].data = &tables;
+ gchar * basename = NULL;
+ var[0].data = &basename;
+ var[1].data = &dir;
gts_file_assign_variables (fp, var);
if (fp->type == GTS_ERROR)
return;
- if (var[0].set) { g_free (t->host); t->host = host; }
- if (var[1].set) { g_free (t->dbname); t->dbname = dbname; }
- if (var[2].set) { g_free (t->user); t->user = user; }
- if (var[3].set) { g_free (t->passwd); t->passwd = passwd; }
- if (var[4].set) { g_free (t->tables); t->tables = tables; }
+ if (var[0].set) { g_free (t->basename); t->basename = basename; }
+ if (!var[1].set) dir = g_strdup (default_dir);
}
-
- if (!postgres_connection (t, fp))
- return;
+ else
+ dir = g_strdup (default_dir);
+
+ if (!strcmp (t->basename, "*")) { /* file globbing */
+ gchar * pattern = g_strconcat (dir, "/*.Data", NULL);
+ glob_t pglob;
+ if (glob (pattern, GLOB_ERR, NULL, &pglob)) {
+ gts_file_error (fp, "cannot find/open RSurface files in path:\n%s", pattern);
+ g_free (pattern);
+ g_free (dir);
+ return;
+ }
+ g_free (pattern);
+ g_free (t->basename);
+ t->basename = NULL;
+ guint i;
+ for (i = 0; i < pglob.gl_pathc; i++) {
+ pglob.gl_pathv[i][strlen (pglob.gl_pathv[i]) - 5] = '\0';
+ t->rs = g_realloc (t->rs, (t->nrs + 1)*sizeof (RSurface *));
+ t->rs[t->nrs] = r_surface_open (pglob.gl_pathv[i], "r", -1);
+ if (!t->rs[t->nrs]) {
+ gts_file_error (fp, "cannot open RSurface `%s'", pglob.gl_pathv[i]);
+ globfree (&pglob);
+ g_free (dir);
+ return;
+ }
+ if (t->basename) {
+ gchar * dirbasename = g_strconcat (t->basename, ",", pglob.gl_pathv[i], NULL);
+ g_free (t->basename);
+ t->basename = dirbasename;
+ }
+ else
+ t->basename = g_strdup (pglob.gl_pathv[i]);
+ t->nrs++;
+ }
+ globfree (&pglob);
+ }
+ else { /* basename is of the form: set1,set2,set3... */
+ gchar * basename = g_strdup (t->basename);
+ if (dir) {
+ g_free (t->basename);
+ t->basename = NULL;
+ }
+ gchar * s = strtok (basename, ",");
+ while (s) {
+ t->rs = g_realloc (t->rs, (t->nrs + 1)*sizeof (RSurface *));
+ if (dir) {
+ gchar * fname = s[0] == '/' ? g_strdup (s) : g_strconcat (dir, "/", s, NULL);
+ t->rs[t->nrs] = r_surface_open (fname, "r", -1);
+ if (t->basename) {
+ gchar * dirbasename = g_strconcat (t->basename, ",", fname, NULL);
+ g_free (t->basename);
+ t->basename = dirbasename;
+ g_free (fname);
+ }
+ else
+ t->basename = fname;
+ }
+ else
+ t->rs[t->nrs] = r_surface_open (s, "r", -1);
+ if (!t->rs[t->nrs]) {
+ gts_file_error (fp, "cannot open RSurface `%s'", s);
+ g_free (basename);
+ g_free (dir);
+ return;
+ }
+ t->nrs++;
+ s = strtok (NULL, ",");
+ }
+ g_free (basename);
+ }
+ g_free (dir);
GfsDomain * domain = GFS_DOMAIN (gfs_object_simulation (*o));
- gchar * name = g_strjoin (NULL, t->name, "0", NULL);
- t->h = gfs_domain_get_or_add_variable (domain, name, "Terrain height");
- g_free (name);
- name = g_strjoin (NULL, t->name, "x", NULL);
- t->hx = gfs_domain_get_or_add_variable (domain, name, "Terrain slope");
- g_free (name);
- name = g_strjoin (NULL, t->name, "y", NULL);
- t->hy = gfs_domain_get_or_add_variable (domain, name, "Terrain slope");
- g_free (name);
- name = g_strjoin (NULL, t->name, "e", NULL);
+ guint i;
+ for (i = 0; i < NM; i++) {
+ gchar * name = g_strdup_printf ("%s%d", t->name, i);
+ t->h[i] = gfs_domain_get_or_add_variable (domain, name, "Terrain height");
+ g_free (name);
+ }
+ gchar * name = g_strjoin (NULL, t->name, "e", NULL);
t->he = gfs_domain_get_or_add_variable (domain, name, "Terrain RMS error");
g_free (name);
name = g_strjoin (NULL, t->name, "n", NULL);
t->hn = gfs_domain_get_or_add_variable (domain, name, "Terrain samples #");
g_free (name);
+ name = g_strjoin (NULL, t->name, "c", NULL);
+ t->cond = gfs_domain_get_or_add_variable (domain, name, "Terrain condition number");
+ g_free (name);
gfs_function_read (t->criterion, domain, fp);
- if (fp->type == GTS_ERROR) {
- PQfinish (t->conn);
- t->conn = NULL;
- }
}
static void refine_terrain_write (GtsObject * o, FILE * fp)
@@ -418,8 +700,7 @@ static void refine_terrain_write (GtsObject * o, FILE * fp)
(* GTS_OBJECT_CLASS (gfs_refine_terrain_class ())->parent_class->write) (o, fp);
fprintf (fp, " %s ", t->name);
gfs_function_write (t->criterion, fp);
- fprintf (fp, " { host = %s dbname = %s user = %s tables = %s }",
- t->host, t->dbname, t->user, t->tables);
+ fprintf (fp, " { basename = %s }", t->basename);
}
static void gfs_refine_terrain_class_init (GfsRefineClass * klass)
@@ -434,11 +715,7 @@ static void gfs_refine_terrain_class_init (GfsRefineClass * klass)
static void gfs_refine_terrain_init (GfsRefineTerrain * t)
{
t->criterion = gfs_function_new (gfs_function_class (), 0.);
- t->host = g_strdup ("localhost");
- t->dbname = g_strdup ("earth");
- t->user = g_strdup ("popinet");
- t->passwd = g_strdup ("CloClo");
- t->tables = g_strdup ("etopo2");
+ t->basename = g_strdup ("*");
}
GfsRefineClass * gfs_refine_terrain_class (void)
@@ -470,6 +747,9 @@ void gfs_module_write (FILE * fp);
const gchar * g_module_check_init (void)
{
+ gchar * dir = getenv ("GFS_RSURFACES_PATH");
+ if (dir)
+ default_dir = dir;
gfs_refine_terrain_class ();
return NULL;
}
diff --git a/modules/xyz2rsurface.c b/modules/xyz2rsurface.c
new file mode 100644
index 0000000..ede3e8a
--- /dev/null
+++ b/modules/xyz2rsurface.c
@@ -0,0 +1,30 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "rsurface.h"
+
+int main (int argc, char** argv)
+{
+ if (argc != 3) {
+ fprintf (stderr, "Usage: %s basename pagesize\n", argv[0]);
+ return -1;
+ }
+
+ RSurface * rs = r_surface_open (argv[1], "w", atoi (argv[2]));
+ if (rs == NULL) {
+ fprintf (stderr, "xyz2rsurface: cannot open files `%s*'\n", argv[1]);
+ return 1;
+ }
+ double x[3];
+ int id, count = 0;
+ while (scanf ("%d %lf %lf %lf", &id, &x[0], &x[1], &x[2]) == 4) {
+ if (!r_surface_insert (rs, x, id)) {
+ fprintf (stderr, "\nxyz2rsurface: error inserting point (%g,%g,%g)\n",
+ x[0], x[1], x[2]);
+ return 1;
+ }
+ count++;
+ }
+ r_surface_close (rs);
+
+ return 0.;
+}
diff --git a/src/utils.c b/src/utils.c
index fa636b2..1c715f3 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -1302,18 +1302,18 @@ void gfs_eigenvalues (gdouble a[FTT_DIMENSION][FTT_DIMENSION],
*
* Replaces @m with its inverse.
*
- * Returns: %FALSE if the inversion encounters a pivot coefficient
- * smaller than or equal to @pivmin (i.e. @m is non-invertible), %TRUE
- * otherwise.
+ * Returns: 0. if the inversion encounters a pivot coefficient smaller
+ * than or equal to @pivmin (i.e. @m is non-invertible), the minimum
+ * absolute value of the pivoting coefficient otherwise.
*/
-gboolean gfs_matrix_inverse (gdouble ** m, guint n, gdouble pivmin)
+gdouble gfs_matrix_inverse (gdouble ** m, guint n, gdouble pivmin)
{
gint * indxc, * indxr, * ipiv;
gint i, icol = 0, irow = 0, j, k, l, ll;
- gdouble big, dum, pivinv, temp;
+ gdouble big, dum, pivinv, temp, minpiv = G_MAXDOUBLE;
- g_return_val_if_fail (m != NULL, FALSE);
- g_return_val_if_fail (pivmin >= 0., FALSE);
+ g_return_val_if_fail (m != NULL, 0.);
+ g_return_val_if_fail (pivmin >= 0., 0.);
indxc = g_malloc (sizeof (gint)*n);
indxr = g_malloc (sizeof (gint)*n);
@@ -1347,8 +1347,10 @@ gboolean gfs_matrix_inverse (gdouble ** m, guint n, gdouble pivmin)
g_free (indxc);
g_free (indxr);
g_free (ipiv);
- return FALSE;
+ return 0.;
}
+ if (fabs (m[icol][icol]) < minpiv)
+ minpiv = fabs (m[icol][icol]);
pivinv = 1.0/m[icol][icol];
m[icol][icol] = 1.0;
for (l = 0; l < n; l++) m[icol][l] *= pivinv;
@@ -1368,7 +1370,7 @@ gboolean gfs_matrix_inverse (gdouble ** m, guint n, gdouble pivmin)
g_free (indxc);
g_free (indxr);
g_free (ipiv);
- return TRUE;
+ return minpiv;
}
/**
diff --git a/src/utils.h b/src/utils.h
index 60d7230..06017c6 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -111,7 +111,7 @@ GtsObjectClass * gfs_object_class_from_name (const gchar * name);
void gfs_eigenvalues (gdouble a[FTT_DIMENSION][FTT_DIMENSION],
gdouble d[FTT_DIMENSION],
gdouble v[FTT_DIMENSION][FTT_DIMENSION]);
-gboolean gfs_matrix_inverse (gdouble ** m,
+gdouble gfs_matrix_inverse (gdouble ** m,
guint n,
gdouble pivmin);
gpointer gfs_matrix_new (guint n,
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list