Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12217) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12218) @@ -572,6 +572,6 @@ ./objects/Bamg/Metric.cpp\ ./objects/Bamg/Metric.h\ - ./objects/Bamg/QuadTree.cpp\ - ./objects/Bamg/QuadTree.h\ + ./objects/Bamg/BamgQuadtree.cpp\ + ./objects/Bamg/BamgQuadtree.h\ ./objects/Bamg/R2.h\ ./objects/Bamg/SetOfE4.cpp\ Index: /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 12217) +++ /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 12218) @@ -92,5 +92,5 @@ //Make Quadtree from background mesh - BTh.MakeQuadTree(); + BTh.MakeBamgQuadtree(); //Bound hmin and hmax Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.cpp (revision 12218) +++ /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.cpp (revision 12218) @@ -0,0 +1,598 @@ +#include +#include +#include + +#include "../objects.h" + +namespace bamg { + + /*MACROS {{{1*/ + /* + * + * J j + * ^ ^ + * | | +--------+--------+ + * | | | | | + * 1X | | | 2 | 3 | + * | | | | | + * | | +--------+--------+ + * | | | | | + * 0X | | | 0 | 1 | + * | | | | | + * | | +--------+--------+ + * | +-----------------------> i + * | + * |----------------------------> I + * X0 X1 + * + * box 0 -> I=0 J=0 IJ=00 = 0 + * box 1 -> I=1 J=0 IJ=01 = 1 + * box 2 -> I=0 J=1 IJ=10 = 2 + * box 3 -> I=1 J=1 IJ=11 = 3 + */ +#define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b))) +#define ABS(i) ((i)<0 ?-(i) :(i)) +#define MAX1(i,j) ((i)>(j) ?(i) :(j)) +#define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2))) + + //IJ(i,j,l) returns the box number of i and j with respect to l + //if !j&l and !i&l -> 0 (box zero: lower left ) + //if !j&l and i&l -> 1 (box one: lower right) + //if j&l and !i&l -> 2 (box two: upper left ) + //if j&l and i&l -> 3 (box three:upper right) +#define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) + + //I_IJ(k,l) returns l if first bit of k is 1, else 0 +#define I_IJ(k,l) ((k&1) ? l:0) + //J_IJ(k,l) returns l if second bit of k is 1, else 0 +#define J_IJ(k,l) ((k&2) ? l:0) + /*}}}*/ + /*DOCUMENTATION What is a BamgQuadtree? {{{1 + * A Quadtree is a very simple way to group vertices according + * to their locations. A square that holds all the points of the mesh + * (or the geometry) is divided into 4 boxes. As soon as one box + * hold more than 4 vertices, it is divided into 4 new boxes, etc... + * There cannot be more than MAXDEEP (=30) subdivision. + * This process is like a Dichotomy in dimension 2 + * + * + - - - - - - - - + - - + - + - + - - - - + + * | | | | X | | + * + - + - + + * | | | | | | + * + - - + - + - + + + * | | | | | + * + * | | | | | + * + - - - - - - - - + - - + - - + - - - - + + * | | | | + * + * | | | | + * + * | | | | + * | | | | + * + - - - - - - - - + - - - - + - - - - + + * | | | + * + * | | | + * + * | | | + * + * | | | + * | | | + * | | | + * | | | + * | | | + * + - - - - - - - - + - - - - - - - - + + * + * The coordinate system used in a quadtree are integers to avoid + * round-off errors. The vertex in the lower left box has the coordinates + * (0 0) + * The upper right vertex has the follwing coordinates: + * 2^30 -1 2^30 -1 in decimal + * 0 1 1 1 .... 1 0 1 1 1 .... 1 in binary + * \-- 29 --/ \-- 29 --/ + * Using binaries is therefore very easy to locate a vertex in a box: + * we just need to look at the bits from the left to the right (See ::Add) + }}}1*/ + + /*Constructors/Destructors*/ + /*FUNCTION BamgQuadtree::BamgQuadtree(){{{1*/ + BamgQuadtree::BamgQuadtree(){ + + /*Number of boxes and vertices*/ + NbBamgQuadtreeBox=0; + NbVertices=0; + + /*Create container*/ + boxcontainer=new DataSet(); + + /*Create Root, pointer toward the main box*/ + root=NewBamgQuadtreeBox(); + + } + /*}}}1*/ + /*FUNCTION BamgQuadtree::BamgQuadtree(Mesh * t,long nbv){{{1*/ + BamgQuadtree::BamgQuadtree(Mesh * t,long nbv){ + + /*Number of boxes and vertices*/ + NbBamgQuadtreeBox=0; + NbVertices=0; + + /*Create container*/ + boxcontainer=new DataSet(); + + /*Create Root, pointer toward the main box*/ + root=NewBamgQuadtreeBox(); + + /*Check Sizes*/ + _assert_(MaxISize>MaxICoor); + + /*Add all vertices of the mesh*/ + if (nbv==-1) nbv=t->nbv; + for (int i=0;ivertices[i]); + + } + /*}}}1*/ + /*FUNCTION BamgQuadtree::~BamgQuadtree(){{{1*/ + BamgQuadtree::~BamgQuadtree() { + delete boxcontainer; + root=NULL; + } + /*}}}1*/ + + /*Methods*/ + /*FUNCTION BamgQuadtree::Add{{{1*/ + void BamgQuadtree::Add(BamgVertex &w){ + /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/Add)*/ + BamgQuadtreeBox** pb=NULL; + BamgQuadtreeBox* b=NULL; + + /*Get integer coodinate of current point w*/ + register long i=w.i.x, j=w.i.y; + + /*Initialize level*/ + register long level=MaxISize; + + /*Get inital box (the largest)*/ + pb = &root; + + /*Find the smallest box where w is located*/ + while((b=*pb) && (b->nbitems<0)){ + + //shift b->nbitems by -1 + b->nbitems--; + + //shifted righ by one bit: level=00000010 -> 00000001 + level >>= 1; + + //Get next subbox according to the bit value (level) + pb = &b->b[IJ(i,j,level)]; + } + + /*OK, we have found b, a Subbox holding vertices (might be full) + check that the vertex is not already in the box*/ + if (b){ + if (b->nbitems > 3 && b->v[3] == &w) return; + if (b->nbitems > 2 && b->v[2] == &w) return; + if (b->nbitems > 1 && b->v[1] == &w) return; + if (b->nbitems > 0 && b->v[0] == &w) return; + } + + /*check that l is not 0 (this should not happen as MaxDepth = 30)*/ + _assert_(level>0); + + /*Now, try to add the vertex, if the subbox is full (nbitems=4), we have to divide it + in 4 new subboxes*/ + while ((b= *pb) && (b->nbitems == 4)){ // the BamgQuadtreeBox is full + + /*Copy the 4 vertices in the current BamgQuadtreebox*/ + BamgVertex* v4[4]; + v4[0]= b->v[0]; + v4[1]= b->v[1]; + v4[2]= b->v[2]; + v4[3]= b->v[3]; + + /*set nbitems as negative + * (box full -> holds 4 pointers toward subboxes and not 4 vertices)*/ + b->nbitems = -b->nbitems; + + /*Initialize the 4 pointers toward the 4 subboxes*/ + b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL; + + /*level = 0010000 -> 0001000*/ + level >>= 1; + + /*Put the four vertices in the new boxes*/ + for (int k=0;k<4;k++){ + + int ij; + /*bb is the new "sub"box of b where v4[k] is located*/ + BamgQuadtreeBox *bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)]; + + // alloc the BamgQuadtreeBox + if (!bb) bb=b->b[ij]=NewBamgQuadtreeBox(); + + /*Copy the current vertex*/ + bb->v[bb->nbitems++] = v4[k]; + } + + /*Get the subbox where w (i,j) is located*/ + pb = &b->b[IJ(i,j,level)]; + } + + /*alloc the BamgQuadtreeBox if necessary*/ + if (!(b=*pb)) b=*pb= NewBamgQuadtreeBox(); + + /*Add w*/ + b->v[b->nbitems++]=&w; + + //Increase NbVertices by one (we have one new vertex) + NbVertices++; + } + /*}}}1*/ + /*FUNCTION BamgQuadtree::NearestVertex{{{1*/ + BamgVertex* BamgQuadtree::NearestVertex(Icoor1 i,Icoor1 j) { + /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/NearestVertex)*/ + + /*Intermediaries*/ + BamgQuadtreeBox *pb[MaxDepth]; + int pi[MaxDepth]; + Icoor1 ii[MaxDepth]; + Icoor1 jj[MaxDepth]; + int level; + long n0; + BamgQuadtreeBox *b; + long h0; + long h = MaxISize; + long hb= MaxISize; + Icoor1 i0=0,j0=0; + + /*initial output as NULL (no vertex found)*/ + BamgVertex* nearest_v=NULL; + + /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/ + Icoor1 iplus( inbitems) return nearest_v; + + /*else, find the smallest non-empty BamgQuadtreeBox containing the point (i,j)*/ + while((n0=b->nbitems)<0){ + + Icoor1 hb2 = hb >> 1; //size of the current box + int k = IJ(iplus,jplus,hb2); //box number (0,1,2 or 3) + BamgQuadtreeBox *b0 = b->b[k]; //pointer toward current box + + /* break if NULL box or empty (Keep previous box b)*/ + if (( b0 == NULL) || (b0->nbitems == 0)) break; + + /*Get next Quadtree box*/ + b=b0; + i0 += I_IJ(k,hb2); // i orign of BamgQuadtreeBox (macro) + j0 += J_IJ(k,hb2); // j orign of BamgQuadtreeBox + hb = hb2; // size of the box (in Int) + } + + /*The box b, is the smallest box containing the point (i,j) and + * has the following properties: + * - n0: number of items (>0 if vertices, else boxes) + * - hb: box size (int) + * - i0: x coordinate of the lower left corner + * - j0: y coordinate of the lower left corner*/ + + /* if the current subbox is holding vertices, we are almost done*/ + if (n0>0){ + //loop over the vertices of the box and find the closest vertex + for(int k=0;kv[k]->i; + + /*Compute norm with w*/ + h0=NORM(iplus,i2.x,jplus,i2.y); + + /*is it smaller than previous value*/ + if (h0v[k]; + } + } + /*return closest vertex*/ + return nearest_v; + } + + /* general case: the current box is empty, we have to go backwards + and find the closest not-empty box and find the closest vertex*/ + + /*Initialize search variables*/ + pb[0]=b; //pointer toward the box b + pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b + ii[0]=i0; //i coordinate of the box lowest left corner + jj[0]=j0; //j coordinate of the box lowest left corner + + /*initialize h: smallest box size, containing a vertex close to w*/ + h=hb; + + /*Main loop*/ + level=0; + do { + + /*get current box*/ + b= pb[level]; + + /*Loop over the items in current box (if not empty!)*/ + while (pi[level]){ + + /*We are looping now over the items of b. k is the current index (in [0 3])*/ + pi[level]--; + int k=pi[level]; + + /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/ + if (b->nbitems>0){ + I2 i2=b->v[k]->i; + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0v[k]; + } + } + /*else: current box b is pointing toward 4 boxes + * test sub-box k and go deeper into the tree if it is non empty + * and contains the point w modulo a size h that is either the size of the smallest + * non empty box containing w, or the closest point to w (so far) */ + else{ + BamgQuadtreeBox* b0=b; + + /*if the next box exists:*/ + if (b=b->b[k]){ + + /*Get size (hb) and coordinates of the current sub-box lowest left corner*/ + hb>>=1; + Icoor1 iii = ii[level]+I_IJ(k,hb); + Icoor1 jjj = jj[level]+J_IJ(k,hb); + + /*if the current point (iplus,jplus) is in b (modulo h), this box is good: + * it is holding vertices that are close to w */ + if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){ + level++; + pb[level]= b; + pi[level]= b->nbitems>0 ?(int) b->nbitems : 4 ; + ii[level]= iii; + jj[level]= jjj; + } + + //else go backwards + else{ + //shifted righ by one bit: hb=001000000 -> 01000000 + b=b0; + hb<<=1; + } + } + else{ + /*Current box is NULL, go to next subbox of b (k=k-1)*/ + b=b0; + } + } + } + + /*We have found a vertex, now, let's try the other boxes of the previous level + * in case there is a vertex closest to w that has not yet been tested*/ + hb <<= 1; + } while (level--); + + /*return nearest_v, nearest vertex*/ + return nearest_v; + + } + /*}}}1*/ + /*FUNCTION BamgQuadtree::NearestVertexWithNormal{{{1*/ + BamgVertex* BamgQuadtree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) { + /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/NearestVertexWithNormal)*/ + + BamgQuadtreeBox * pb[ MaxDepth ]; + int pi[ MaxDepth ]; + Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; + int l; // level + BamgQuadtreeBox * b; + long h =MaxISize,h0; + long hb=MaxISize; + Icoor1 i0=0,j0=0; + Icoor1 iplus( inbitems) + return vn; // empty tree + + while( (n0 = b->nbitems) < 0) + { + // search the non empty + // BamgQuadtreeBox containing the point (i,j) + register Icoor1 hb2 = hb >> 1 ; + register int k = IJ(iplus,jplus,hb2);// BamgQuadtreeBox number of size hb2 contening i;j + register BamgQuadtreeBox * b0= b->b[k]; + if ( ( b0 == 0) || (b0->nbitems == 0) ) + break; // null box or empty => break + b=b0; + i0 += I_IJ(k,hb2); // i orign of BamgQuadtreeBox + j0 += J_IJ(k,hb2); // j orign of BamgQuadtreeBox + hb = hb2; + } + + + if ( n0 > 0) + { + for(register int k=0;kv[k]->i; + // try if is in the right direction -- + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k];} + } + if (vn) return vn; + } + // general case ----- + // INITIALISATION OF THE HEAP + l =0; // level + pb[0]= b; + pi[0]=b->nbitems>0 ?(int) b->nbitems : 4 ; + ii[0]=i0; + jj[0]=j0; + h=hb; + do { // walk on the tree + b= pb[l]; + while (pi[l]--) // loop on 4 element of the box + { + int k = pi[l]; + + if (b->nbitems>0) // BamgVertex BamgQuadtreeBox none empty + { + I2 i2 = b->v[k]->i; + // if good direction when try -- + + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k]; + } + } + else // Pointer BamgQuadtreeBox + { + register BamgQuadtreeBox *b0=b; + if ((b=b->b[k])) + { + hb >>=1 ; // div by 2 + register Icoor1 iii = ii[l]+I_IJ(k,hb); + register Icoor1 jjj = jj[l]+J_IJ(k,hb); + + if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)) + { + pb[++l]= b; + pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; + ii[l]= iii; + jj[l]= jjj; + + } + else + b=b0, hb <<=1 ; + } + else + b=b0; + } + } + hb <<= 1; // mul by 2 + } while (l--); + + return vn; + } + /*}}}1*/ + /*FUNCTION BamgQuadtree::NewBamgQuadtreeBox {{{1*/ + BamgQuadtree::BamgQuadtreeBox* BamgQuadtree::NewBamgQuadtreeBox(void){ + + /*Output*/ + BamgQuadtreeBox* newbox=NULL; + + /*Create and initialize a new box*/ + newbox=new BamgQuadtreeBox; + newbox->nbitems=0; + newbox->b[0]=NULL; + newbox->b[1]=NULL; + newbox->b[2]=NULL; + newbox->b[3]=NULL; + + /*Add root to the container*/ + boxcontainer->AddObject(newbox); + + /*Increase counter*/ + NbBamgQuadtreeBox++; + + /*currentbox now points toward next quadtree box*/ + return newbox; + }/*}}}*/ + /*FUNCTION BamgQuadtree::ToClose {{{1*/ + BamgVertex* BamgQuadtree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){ + /*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.cpp/ToClose)*/ + + const Icoor1 i=v.i.x; + const Icoor1 j=v.i.y; + const R2 X(v.r); + const Metric Mx(v.m); + + BamgQuadtreeBox * pb[ MaxDepth ]; + int pi[ MaxDepth ]; + Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; + register int l=0; // level + register BamgQuadtreeBox * b; + Icoor1 h=MaxISize; + Icoor1 hb = MaxISize; + Icoor1 i0=0,j0=0; + + // BamgVertex *vn=0; + + if (!root->nbitems) + return 0; // empty tree + + // general case ----- + pb[0]=root; + pi[0]=root->nbitems>0 ?(int) root->nbitems : 4 ; + ii[0]=i0; + jj[0]=j0; + h=hb; + do { + b= pb[l]; + while (pi[l]--){ + register int k = pi[l]; + + if (b->nbitems>0){ // BamgVertex BamgQuadtreeBox none empty + I2 i2 = b->v[k]->i; + if ( ABS(i-i2.x) v[k]->r); + double dd; + if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ){ + return b->v[k]; + } + } + } + else{ // Pointer BamgQuadtreeBox + register BamgQuadtreeBox *b0=b; + if ((b=b->b[k])){ + hb >>=1 ; // div by 2 + register long iii = ii[l]+I_IJ(k,hb); + register long jjj = jj[l]+J_IJ(k,hb); + + if (INTER_SEG(iii,iii+hb,i-hx,i+hx) && INTER_SEG(jjj,jjj+hb,j-hy,j+hy)){ + pb[++l]= b; + pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; + ii[l]= iii; + jj[l]= jjj; + + } + else{ + b=b0; + hb <<=1 ; + } + } + else{ + b=b0; + } + } + } + hb <<= 1; // mul by 2 + } while (l--); + + return 0; + } + /*}}}1*/ +} Index: /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h (revision 12218) +++ /issm/trunk-jpl/src/c/objects/Bamg/BamgQuadtree.h (revision 12218) @@ -0,0 +1,60 @@ +/*Original code from Frederic Hecht (BAMG v1.01, BamgQuadtree.h)*/ +#ifndef _QUADTREE_H +#define _QUADTREE_H + +#include "./include.h" + +namespace bamg { + + const int MaxDepth = 30; + const long MaxISize = ( 1L << MaxDepth); // = 2^30 : 010000000000..000 (bitwise operation) + + class BamgVertex; + + class BamgQuadtree{ + + private: + + /*A quadtree box contains a maximum of 4 vertices. 4 other quadtree boxes are + * created if a fifth vertex is added to the same box. A Quadtree box is therefore + * composed of EITHER: + * - up to 4 vertices + * - 4 "sub" quadtree boxes*/ + class BamgQuadtreeBox: public Object{ + public: + int nbitems; // number of current vertices in the box + union{ + BamgQuadtreeBox* b[4]; + BamgVertex* v[4]; + }; + /*Object functions*/ + void Echo() {_error_("not implemented yet"); }; + void DeepEcho() {_error_("not implemented yet"); }; + int Id() {_error_("not implemented yet"); }; + int MyRank() {_error_("not implemented yet"); }; + int ObjectEnum() {_error_("not implemented yet"); }; + Object *copy() {_error_("not implemented yet"); }; + }; + + /*BamgQuadtree private Fields*/ + DataSet* boxcontainer; + + public: + + /*BamgQuadtree public Fields*/ + BamgQuadtreeBox* root; + long NbBamgQuadtreeBox; + long NbVertices; + + BamgQuadtree(); + BamgQuadtree(Mesh *t,long nbv=-1); + ~BamgQuadtree(); + + BamgVertex *NearestVertex(Icoor1 i,Icoor1 j); + BamgVertex *NearestVertexWithNormal(Icoor1 i,Icoor1 j); + BamgQuadtreeBox *NewBamgQuadtreeBox(void); + BamgVertex *ToClose(BamgVertex &,double ,Icoor1,Icoor1); + void Add(BamgVertex &w); + }; +} +#endif Index: /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp (revision 12217) +++ /issm/trunk-jpl/src/c/objects/Bamg/Geometry.cpp (revision 12218) @@ -505,5 +505,5 @@ float *eangle = new float[nbe]; double eps = 1e-20; - QuadTree quadtree; // build quadtree to find duplicates + BamgQuadtree quadtree; // build quadtree to find duplicates BamgVertex *v0 = vertices; GeomVertex *v0g = (GeomVertex*) (void*)v0; Index: /issm/trunk-jpl/src/c/objects/Bamg/Geometry.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/Geometry.h (revision 12217) +++ /issm/trunk-jpl/src/c/objects/Bamg/Geometry.h (revision 12218) @@ -12,5 +12,5 @@ class Triangle; - class QuadTree; + class BamgQuadtree; class GeomSubDomain; class Edge; @@ -20,17 +20,17 @@ public: - long NbRef; // counter of ref on the this class if 0 we can delete - long nbv; // number of vertices - long nbe; // number of edges - long nbsubdomains; - long nbcurves; + long NbRef; // counter of ref on the this class if 0 we can delete + long nbv; // number of vertices + long nbe; // number of edges + long nbsubdomains; + long nbcurves; GeomVertex *vertices; GeomEdge *edges; - QuadTree *quadtree; + BamgQuadtree *quadtree; GeomSubDomain *subdomains; - Curve *curves; - R2 pmin,pmax; // domain extrema coordinates - double coefIcoor; // coef to integer Icoor1; - double MaxCornerAngle; + Curve *curves; + R2 pmin,pmax; // domain extrema coordinates + double coefIcoor; // coef to integer Icoor1; + double MaxCornerAngle; //Constructor/Destructors @@ -44,5 +44,5 @@ GeomVertex &operator[](long i) { return vertices[i]; }; const GeomEdge &operator()(long i) const { return edges[i]; }; - GeomEdge &operator()(long i) { return edges[i]; }; + GeomEdge &operator()(long i) { return edges[i]; }; //Methods Index: /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp (revision 12217) +++ /issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp (revision 12218) @@ -2885,5 +2885,5 @@ //build quadtree - if (!quadtree) quadtree = new QuadTree(this,0); + if (!quadtree) quadtree = new BamgQuadtree(this,0); quadtree->Add(*v0); quadtree->Add(*v1); @@ -3107,10 +3107,10 @@ } /*}}}1*/ - /*FUNCTION Mesh::MakeQuadTree{{{1*/ - void Mesh::MakeQuadTree() { - /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/MakeQuadTree)*/ + /*FUNCTION Mesh::MakeBamgQuadtree{{{1*/ + void Mesh::MakeBamgQuadtree() { + /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/MakeBamgQuadtree)*/ long int verbose=0; - if ( !quadtree ) quadtree = new QuadTree(this); + if ( !quadtree ) quadtree = new BamgQuadtree(this); } @@ -3622,5 +3622,5 @@ if (!quadtree) delete quadtree; //ReInitialise; - quadtree = new QuadTree(this,0); + quadtree = new BamgQuadtree(this,0); quadtree->Add(*v0); quadtree->Add(*v1); @@ -3900,5 +3900,5 @@ if (quadtree){ delete quadtree; - quadtree = new QuadTree(this); + quadtree = new BamgQuadtree(this); } @@ -4116,5 +4116,5 @@ delete [] tstart; - if (quadtree) quadtree= new QuadTree(this); + if (quadtree) quadtree= new BamgQuadtree(this); } /*}}}1*/ Index: /issm/trunk-jpl/src/c/objects/Bamg/Mesh.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/Mesh.h (revision 12217) +++ /issm/trunk-jpl/src/c/objects/Bamg/Mesh.h (revision 12218) @@ -17,5 +17,5 @@ class Geometry; class CrackedEdge; - class QuadTree; + class BamgQuadtree; class SubDomain; @@ -29,5 +29,5 @@ Triangle *triangles; Edge *edges; - QuadTree *quadtree; + BamgQuadtree *quadtree; BamgVertex **orderedvertices; SubDomain *subdomains; @@ -94,5 +94,5 @@ void MakeQuadrangles(double costheta); int SplitElement(int choice); - void MakeQuadTree(); + void MakeBamgQuadtree(); void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); long InsertNewPoints(long nbvold,long & NbTSwap) ; Index: sm/trunk-jpl/src/c/objects/Bamg/QuadTree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/QuadTree.cpp (revision 12217) +++ (revision ) @@ -1,598 +1,0 @@ -#include -#include -#include - -#include "../objects.h" - -namespace bamg { - - /*MACROS {{{1*/ - /* - * - * J j - * ^ ^ - * | | +--------+--------+ - * | | | | | - * 1X | | | 2 | 3 | - * | | | | | - * | | +--------+--------+ - * | | | | | - * 0X | | | 0 | 1 | - * | | | | | - * | | +--------+--------+ - * | +-----------------------> i - * | - * |----------------------------> I - * X0 X1 - * - * box 0 -> I=0 J=0 IJ=00 = 0 - * box 1 -> I=1 J=0 IJ=01 = 1 - * box 2 -> I=0 J=1 IJ=10 = 2 - * box 3 -> I=1 J=1 IJ=11 = 3 - */ -#define INTER_SEG(a,b,x,y) (((y) > (a)) && ((x) <(b))) -#define ABS(i) ((i)<0 ?-(i) :(i)) -#define MAX1(i,j) ((i)>(j) ?(i) :(j)) -#define NORM(i1,j1,i2,j2) MAX1(ABS((i1)-(j1)),ABS((i2)-(j2))) - - //IJ(i,j,l) returns the box number of i and j with respect to l - //if !j&l and !i&l -> 0 (box zero: lower left ) - //if !j&l and i&l -> 1 (box one: lower right) - //if j&l and !i&l -> 2 (box two: upper left ) - //if j&l and i&l -> 3 (box three:upper right) -#define IJ(i,j,l) ((j&l) ? ((i&l) ? 3:2 ) :((i&l) ? 1:0 )) - - //I_IJ(k,l) returns l if first bit of k is 1, else 0 -#define I_IJ(k,l) ((k&1) ? l:0) - //J_IJ(k,l) returns l if second bit of k is 1, else 0 -#define J_IJ(k,l) ((k&2) ? l:0) - /*}}}*/ - /*DOCUMENTATION What is a QuadTree? {{{1 - * A Quadtree is a very simple way to group vertices according - * to their locations. A square that holds all the points of the mesh - * (or the geometry) is divided into 4 boxes. As soon as one box - * hold more than 4 vertices, it is divided into 4 new boxes, etc... - * There cannot be more than MAXDEEP (=30) subdivision. - * This process is like a Dichotomy in dimension 2 - * - * + - - - - - - - - + - - + - + - + - - - - + - * | | | | X | | - * + - + - + - * | | | | | | - * + - - + - + - + + - * | | | | | - * - * | | | | | - * + - - - - - - - - + - - + - - + - - - - + - * | | | | - * - * | | | | - * - * | | | | - * | | | | - * + - - - - - - - - + - - - - + - - - - + - * | | | - * - * | | | - * - * | | | - * - * | | | - * | | | - * | | | - * | | | - * | | | - * + - - - - - - - - + - - - - - - - - + - * - * The coordinate system used in a quadtree are integers to avoid - * round-off errors. The vertex in the lower left box has the coordinates - * (0 0) - * The upper right vertex has the follwing coordinates: - * 2^30 -1 2^30 -1 in decimal - * 0 1 1 1 .... 1 0 1 1 1 .... 1 in binary - * \-- 29 --/ \-- 29 --/ - * Using binaries is therefore very easy to locate a vertex in a box: - * we just need to look at the bits from the left to the right (See ::Add) - }}}1*/ - - /*Constructors/Destructors*/ - /*FUNCTION QuadTree::QuadTree(){{{1*/ - QuadTree::QuadTree(){ - - /*Number of boxes and vertices*/ - NbQuadTreeBox=0; - NbVertices=0; - - /*Create container*/ - boxcontainer=new DataSet(); - - /*Create Root, pointer toward the main box*/ - root=NewQuadTreeBox(); - - } - /*}}}1*/ - /*FUNCTION QuadTree::QuadTree(Mesh * t,long nbv){{{1*/ - QuadTree::QuadTree(Mesh * t,long nbv){ - - /*Number of boxes and vertices*/ - NbQuadTreeBox=0; - NbVertices=0; - - /*Create container*/ - boxcontainer=new DataSet(); - - /*Create Root, pointer toward the main box*/ - root=NewQuadTreeBox(); - - /*Check Sizes*/ - _assert_(MaxISize>MaxICoor); - - /*Add all vertices of the mesh*/ - if (nbv==-1) nbv=t->nbv; - for (int i=0;ivertices[i]); - - } - /*}}}1*/ - /*FUNCTION QuadTree::~QuadTree(){{{1*/ - QuadTree::~QuadTree() { - delete boxcontainer; - root=NULL; - } - /*}}}1*/ - - /*Methods*/ - /*FUNCTION QuadTree::Add{{{1*/ - void QuadTree::Add(BamgVertex &w){ - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/Add)*/ - QuadTreeBox** pb=NULL; - QuadTreeBox* b=NULL; - - /*Get integer coodinate of current point w*/ - register long i=w.i.x, j=w.i.y; - - /*Initialize level*/ - register long level=MaxISize; - - /*Get inital box (the largest)*/ - pb = &root; - - /*Find the smallest box where w is located*/ - while((b=*pb) && (b->nbitems<0)){ - - //shift b->nbitems by -1 - b->nbitems--; - - //shifted righ by one bit: level=00000010 -> 00000001 - level >>= 1; - - //Get next subbox according to the bit value (level) - pb = &b->b[IJ(i,j,level)]; - } - - /*OK, we have found b, a Subbox holding vertices (might be full) - check that the vertex is not already in the box*/ - if (b){ - if (b->nbitems > 3 && b->v[3] == &w) return; - if (b->nbitems > 2 && b->v[2] == &w) return; - if (b->nbitems > 1 && b->v[1] == &w) return; - if (b->nbitems > 0 && b->v[0] == &w) return; - } - - /*check that l is not 0 (this should not happen as MaxDepth = 30)*/ - _assert_(level>0); - - /*Now, try to add the vertex, if the subbox is full (nbitems=4), we have to divide it - in 4 new subboxes*/ - while ((b= *pb) && (b->nbitems == 4)){ // the QuadTreeBox is full - - /*Copy the 4 vertices in the current QuadTreebox*/ - BamgVertex* v4[4]; - v4[0]= b->v[0]; - v4[1]= b->v[1]; - v4[2]= b->v[2]; - v4[3]= b->v[3]; - - /*set nbitems as negative - * (box full -> holds 4 pointers toward subboxes and not 4 vertices)*/ - b->nbitems = -b->nbitems; - - /*Initialize the 4 pointers toward the 4 subboxes*/ - b->b[0]=b->b[1]=b->b[2]=b->b[3]=NULL; - - /*level = 0010000 -> 0001000*/ - level >>= 1; - - /*Put the four vertices in the new boxes*/ - for (int k=0;k<4;k++){ - - int ij; - /*bb is the new "sub"box of b where v4[k] is located*/ - QuadTreeBox *bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,level)]; - - // alloc the QuadTreeBox - if (!bb) bb=b->b[ij]=NewQuadTreeBox(); - - /*Copy the current vertex*/ - bb->v[bb->nbitems++] = v4[k]; - } - - /*Get the subbox where w (i,j) is located*/ - pb = &b->b[IJ(i,j,level)]; - } - - /*alloc the QuadTreeBox if necessary*/ - if (!(b=*pb)) b=*pb= NewQuadTreeBox(); - - /*Add w*/ - b->v[b->nbitems++]=&w; - - //Increase NbVertices by one (we have one new vertex) - NbVertices++; - } - /*}}}1*/ - /*FUNCTION QuadTree::NearestVertex{{{1*/ - BamgVertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/NearestVertex)*/ - - /*Intermediaries*/ - QuadTreeBox *pb[MaxDepth]; - int pi[MaxDepth]; - Icoor1 ii[MaxDepth]; - Icoor1 jj[MaxDepth]; - int level; - long n0; - QuadTreeBox *b; - long h0; - long h = MaxISize; - long hb= MaxISize; - Icoor1 i0=0,j0=0; - - /*initial output as NULL (no vertex found)*/ - BamgVertex* nearest_v=NULL; - - /*Project w coordinates (i,j) onto [0,MaxISize-1] x [0,MaxISize-1] -> (iplus,jplus)*/ - Icoor1 iplus( inbitems) return nearest_v; - - /*else, find the smallest non-empty QuadTreeBox containing the point (i,j)*/ - while((n0=b->nbitems)<0){ - - Icoor1 hb2 = hb >> 1; //size of the current box - int k = IJ(iplus,jplus,hb2); //box number (0,1,2 or 3) - QuadTreeBox *b0 = b->b[k]; //pointer toward current box - - /* break if NULL box or empty (Keep previous box b)*/ - if (( b0 == NULL) || (b0->nbitems == 0)) break; - - /*Get next Quadtree box*/ - b=b0; - i0 += I_IJ(k,hb2); // i orign of QuadTreeBox (macro) - j0 += J_IJ(k,hb2); // j orign of QuadTreeBox - hb = hb2; // size of the box (in Int) - } - - /*The box b, is the smallest box containing the point (i,j) and - * has the following properties: - * - n0: number of items (>0 if vertices, else boxes) - * - hb: box size (int) - * - i0: x coordinate of the lower left corner - * - j0: y coordinate of the lower left corner*/ - - /* if the current subbox is holding vertices, we are almost done*/ - if (n0>0){ - //loop over the vertices of the box and find the closest vertex - for(int k=0;kv[k]->i; - - /*Compute norm with w*/ - h0=NORM(iplus,i2.x,jplus,i2.y); - - /*is it smaller than previous value*/ - if (h0v[k]; - } - } - /*return closest vertex*/ - return nearest_v; - } - - /* general case: the current box is empty, we have to go backwards - and find the closest not-empty box and find the closest vertex*/ - - /*Initialize search variables*/ - pb[0]=b; //pointer toward the box b - pi[0]=b->nbitems>0?(int)b->nbitems:4;//number of boxes in b - ii[0]=i0; //i coordinate of the box lowest left corner - jj[0]=j0; //j coordinate of the box lowest left corner - - /*initialize h: smallest box size, containing a vertex close to w*/ - h=hb; - - /*Main loop*/ - level=0; - do { - - /*get current box*/ - b= pb[level]; - - /*Loop over the items in current box (if not empty!)*/ - while (pi[level]){ - - /*We are looping now over the items of b. k is the current index (in [0 3])*/ - pi[level]--; - int k=pi[level]; - - /*if the current subbox is holding vertices (b->nbitems<0 is subboxes)*/ - if (b->nbitems>0){ - I2 i2=b->v[k]->i; - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0v[k]; - } - } - /*else: current box b is pointing toward 4 boxes - * test sub-box k and go deeper into the tree if it is non empty - * and contains the point w modulo a size h that is either the size of the smallest - * non empty box containing w, or the closest point to w (so far) */ - else{ - QuadTreeBox* b0=b; - - /*if the next box exists:*/ - if (b=b->b[k]){ - - /*Get size (hb) and coordinates of the current sub-box lowest left corner*/ - hb>>=1; - Icoor1 iii = ii[level]+I_IJ(k,hb); - Icoor1 jjj = jj[level]+J_IJ(k,hb); - - /*if the current point (iplus,jplus) is in b (modulo h), this box is good: - * it is holding vertices that are close to w */ - if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){ - level++; - pb[level]= b; - pi[level]= b->nbitems>0 ?(int) b->nbitems : 4 ; - ii[level]= iii; - jj[level]= jjj; - } - - //else go backwards - else{ - //shifted righ by one bit: hb=001000000 -> 01000000 - b=b0; - hb<<=1; - } - } - else{ - /*Current box is NULL, go to next subbox of b (k=k-1)*/ - b=b0; - } - } - } - - /*We have found a vertex, now, let's try the other boxes of the previous level - * in case there is a vertex closest to w that has not yet been tested*/ - hb <<= 1; - } while (level--); - - /*return nearest_v, nearest vertex*/ - return nearest_v; - - } - /*}}}1*/ - /*FUNCTION QuadTree::NearestVertexWithNormal{{{1*/ - BamgVertex* QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) { - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/NearestVertexWithNormal)*/ - - QuadTreeBox * pb[ MaxDepth ]; - int pi[ MaxDepth ]; - Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; - int l; // level - QuadTreeBox * b; - long h =MaxISize,h0; - long hb=MaxISize; - Icoor1 i0=0,j0=0; - Icoor1 iplus( inbitems) - return vn; // empty tree - - while( (n0 = b->nbitems) < 0) - { - // search the non empty - // QuadTreeBox containing the point (i,j) - register Icoor1 hb2 = hb >> 1 ; - register int k = IJ(iplus,jplus,hb2);// QuadTreeBox number of size hb2 contening i;j - register QuadTreeBox * b0= b->b[k]; - if ( ( b0 == 0) || (b0->nbitems == 0) ) - break; // null box or empty => break - b=b0; - i0 += I_IJ(k,hb2); // i orign of QuadTreeBox - j0 += J_IJ(k,hb2); // j orign of QuadTreeBox - hb = hb2; - } - - - if ( n0 > 0) - { - for(register int k=0;kv[k]->i; - // try if is in the right direction -- - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k];} - } - if (vn) return vn; - } - // general case ----- - // INITIALISATION OF THE HEAP - l =0; // level - pb[0]= b; - pi[0]=b->nbitems>0 ?(int) b->nbitems : 4 ; - ii[0]=i0; - jj[0]=j0; - h=hb; - do { // walk on the tree - b= pb[l]; - while (pi[l]--) // loop on 4 element of the box - { - int k = pi[l]; - - if (b->nbitems>0) // BamgVertex QuadTreeBox none empty - { - I2 i2 = b->v[k]->i; - // if good direction when try -- - - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k]; - } - } - else // Pointer QuadTreeBox - { - register QuadTreeBox *b0=b; - if ((b=b->b[k])) - { - hb >>=1 ; // div by 2 - register Icoor1 iii = ii[l]+I_IJ(k,hb); - register Icoor1 jjj = jj[l]+J_IJ(k,hb); - - if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)) - { - pb[++l]= b; - pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; - ii[l]= iii; - jj[l]= jjj; - - } - else - b=b0, hb <<=1 ; - } - else - b=b0; - } - } - hb <<= 1; // mul by 2 - } while (l--); - - return vn; - } - /*}}}1*/ - /*FUNCTION QuadTree::NewQuadTreeBox {{{1*/ - QuadTree::QuadTreeBox* QuadTree::NewQuadTreeBox(void){ - - /*Output*/ - QuadTreeBox* newbox=NULL; - - /*Create and initialize a new box*/ - newbox=new QuadTreeBox; - newbox->nbitems=0; - newbox->b[0]=NULL; - newbox->b[1]=NULL; - newbox->b[2]=NULL; - newbox->b[3]=NULL; - - /*Add root to the container*/ - boxcontainer->AddObject(newbox); - - /*Increase counter*/ - NbQuadTreeBox++; - - /*currentbox now points toward next quadtree box*/ - return newbox; - }/*}}}*/ - /*FUNCTION QuadTree::ToClose {{{1*/ - BamgVertex* QuadTree::ToClose(BamgVertex & v,double seuil,Icoor1 hx,Icoor1 hy){ - /*Original code from Frederic Hecht (BAMG v1.01, QuadTree.cpp/ToClose)*/ - - const Icoor1 i=v.i.x; - const Icoor1 j=v.i.y; - const R2 X(v.r); - const Metric Mx(v.m); - - QuadTreeBox * pb[ MaxDepth ]; - int pi[ MaxDepth ]; - Icoor1 ii[ MaxDepth ], jj [ MaxDepth]; - register int l=0; // level - register QuadTreeBox * b; - Icoor1 h=MaxISize; - Icoor1 hb = MaxISize; - Icoor1 i0=0,j0=0; - - // BamgVertex *vn=0; - - if (!root->nbitems) - return 0; // empty tree - - // general case ----- - pb[0]=root; - pi[0]=root->nbitems>0 ?(int) root->nbitems : 4 ; - ii[0]=i0; - jj[0]=j0; - h=hb; - do { - b= pb[l]; - while (pi[l]--){ - register int k = pi[l]; - - if (b->nbitems>0){ // BamgVertex QuadTreeBox none empty - I2 i2 = b->v[k]->i; - if ( ABS(i-i2.x) v[k]->r); - double dd; - if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ){ - return b->v[k]; - } - } - } - else{ // Pointer QuadTreeBox - register QuadTreeBox *b0=b; - if ((b=b->b[k])){ - hb >>=1 ; // div by 2 - register long iii = ii[l]+I_IJ(k,hb); - register long jjj = jj[l]+J_IJ(k,hb); - - if (INTER_SEG(iii,iii+hb,i-hx,i+hx) && INTER_SEG(jjj,jjj+hb,j-hy,j+hy)){ - pb[++l]= b; - pi[l]= b->nbitems>0 ?(int) b->nbitems : 4 ; - ii[l]= iii; - jj[l]= jjj; - - } - else{ - b=b0; - hb <<=1 ; - } - } - else{ - b=b0; - } - } - } - hb <<= 1; // mul by 2 - } while (l--); - - return 0; - } - /*}}}1*/ -} Index: sm/trunk-jpl/src/c/objects/Bamg/QuadTree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Bamg/QuadTree.h (revision 12217) +++ (revision ) @@ -1,61 +1,0 @@ -/*Original code from Frederic Hecht (BAMG v1.01, QuadTree.h)*/ -#ifndef _QUADTREE_H -#define _QUADTREE_H - -#include "./include.h" - -namespace bamg { - - const int MaxDepth = 30; - const long MaxISize = ( 1L << MaxDepth); // = 2^30 : 010000000000..000 (bitwise operation) - - class BamgVertex; - - class QuadTree{ - - private: - - /*A quadtree box contains a maximum of 4 vertices. 4 other quadtree boxes are - * created if a fifth vertex is added to the same box. A Quadtree box is therefore - * composed of EITHER: - * - up to 4 vertices - * - 4 "sub" quadtree boxes*/ - class QuadTreeBox: public Object{ - public: - int nbitems; // number of current vertices in the box - union{ - QuadTreeBox* b[4]; - BamgVertex* v[4]; - }; - /*Object functions*/ - void Echo(){_error_("not implemented yet");}; - void DeepEcho(){_error_("not implemented yet");}; - int Id(){_error_("not implemented yet");}; - int MyRank(){_error_("not implemented yet");}; - int ObjectEnum(){_error_("not implemented yet");}; - Object* copy(){_error_("not implemented yet");}; - }; - - /*QuadTree private Fields*/ - DataSet* boxcontainer; - - public: - - /*QuadTree public Fields*/ - QuadTreeBox* root; - long NbQuadTreeBox; - long NbVertices; - - QuadTree(); - QuadTree(Mesh *t,long nbv=-1); - ~QuadTree(); - - BamgVertex* NearestVertex(Icoor1 i,Icoor1 j); - BamgVertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j); - QuadTreeBox* NewQuadTreeBox(void); - BamgVertex* ToClose(BamgVertex & ,double ,Icoor1,Icoor1); - void Add( BamgVertex & w); - - }; -} -#endif Index: /issm/trunk-jpl/src/c/objects/objects.h =================================================================== --- /issm/trunk-jpl/src/c/objects/objects.h (revision 12217) +++ /issm/trunk-jpl/src/c/objects/objects.h (revision 12218) @@ -170,5 +170,5 @@ #include "./Bamg/Mesh.h" #include "./Bamg/Geometry.h" -#include "./Bamg/QuadTree.h" +#include "./Bamg/BamgQuadtree.h" #include "./Bamg/SetOfE4.h"