[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:55:46 UTC 2009
The following commit has been merged in the upstream branch:
commit 67b35b790f4935d2b8f0d7b9c4fdba875ac73584
Author: Stephane Popinet <popinet at users.sf.net>
Date: Thu Oct 2 11:53:28 2008 +1000
New function r_surface_query_region_sum()
Efficiently computes the statistics for a given region using the hierarchical
statistics stored in the R*-tree.
darcs-hash:20081002015328-d4795-b0a68fccbe1d91b3e649692391385d9037c74482.gz
diff --git a/modules/RStarTree/RSTBase.h b/modules/RStarTree/RSTBase.h
index d4cb113..9a44dd3 100644
--- a/modules/RStarTree/RSTBase.h
+++ b/modules/RStarTree/RSTBase.h
@@ -60,8 +60,9 @@ typedef enum {
} Side;
typedef struct {
- typrect rect;
- int ptrtosub;
+ typrect rect;
+ int ptrtosub;
+ typdirinfo info;
} typDIRent, *refDIRent; /* inner entry */
typedef struct {
typrect rect;
diff --git a/modules/RStarTree/RSTInterUtil.c b/modules/RStarTree/RSTInterUtil.c
index f23d52e..9e68ff0 100644
--- a/modules/RStarTree/RSTInterUtil.c
+++ b/modules/RStarTree/RSTInterUtil.c
@@ -213,6 +213,7 @@ void SetCheckDir(RSTREE R, boolean creation)
if (creation) {
(*par).direntrylen= sizeof(typDIRent); /* set */
PACKEDdirentrylen= sizeof(typrect) + sizeof(int);
+#if 0
if (PACKEDdirentrylen != (*par).direntrylen) {
printf("\n%s\n"," ----- WARNING -----");
printf("%s\n","Directory entries are not packed!");
@@ -220,6 +221,7 @@ void SetCheckDir(RSTREE R, boolean creation)
printf("%s %d\n"," Real space needed:",(*par).direntrylen);
printf("%s\n",/* implicitly */"Applying the latter!");
}
+#endif
}
SIZE_DIRnodeOf3= sizeof(typDIRnodeOf3);
diff --git a/modules/RStarTree/RSTQuery.c b/modules/RStarTree/RSTQuery.c
index 5ebb949..eb99d40 100644
--- a/modules/RStarTree/RSTQuery.c
+++ b/modules/RStarTree/RSTQuery.c
@@ -338,3 +338,152 @@ void All(RSTREE R,
}
/************************************************************************/
+
+static void dirinfo_init (typdirinfo * info)
+{
+ memset (info, 0, sizeof (typdirinfo));
+ info->Hmin = 1e30;
+ info->Hmax = - 1e30;
+}
+
+static void dirinfo_add_dir (typdirinfo * info, typdirinfo * a)
+{
+ info->m01 += a->m01;
+ info->m02 += a->m02;
+ info->m03 += a->m03;
+ info->m11 += a->m11;
+ info->m13 += a->m13;
+ info->m22 += a->m22;
+ info->m23 += a->m23;
+ info->m33 += a->m33;
+ info->H0 += a->H0;
+ info->H1 += a->H1;
+ info->H2 += a->H2;
+ info->H3 += a->H3;
+ info->H4 += a->H4;
+ info->n += a->n;
+ if (a->Hmin < info->Hmin) info->Hmin = a->Hmin;
+ if (a->Hmax > info->Hmax) info->Hmax = a->Hmax;
+}
+
+static void dirinfo_add_data (typdirinfo * info, typrect rect, refinfo data)
+{
+ double p[3];
+ p[0] = rect[0].l; p[1] = rect[1].l; p[2] = data->height;
+ info->m01 += p[0];
+ info->m02 += p[1];
+ info->m03 += p[0]*p[1];
+ info->m11 += p[0]*p[0];
+ info->m13 += p[0]*p[0]*p[1];
+ info->m22 += p[1]*p[1];
+ info->m23 += p[0]*p[1]*p[1];
+ info->m33 += p[0]*p[0]*p[1]*p[1];
+ info->H0 += p[2];
+ info->H1 += p[0]*p[2];
+ info->H2 += p[1]*p[2];
+ info->H3 += p[0]*p[1]*p[2];
+ info->H4 += p[2]*p[2];
+ info->n++;
+ if (p[2] < info->Hmin) info->Hmin = p[2];
+ if (p[2] > info->Hmax) info->Hmax = p[2];
+}
+
+void UpdateAll(RSTREE R,
+ int depth,
+ typdirinfo * info)
+
+{
+ refDIRnode DIN;
+ refDATAnode DAN;
+ typrect rectfound;
+ int i;
+
+ if (depth != (*R).parameters._.height) {
+
+ DIN= &(*(*R).N[depth]).DIR;
+
+ for (i= 0; i < (*DIN).nofentries; i++) {
+ (*R).E[depth]= i;
+ if ((*DIN).entries[i].ptrtosub != (*R).P[depth+1]) {
+ NewNode(R,depth+1);
+ }
+ dirinfo_init (&(*DIN).entries[i].info);
+ (*R).Nmodified[depth] = 1;
+ UpdateAll(R,depth+1,&(*DIN).entries[i].info);
+ dirinfo_add_dir (info, &(*DIN).entries[i].info);
+ }
+ }
+ else {
+
+ DAN= &(*(*R).N[depth]).DATA;
+
+ for (i= 0; i < (*DAN).nofentries; i++) {
+ (*R).E[depth]= i;
+
+ CopyRect(R,(*DAN).entries[i].rect,rectfound); /* avoid modification */
+ dirinfo_add_data (info, rectfound, &(*DAN).entries[i].info);
+ }
+ }
+}
+
+/************************************************************************/
+
+void RgnQueryInfo(RSTREE R,
+ int depth,
+ Check includes,
+ Check intersects,
+ void * data,
+ typdirinfo * info)
+
+{
+ refcount c;
+ refDIRnode DIN;
+ refDATAnode DAN;
+ boolean istoread;
+ typrect rectfound;
+ int i;
+
+ if (depth != (*R).parameters._.height) {
+
+ c= &(*R).count;
+ if ((*c).countflag) {
+ (*c).dirvisitcount++;
+ }
+
+ DIN= &(*(*R).N[depth]).DIR;
+
+ for (i= 0; i < (*DIN).nofentries; i++) {
+ if ((* includes) ((*DIN).entries[i].rect,data))
+ dirinfo_add_dir (info, &(*DIN).entries[i].info);
+ else if ((* intersects) ((*DIN).entries[i].rect,data)) {
+ (*R).E[depth]= i;
+ istoread= (*DIN).entries[i].ptrtosub != (*R).P[depth+1];
+ if (istoread) {
+ NewNode(R,depth+1);
+ }
+ RgnQueryInfo(R,depth+1,includes,intersects,data,info);
+ }
+ }
+ }
+ else {
+
+ c= &(*R).count;
+ if ((*c).countflag) {
+ (*c).datavisitcount++;
+ }
+
+ DAN= &(*(*R).N[depth]).DATA;
+
+ for (i= 0; i < (*DAN).nofentries; i++) {
+
+ if ((* includes) ((*DAN).entries[i].rect,data)) {
+ (*R).E[depth]= i;
+
+ CopyRect(R,(*DAN).entries[i].rect,rectfound); /* avoid modification */
+ dirinfo_add_data (info, rectfound, &(*DAN).entries[i].info);
+ }
+ }
+ }
+}
+
+/************************************************************************/
diff --git a/modules/RStarTree/RSTQuery.h b/modules/RStarTree/RSTQuery.h
index dfa1d72..53ba483 100644
--- a/modules/RStarTree/RSTQuery.h
+++ b/modules/RStarTree/RSTQuery.h
@@ -42,13 +42,22 @@ void RgnQuery(RSTREE R,
void *buf,
boolean *finish);
-
void All(RSTREE R,
int depth,
QueryManageProc Manage,
void *buf,
boolean *finish);
+void UpdateAll(RSTREE R,
+ int depth,
+ typdirinfo * info);
+
+void RgnQueryInfo(RSTREE R,
+ int depth,
+ Check includes,
+ Check intersects,
+ void * data,
+ typdirinfo * info);
#endif /* !__RSTQuery_h */
diff --git a/modules/RStarTree/RStarTree.c b/modules/RStarTree/RStarTree.c
index 5fc4bc6..ac33da8 100644
--- a/modules/RStarTree/RStarTree.c
+++ b/modules/RStarTree/RStarTree.c
@@ -303,13 +303,12 @@ boolean ExistsRegion(RSTREE R,
boolean *regionfound)
{
- int i;
-
if (R == NULL) {
*regionfound= FALSE;
return FALSE;
}
- /**/
+ /*
+ int i;
for (i= 2; i <= (*R).parameters._.height; i++) {
if ((*R).Nmodified[i]) {
PutNode(R,(*R).N[i],(*R).P[i],i);
@@ -317,7 +316,7 @@ boolean ExistsRegion(RSTREE R,
}
(*R).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
(*R).RSTDone= TRUE;
@@ -336,13 +335,12 @@ boolean RegionCount(RSTREE R,
int *recordcount)
{
- int i;
-
if (R == NULL) {
*recordcount= 0;
return FALSE;
}
- /**/
+ /*
+ int i;
for (i= 2; i <= (*R).parameters._.height; i++) {
if ((*R).Nmodified[i]) {
PutNode(R,(*R).N[i],(*R).P[i],i);
@@ -350,7 +348,7 @@ boolean RegionCount(RSTREE R,
}
(*R).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
(*R).RSTDone= TRUE;
@@ -371,12 +369,12 @@ boolean RegionQuery(RSTREE R,
{
boolean finish;
- int i;
if (R == NULL) {
return FALSE;
}
- /**/
+ /*
+ int i;
for (i= 2; i <= (*R).parameters._.height; i++) {
if ((*R).Nmodified[i]) {
PutNode(R,(*R).N[i],(*R).P[i],i);
@@ -384,7 +382,7 @@ boolean RegionQuery(RSTREE R,
}
(*R).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
(*R).RSTDone= TRUE;
@@ -396,18 +394,47 @@ boolean RegionQuery(RSTREE R,
/************************************************************************/
+boolean RegionQueryInfo(RSTREE R,
+ Check includes,
+ Check intersects,
+ void * data,
+ typdirinfo * info)
+
+{
+ if (R == NULL) {
+ return FALSE;
+ }
+ /*
+ int i;
+ for (i= 2; i <= (*R).parameters._.height; i++) {
+ if ((*R).Nmodified[i]) {
+ PutNode(R,(*R).N[i],(*R).P[i],i);
+ (*R).Nmodified[i]= FALSE;
+ }
+ (*R).P[i]= 0;
+ }
+ *//* to be inserted, if main memory path shall be initialized
+ for test purpose */
+
+ (*R).RSTDone= TRUE;
+ RgnQueryInfo(R,1,includes,intersects,data,info);
+ return (*R).RSTDone;
+}
+
+/************************************************************************/
+
boolean AllQuery(RSTREE R,
QueryManageProc ManageProc,
void *buf)
{
boolean finish;
- int i;
if (R == NULL) {
return FALSE;
}
- /**/
+ /*
+ int i;
for (i= 2; i <= (*R).parameters._.height; i++) {
if ((*R).Nmodified[i]) {
PutNode(R,(*R).N[i],(*R).P[i],i);
@@ -415,7 +442,7 @@ boolean AllQuery(RSTREE R,
}
(*R).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
(*R).RSTDone= TRUE;
@@ -423,6 +450,34 @@ boolean AllQuery(RSTREE R,
All(R,1,ManageProc,buf,&finish);
return (*R).RSTDone;
}
+
+/************************************************************************/
+
+boolean Update(RSTREE R)
+
+{
+ typdirinfo info;
+
+ if (R == NULL) {
+ return FALSE;
+ }
+ /*
+ int i;
+ for (i= 2; i <= (*R).parameters._.height; i++) {
+ if ((*R).Nmodified[i]) {
+ PutNode(R,(*R).N[i],(*R).P[i],i);
+ (*R).Nmodified[i]= FALSE;
+ }
+ (*R).P[i]= 0;
+ }
+ *//* to be inserted, if main memory path shall be initialized
+ for test purpose */
+
+ (*R).RSTDone= TRUE;
+ UpdateAll(R,1,&info);
+ return (*R).RSTDone;
+}
+
/************************************************************************/
boolean JoinCountNv(RSTREE R1, RSTREE R2,
@@ -463,7 +518,7 @@ boolean JoinCountNv(RSTREE R1, RSTREE R2,
abort();
}
}
- /**/
+ /*
for (i= 2; i <= (*R1).parameters._.height; i++) {
if ((*R1).Nmodified[i]) {
PutNode(R1,(*R1).N[i],(*R1).P[i],i);
@@ -478,7 +533,7 @@ boolean JoinCountNv(RSTREE R1, RSTREE R2,
}
(*R2).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
if (! ((*R1).RSTDone && (*R2).RSTDone) ) {
@@ -541,7 +596,7 @@ boolean JoinNv(RSTREE R1, RSTREE R2,
abort();
}
}
- /**/
+ /*
for (i= 2; i <= (*R1).parameters._.height; i++) {
if ((*R1).Nmodified[i]) {
PutNode(R1,(*R1).N[i],(*R1).P[i],i);
@@ -556,7 +611,7 @@ boolean JoinNv(RSTREE R1, RSTREE R2,
}
(*R2).P[i]= 0;
}
- /**//* to be inserted, if main memory path shall be initialized
+ *//* to be inserted, if main memory path shall be initialized
for test purpose */
if (! ((*R1).RSTDone && (*R2).RSTDone) ) {
diff --git a/modules/RStarTree/RStarTree.h b/modules/RStarTree/RStarTree.h
index eb55ef0..0c6dcd3 100644
--- a/modules/RStarTree/RStarTree.h
+++ b/modules/RStarTree/RStarTree.h
@@ -83,6 +83,19 @@ typedef struct {
associated with a data record.
RESTRICTION: sizeof(typinfo) >= sizeof(int) must hold! */
+typedef struct {
+ double m01, m02, m03;
+ double m11, m13;
+ double m22, m23, m33;
+ double H0, H1, H2, H3, H4;
+ float Hmin, Hmax;
+ int n;
+} typdirinfo;
+
+/* A typdirinfo is a struct which may contain arbitrary information
+ associated with a directory record. */
+
+typedef int (* Check) (typrect rect, void * data);
/* ------------------------- private includes -------------------------- */
@@ -432,6 +445,11 @@ boolean RegionQuery(RSTREE rst,
pointer:
Arbitrary pointer passed through to the procedure Manage. */
+boolean RegionQueryInfo(RSTREE R,
+ Check includes,
+ Check intersects,
+ void * data,
+ typdirinfo * info);
boolean AllQuery(RSTREE rst,
QueryManageProc Manage,
@@ -459,6 +477,9 @@ boolean AllQuery(RSTREE rst,
See also RegionQuery. */
+boolean Update(RSTREE rst);
+
+ /* Updates the directory nodes typdirinfo. */
boolean JoinCountNv(RSTREE rst1,
RSTREE rst2,
diff --git a/modules/rsurface.c b/modules/rsurface.c
index f26d41a..e98b10a 100644
--- a/modules/rsurface.c
+++ b/modules/rsurface.c
@@ -106,7 +106,28 @@ void r_surface_query_region (RSurface * rt,
RegionQuery (rt->t, rect, unused, Intersects, Intersects, ManageQuery, data);
}
+void r_surface_sum_init (RSurfaceSum * sum)
+{
+ memset (sum, 0, sizeof (RSurfaceSum));
+ sum->Hmin = 1e30;
+ sum->Hmax = - 1e30;
+}
+
+void r_surface_query_region_sum (RSurface * rt,
+ RSurfaceCheck includes,
+ RSurfaceCheck intersects,
+ void * data,
+ RSurfaceSum * sum)
+{
+ RegionQueryInfo (rt->t, (Check) includes, (Check) intersects, data, (typdirinfo *) sum);
+}
+
const char * r_surface_name (RSurface * rt)
{
return rt->name;
}
+
+void r_surface_update (RSurface * rt)
+{
+ Update (rt->t);
+}
diff --git a/modules/rsurface.h b/modules/rsurface.h
index 0112937..5746fdb 100644
--- a/modules/rsurface.h
+++ b/modules/rsurface.h
@@ -1,7 +1,25 @@
typedef struct _RSurface RSurface;
+typedef struct { /* needs to be identical to typdirinfo in RStarTree.h */
+ double m01, m02, m03;
+ double m11, m13;
+ double m22, m23, m33;
+ double H0, H1, H2, H3, H4;
+ float Hmin, Hmax;
+ int n;
+} RSurfaceSum;
+
+typedef struct {
+ float l, h;
+} RSurfaceInterval;
+
+typedef RSurfaceInterval RSurfaceRect[2];
+
+typedef int (* RSurfaceCheck) (RSurfaceRect rect, void * data);
+
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_update (RSurface * rt);
void r_surface_close (RSurface * rt);
int r_surface_insert (RSurface * rt, double p[3], int id);
@@ -12,4 +30,10 @@ void r_surface_query (RSurface * rt,
void r_surface_query_region (RSurface * rt,
double min[2], double max[2],
RSurfaceQuery q, void * user_data);
+void r_surface_sum_init (RSurfaceSum * sum);
+void r_surface_query_region_sum (RSurface * rt,
+ RSurfaceCheck includes,
+ RSurfaceCheck intersects,
+ void * data,
+ RSurfaceSum * sum);
const char * r_surface_name (RSurface * rt);
diff --git a/modules/xyz2rsurface.c b/modules/xyz2rsurface.c
index c321cd6..79c54e0 100644
--- a/modules/xyz2rsurface.c
+++ b/modules/xyz2rsurface.c
@@ -100,6 +100,7 @@ int main (int argc, char** argv)
if (verbose && (count % 1000) == 0)
fprintf (stderr, "\rxyz2rsurface: %9d points inserted", count);
}
+ r_surface_update (rs);
r_surface_close (rs);
if (verbose)
fprintf (stderr, "\rxyz2rsurface: %9d points inserted\n", count);
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list