Index: sm/trunk/src/c/Bamgx/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.cpp (revision 2840) +++ (revision ) @@ -1,456 +1,0 @@ -// -*- Mode : c++ -*- -// -// SUMMARY : -// USAGE : -// ORG : -// AUTHOR : Frederic Hecht -// E-MAIL : hecht@ann.jussieu.fr -// - -/* - - This file is part of Freefem++ - - Freefem++ is free software; you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation; either version 2.1 of the License, or - (at your option) any later version. - - Freefem++ is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with Freefem++; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - */ - -#include -//#include -#include -#include - -#include "Mesh2.h" -#include "QuadTree.h" - -namespace bamg { - -#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))) - -#define IJ(i,j,l) ( ( j & l) ? (( i & l) ? 3 : 2 ) :( ( i & l)? 1 : 0 )) -#define I_IJ(k,l) (( k&1) ? l : 0) -#define J_IJ(k,l) (( k&2) ? l : 0) - - - -Vertex * QuadTree::NearestVertex(Icoor1 i,Icoor1 j) -{ - QuadTreeBox * pb[ MaxDeep ]; - int pi[ MaxDeep ]; - Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; - register int l=0; // level - register QuadTreeBox * b; - IntQuad h=MaxISize,h0; - IntQuad hb = MaxISize; - Icoor1 i0=0,j0=0; - Icoor1 iplus( in) - return vn; // empty tree - - while( (n0 = b->n) < 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->n == 0) ) - break; // null box or empty => break - NbQuadTreeBoxSearch++; - 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; - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k];} - NbVerticesSearch++; - } - return vn; - } - // general case ----- - pb[0]= b; - pi[0]=b->n>0 ?(int) b->n : 4 ; - ii[0]=i0; - jj[0]=j0; - h=hb; - do { - b= pb[l]; - while (pi[l]--) - { - register int k = pi[l]; - - if (b->n>0) // Vertex QuadTreeBox none empty - { - NbVerticesSearch++; - I2 i2 = b->v[k]->i; - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k]; - } - } - else // Pointer QuadTreeBox - { - register QuadTreeBox *b0=b; - NbQuadTreeBoxSearch++; - 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->n>0 ?(int) b->n : 4 ; - ii[l]= iii; - jj[l]= jjj; - - } - else - b=b0, hb <<=1 ; - } - else - b=b0; - } - } - hb <<= 1; // mul by 2 - } while (l--); - - return vn; -} - - - -Vertex * QuadTree::ToClose(Vertex & v,Real8 seuil,Icoor1 hx,Icoor1 hy) -{ - const Icoor1 i=v.i.x; - const Icoor1 j=v.i.y; - const R2 X(v.r); - const Metric Mx(v.m); - - QuadTreeBox * pb[ MaxDeep ]; - int pi[ MaxDeep ]; - Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; - register int l=0; // level - register QuadTreeBox * b; - Icoor1 h=MaxISize; - Icoor1 hb = MaxISize; - Icoor1 i0=0,j0=0; - - // Vertex *vn=0; - - if (!root->n) - return 0; // empty tree - - // general case ----- - pb[0]=root; - pi[0]=root->n>0 ?(int) root->n : 4 ; - ii[0]=i0; - jj[0]=j0; - h=hb; - do { - b= pb[l]; - while (pi[l]--) - { - register int k = pi[l]; - - if (b->n>0) // Vertex QuadTreeBox none empty - { - NbVerticesSearch++; - I2 i2 = b->v[k]->i; - if ( ABS(i-i2.x) v[k]->r); - Real8 dd; - // old code if( Mx(XY) + b->v[k]->m(XY) < seuil ) - if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ) - { - // cout << CurrentTh->Number(v) << "is To Close " - // << CurrentTh->Number( b->v[k]) << " l=" <v[k]; - } - } - } - else // Pointer QuadTreeBox - { - register QuadTreeBox *b0=b; - NbQuadTreeBoxSearch++; - 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->n>0 ?(int) b->n : 4 ; - ii[l]= iii; - jj[l]= jjj; - - } - else - b=b0, hb <<=1 ; - } - else - b=b0; - } - } - hb <<= 1; // mul by 2 - } while (l--); - - return 0; -} - - -void QuadTree::Add( Vertex & w) -{ - QuadTreeBox ** pb , *b; - register long i=w.i.x, j=w.i.y,l=MaxISize; - pb = &root; - // cout << pb << " " << &root << endl; - while( (b=*pb) && (b->n<0)) - { - b->n--; - l >>= 1; - pb = &b->b[IJ(i,j,l)]; - } - if (b) { - if (b->n > 3 && b->v[3] == &w) return; - if (b->n > 2 && b->v[2] == &w) return; - if (b->n > 1 && b->v[1] == &w) return; - if (b->n > 0 && b->v[0] == &w) return; - } - assert(l); - while ((b= *pb) && (b->n == 4)) // the QuadTreeBox is full - { - Vertex *v4[4]; // copy of the QuadTreeBox vertices - - v4[0]= b->v[0]; - v4[1]= b->v[1]; - v4[2]= b->v[2]; - v4[3]= b->v[3]; - b->n = -b->n; // mark is pointer QuadTreeBox - b->b[0]=b->b[1]=b->b[2]=b->b[3]=0; // set empty QuadTreeBox ptr - l >>= 1; // div the size by 2 - for (register int k=0;k<4;k++) // for the 4 vertices find the sub QuadTreeBox ij - { - register int ij; - register QuadTreeBox * bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,l)]; - if (!bb) - bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox - // cout << bb << " " << k << " " << ij << endl; - bb->v[bb->n++] = v4[k]; - } - pb = &b->b[IJ(i,j,l)]; - } - if (!(b = *pb)) - b=*pb= NewQuadTreeBox(); // alloc the QuadTreeBox - // cout << b << " " << b->n << endl; - b->v[b->n++]=&w; // we add the vertex - NbVertices++; -} - -QuadTree::QuadTree(Triangles * t,long nbv) : - lenStorageQuadTreeBox(t->nbvx/8+10), - th(t), - NbQuadTreeBox(0), - NbVertices(0), - NbQuadTreeBoxSearch(0), - NbVerticesSearch(0) -{ - if (nbv == -1) nbv = t->nbv; - sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); - root=NewQuadTreeBox(); - assert( MaxISize > MaxICoor); - for (Int4 i=0;ivertices[i]); -} - -QuadTree::QuadTree() : - lenStorageQuadTreeBox(100), - th(0), - NbQuadTreeBox(0), - NbVertices(0), - NbQuadTreeBoxSearch(0), - NbVerticesSearch(0) -{ - sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); - root=NewQuadTreeBox(); -} -QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) -{ - len = ll; - n = nn; - b = new QuadTreeBox[ll]; - for (int i = 0; i n) - return vn; // empty tree - - while( (n0 = b->n) < 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->n == 0) ) - break; // null box or empty => break - NbQuadTreeBoxSearch++; - 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 sens -- - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k];} - NbVerticesSearch++; - } - if (vn) return vn; - } - // general case ----- - // INITIALISATION OF THE HEAP - l =0; // level - pb[0]= b; - pi[0]=b->n>0 ?(int) b->n : 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->n>0) // Vertex QuadTreeBox none empty - { - NbVerticesSearch++; - I2 i2 = b->v[k]->i; - // if good sens when try -- - - h0 = NORM(iplus,i2.x,jplus,i2.y); - if (h0 v[k]; - } - } - else // Pointer QuadTreeBox - { - register QuadTreeBox *b0=b; - NbQuadTreeBoxSearch++; - 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->n>0 ?(int) b->n : 4 ; - ii[l]= iii; - jj[l]= jjj; - - } - else - b=b0, hb <<=1 ; - } - else - b=b0; - } - } - hb <<= 1; // mul by 2 - } while (l--); - - return vn; -} - - -} // end of namespace bamg - Index: /issm/trunk/src/c/Bamgx/QuadTree.h =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2840) +++ /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2841) @@ -1,30 +1,2 @@ -// -*- Mode : c++ -*- -// -// SUMMARY : -// USAGE : -// ORG : -// AUTHOR : Frederic Hecht -// E-MAIL : hecht@ann.jussieu.fr -// - -/* - - This file is part of Freefem++ - - Freefem++ is free software; you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation; either version 2.1 of the License, or - (at your option) any later version. - - Freefem++ is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with Freefem++; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - */ - namespace bamg { @@ -106,3 +78,2 @@ }; } -//#undef IJ Index: /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 2840) +++ /issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp (revision 2841) @@ -15,5 +15,5 @@ /*Constructor/Destructor*/ - /*Others*/ + /*Methods*/ /*FUNCTION GeometricalEdge::R1tg{{{1*/ Real8 GeometricalEdge::R1tg(Real8 theta,R2 & t) const // 1/R of radius of cuvature Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2840) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2841) @@ -17,4 +17,51 @@ /*Constructors/Destructors*/ + /*FUNCTION Geometry::Geometry(const Geometry & Gh){{{1*/ + Geometry::Geometry(const Geometry & Gh) { + Int4 i; + *this = Gh; + NbRef =0; + quadtree=0; + name = new char[strlen(Gh.name)+4]; + strcpy(name,"cp:"); + strcat(name,Gh.name); + vertices = nbv ? new GeometricalVertex[nbv] : NULL; + triangles = nbt ? new Triangle[nbt]:NULL; + edges = nbe ? new GeometricalEdge[nbe]:NULL; + curves= NbOfCurves ? new Curve[NbOfCurves]:NULL; + subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL; + for (i=0;i9) + cout << "DELETE ~Geometry "<< this << endl; + if(vertices) delete [] vertices;vertices=0; + if(edges) delete [] edges;edges=0; + // if(edgescomponante) delete [] edgescomponante; edgescomponante=0; + if(triangles) delete [] triangles;triangles=0; + if(quadtree) delete quadtree;quadtree=0; + if(curves) delete []curves;curves=0;NbOfCurves=0; + if(name) delete [] name;name=0; + if(subdomains) delete [] subdomains;subdomains=0; + // if(ordre) delete [] ordre; + EmptyGeometry(); + } + /*}}}1*/ /*IO*/ @@ -433,599 +480,550 @@ /*}}}1*/ - /*Other*/ - -void Geometry::EmptyGeometry() // empty geometry - { - OnDisk=0; - NbRef=0; - name =0; - quadtree=0; - curves=0; - // edgescomponante=0; - triangles=0; - edges=0; - vertices=0; - NbSubDomains=0; - // nbtf=0; -// BeginOfCurve=0; - nbiv=nbv=nbvx=0; - nbe=nbt=nbtx=0; - NbOfCurves=0; -// BeginOfCurve=0; - subdomains=0; - MaximalAngleOfCorner = 10*Pi/180; - } - - - -Geometry::Geometry(const Geometry & Gh) - { Int4 i; - *this = Gh; - NbRef =0; - quadtree=0; - name = new char[strlen(Gh.name)+4]; - strcpy(name,"cp:"); - strcat(name,Gh.name); - vertices = nbv ? new GeometricalVertex[nbv] : NULL; - triangles = nbt ? new Triangle[nbt]:NULL; - edges = nbe ? new GeometricalEdge[nbe]:NULL; - curves= NbOfCurves ? new Curve[NbOfCurves]:NULL; - subdomains = NbSubDomains ? new GeometricalSubDomain[NbSubDomains]:NULL; - for (i=0;iAdj[0]; - else if ( (AB,BP) > 0) - on = on->Adj[1]; - else - return on; - } - return on; -} -GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const - { - Real8 save_s=s; - int NbTry=0; + /*Methods*/ + /*FUNCTION Geometry::EmptyGeometry(){{{1*/ + void Geometry::EmptyGeometry() { + OnDisk=0; + NbRef=0; + name =0; + quadtree=0; + curves=0; + // edgescomponante=0; + triangles=0; + edges=0; + vertices=0; + NbSubDomains=0; + // nbtf=0; + // BeginOfCurve=0; + nbiv=nbv=nbvx=0; + nbe=nbt=nbtx=0; + NbOfCurves=0; + // BeginOfCurve=0; + subdomains=0; + MaximalAngleOfCorner = 10*Pi/180; + } + /*}}}1*/ + /*FUNCTION Geometry::Contening{{{1*/ + GeometricalEdge* Geometry::Contening(const R2 P, GeometricalEdge * start) const { + GeometricalEdge* on =start,* pon=0; + // walk with the cos on geometry + // cout << P ; + int k=0; + while(pon != on) + { + pon = on; + assert(k++<100); + R2 A= (*on)[0]; + R2 B= (*on)[1]; + R2 AB = B-A; + R2 AP = P-A; + R2 BP = P-B; + // cout << ":: " << on - edges << " " << AB*AP << " " << AB*BP << " " << A << B << endl; + if ( (AB,AP) < 0) + on = on->Adj[0]; + else if ( (AB,BP) > 0) + on = on->Adj[1]; + else + return on; + } + return on; + } + /*}}}1*/ + /*FUNCTION Geometry::Geometry::ProjectOnCurve {{{1*/ + GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const { + Real8 save_s=s; + int NbTry=0; retry: - s=save_s; - GeometricalEdge * on = e.on; - assert(on); - assert( e[0].on && e[1].on); - const Vertex &v0=e[0],&v1=e[1]; - V.m = Metric(1.0-s, v0,s, v1); + s=save_s; + GeometricalEdge * on = e.on; + assert(on); + assert( e[0].on && e[1].on); + const Vertex &v0=e[0],&v1=e[1]; + V.m = Metric(1.0-s, v0,s, v1); #define MXE__LINE __LINE__+1 - const int mxe =100; - GeometricalEdge *ge[mxe+1]; - int sensge[mxe+1]; - Real8 lge[mxe+1]; - int bge=mxe/2,tge=bge; - ge[bge] = e.on; - sensge[bge]=1; - - R2 V0 = v0,V1=v1,V01=V1-V0; - VertexOnGeom vg0= *v0.on, vg1=*v1.on; - if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl; - - // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; - GeometricalEdge * eg0=on, *eg1=on; - R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; - if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) < " << v1.r << endl - << "sg 0 = " << vg0 - << " on = " << Number(on) << ":" << Ag << Bg << "; " - << " sg 1= " << vg1 - << "--------------------------------------------" << endl; - while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0) - { - if (bge<=0) { - // int kkk; - // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ; - if(NbTry) - { - cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; - cerr << " The mesh of the Geometry is to fine: "; - cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; - cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; - cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; - cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; - MeshError(222); + const int mxe =100; + GeometricalEdge *ge[mxe+1]; + int sensge[mxe+1]; + Real8 lge[mxe+1]; + int bge=mxe/2,tge=bge; + ge[bge] = e.on; + sensge[bge]=1; + + R2 V0 = v0,V1=v1,V01=V1-V0; + VertexOnGeom vg0= *v0.on, vg1=*v1.on; + if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl; + + // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; + GeometricalEdge * eg0=on, *eg1=on; + R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; + if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) < " << v1.r << endl + << "sg 0 = " << vg0 + << " on = " << Number(on) << ":" << Ag << Bg << "; " + << " sg 1= " << vg1 + << "--------------------------------------------" << endl; + while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0) + { + if (bge<=0) { + // int kkk; + // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ; + if(NbTry) + { + cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; + cerr << " The mesh of the Geometry is to fine: "; + cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; + cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; + cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; + cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; + MeshError(222); + } + NbTry++; + goto retry;} + GeometricalEdge* tmpge = eg0; + if(NbTry) + cout << "bug: --Edge @" << Number(tmpge) << " = "<< Number(eg0) << ":" <Adj[0]) << "," << + Number(eg0->Adj[1]) <<"," ; + ge[--bge] =eg0 = eg0->Adj[sens0]; + assert(bge>=0 && bge <= mxe); + sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); + if(NbTry) + cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S " + << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << "," + << Number(eg0->Adj[1]) <<"," << endl + <Adj[0]) << "," << - Number(eg0->Adj[1]) <<"," ; - ge[--bge] =eg0 = eg0->Adj[sens0]; - assert(bge>=0 && bge <= mxe); - sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); - if(NbTry) - cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S " - << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << "," - << Number(eg0->Adj[1]) <<"," << endl - <=mxe ) { - cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; - NbTry++; - if (NbTry<2) goto retry; - cerr << " The mesh of the Geometry is to fine:" ; - cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; - cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; - cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; - cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; - MeshError(223); + if(NbTry) cout << Number((GeometricalEdge*) vg1) << " " << Number((GeometricalVertex*) vg1) << endl; + while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) + { + if(tge>=mxe ) { + cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; + NbTry++; + if (NbTry<2) goto retry; + cerr << " The mesh of the Geometry is to fine:" ; + cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; + cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; + cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; + cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; + MeshError(223); + } + + GeometricalEdge* tmpge = eg1; + if(NbTry) + cout << "++Edge @" << tmpge << " = " << Number(eg1) <<"%" << Number(eg1->Adj[0]) << "," + << Number(eg1->Adj[1]) <<"," ; + ge[++tge] =eg1 = eg1->Adj[sens1]; + sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; + assert(tge>=0 && tge <= mxe); + if(NbTry) + cout << " Edge " << Number(eg1) << " " << sens1 << " S " + <Adj[0]) << "," << Number(eg1->Adj[1]) <<"," + <=0 && i <= mxe); + BB = (*ge[i])[sensge[i]]; + lge[i]=ll += Norme2(AA-BB); + // cout << " ll " << i << BB << ll << " " <Adj[0]) << "," - << Number(eg1->Adj[1]) <<"," ; - ge[++tge] =eg1 = eg1->Adj[sens1]; - sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; - assert(tge>=0 && tge <= mxe); - if(NbTry) - cout << " Edge " << Number(eg1) << " " << sens1 << " S " - <Adj[0]) << "," << Number(eg1->Adj[1]) <<"," - <=0 && i <= mxe); - BB = (*ge[i])[sensge[i]]; - lge[i]=ll += Norme2(AA-BB); - // cout << " ll " << i << BB << ll << " " <r - vertices[i]) < eps ) - { // link v & vertices[i] - // vieille ruse pour recuperer j - GeometricalVertex * vg = (GeometricalVertex *) (void *) v; - int j = vg-v0g; - assert( v == & (Vertex &) vertices[j]); - vertices[i].link = vertices + j; - k++; - } - else quadtree.Add(vertices[i]); - } - if (k) { - cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl; - //if (verbosity>10) - { - cout << " The duplicate vertex " << endl; - for (i=0;i20) + cout << "Geometry::AfterRead()" << nbv << " " << nbe << endl; + Int4 i,k=0; ; + int jj; // jj in [0,1] + Int4 * hv = new Int4 [ nbv]; + Int4 * ev = new Int4 [ 2 * nbe ]; + float * eangle = new float[ nbe ]; + { + double eps = 1e-20; + QuadTree quadtree; // to find same vertices + Vertex * v0 = vertices; + GeometricalVertex * v0g = (GeometricalVertex *) (void *) v0; + int k=0; + for (i=0;ir - vertices[i]) < eps ) + { // link v & vertices[i] + // vieille ruse pour recuperer j + GeometricalVertex * vg = (GeometricalVertex *) (void *) v; + int j = vg-v0g; + assert( v == & (Vertex &) vertices[j]); + vertices[i].link = vertices + j; + k++; + } + else quadtree.Add(vertices[i]); + } + if (k) { + cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl; + //if (verbosity>10) + { + cout << " The duplicate vertex " << endl; + for (i=0;i7) + for (i=0;ir - edges[i].v[0]->r; + Real8 lv10 = Norme2(v10); + if(lv10 == 0) { + cerr << "The length of " <9) + cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<=0) + { + ord++; + register Int4 i1= n /2; + register Int4 j1 = n % 2; + register Int4 *pn = ev + n; + float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; + n = *pn; + if (angleold > angle) // exch to have : po -> pn -> p + exch=1,*pn = *po,*po=*p,*p=n,po = pn; + else // to have : po -> p -> pn + angleold = angle, po = p,p = pn; + } + } // end while (exch) + + if (ord >= 1 ) + { /* + Int4 n = hv[i]; + while ( n >=0) + { Int4 i1 = n/2,j1 = n%2; + //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi; + n = ev[n]; + } + */ + } + if(ord == 2) { // angulare test to find a corner + Int4 n1 = hv[i]; + Int4 n2 = ev[n1]; + Int4 i1 = n1 /2, i2 = n2/2; // edge number + Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge + float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; + float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; + float da12 = Abs(angle2-angle1); + if(verbosity>9) + cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi + << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl; + + if (( da12 >= MaximalAngleOfCorner ) + && (da12 <= 2*Pi -MaximalAngleOfCorner)) { + vertices[i].SetCorner() ; + if(verbosity>7) + cout << " The vertex " << i << " is a corner (angle) " + << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;} + // if the ref a changing then is SetRequired(); + + if (edges[i1].flag != edges[i2].flag || edges[i1].Required()) + { + vertices[i].SetRequired(); + if(verbosity>7) + cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;} + + if (edges[i1].ref != edges[i2].ref) { + vertices[i].SetRequired(); + if(verbosity>7) + cout << " The vertex " << i << " is Required ref" << endl;} + } ; + + if(ord != 2) { + vertices[i].SetCorner(); + if(verbosity>7) + cout << " the vertex " << i << " is a corner ordre = " << ord << endl; + } + // close the liste around the vertex + { Int4 no=-1, ne = hv[i]; + while ( ne >=0) + ne = ev[no=ne]; + if(no>=0) + ev[no] = hv[i]; + } // now the list around the vertex is circular + + } // end for (i=0;i " + << Number(edges[i1].v[j1]) <10) + cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl; + } + + // generation of all the tangente + for (i=0;ir -edges[i].v[0]->r; + Real8 lAB = Norme2(AB); // length of current edge AB + Real8 ltg2[2]; + ltg2[0]=0;ltg2[1]=0; + for (jj=0;jj<2;jj++) { + R2 tg = edges[i].tg[jj]; + Real8 ltg = Norme2(tg); // length of tg + if(ltg == 0) {// no tg + if( ! edges[i].v[jj]->Corner()) { // not a Corner + tg = edges[i].v[1-jj]->r + - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r; + ltg = Norme2(tg); + tg = tg *(lAB/ltg),ltg=lAB; + /* + if(edges[i].ref >=4) + cout << " tg " << tg.x << " "<< tg.y << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = " + << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl; + */ + } + + //else ;// a Corner with no tangent => nothing to do + } // a tg + else + tg = tg *(lAB/ltg),ltg=lAB; + ltg2[jj] = ltg; + if ( (tg,AB) < 0) + tg = -tg; + //if(edges[i].ref >=4) cout << " tg = " << tg << endl; + edges[i].tg[jj] = tg; + } // for (jj=0;jj<2;jj++) + + if (ltg2[0]!=0) edges[i].SetTgA(); + if (ltg2[1]!=0) edges[i].SetTgB(); + } // for (i=0;i7) + for (i=0;iSetMark(); + nbgem++; + e->CurveNumber=NbOfCurves; + if(curves) { + curves[NbOfCurves].ee=e; + curves[NbOfCurves].ke=k1; + } + + GeometricalVertex *b=(*e)(k1); + if (a == b || b->Required() ) break; + k0 = e->SensAdj[k1];// vertex in next edge + e = e->Adj[k1]; // next edge + + }// for(;;) + if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl; + NbOfCurves++; + if(level) { + if(verbosity>4) + cout << " Warning: Curve "<< NbOfCurves << " without required vertex " + << "so the vertex " << Number(a) << " become required " <SetRequired(); + } + + }} + assert(nbgem && nbe); + + if(step==0) { + curves = new Curve[NbOfCurves]; + } + } + for(int i=0;ilink; + //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; + curves[i].master=true; + if(be->Equi() || be->ReverseEqui() ) + { + assert(eqbe); + int nc = eqbe->CurveNumber; + assert(i!=nc); + curves[i].next=curves[nc].next; + curves[i].master=false; + curves[nc].next=curves+i; + if(be->ReverseEqui()) + curves[i].Reverse(); + } + } + + if(verbosity>3) + cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <4) + for(int i=0;i7) - for (i=0;ir - edges[i].v[0]->r; - Real8 lv10 = Norme2(v10); - if(lv10 == 0) { - cerr << "The length of " <9) - cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<=0) - { - ord++; - register Int4 i1= n /2; - register Int4 j1 = n % 2; - register Int4 *pn = ev + n; - float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; - n = *pn; - if (angleold > angle) // exch to have : po -> pn -> p - exch=1,*pn = *po,*po=*p,*p=n,po = pn; - else // to have : po -> p -> pn - angleold = angle, po = p,p = pn; - } - } // end while (exch) - - if (ord >= 1 ) - { /* - Int4 n = hv[i]; - while ( n >=0) - { Int4 i1 = n/2,j1 = n%2; - //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi; - n = ev[n]; - } - */ - } - if(ord == 2) { // angulare test to find a corner - Int4 n1 = hv[i]; - Int4 n2 = ev[n1]; - Int4 i1 = n1 /2, i2 = n2/2; // edge number - Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge - float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; - float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; - float da12 = Abs(angle2-angle1); - if(verbosity>9) - cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi - << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl; - - if (( da12 >= MaximalAngleOfCorner ) - && (da12 <= 2*Pi -MaximalAngleOfCorner)) { - vertices[i].SetCorner() ; - if(verbosity>7) - cout << " The vertex " << i << " is a corner (angle) " - << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;} - // if the ref a changing then is SetRequired(); - - if (edges[i1].flag != edges[i2].flag || edges[i1].Required()) - { - vertices[i].SetRequired(); - if(verbosity>7) - cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;} - - if (edges[i1].ref != edges[i2].ref) { - vertices[i].SetRequired(); - if(verbosity>7) - cout << " The vertex " << i << " is Required ref" << endl;} - } ; - - if(ord != 2) { - vertices[i].SetCorner(); - if(verbosity>7) - cout << " the vertex " << i << " is a corner ordre = " << ord << endl; - } - // close the liste around the vertex - { Int4 no=-1, ne = hv[i]; - while ( ne >=0) - ne = ev[no=ne]; - if(no>=0) - ev[no] = hv[i]; - } // now the list around the vertex is circular - - } // end for (i=0;i " - << Number(edges[i1].v[j1]) <10) - cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl; - } - - // generation of all the tangente - for (i=0;ir -edges[i].v[0]->r; - Real8 lAB = Norme2(AB); // length of current edge AB - Real8 ltg2[2]; - ltg2[0]=0;ltg2[1]=0; - for (jj=0;jj<2;jj++) { - R2 tg = edges[i].tg[jj]; - Real8 ltg = Norme2(tg); // length of tg - if(ltg == 0) {// no tg - if( ! edges[i].v[jj]->Corner()) { // not a Corner - tg = edges[i].v[1-jj]->r - - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r; - ltg = Norme2(tg); - tg = tg *(lAB/ltg),ltg=lAB; - /* - if(edges[i].ref >=4) - cout << " tg " << tg.x << " "<< tg.y << " " << edges[i].v[1-jj]->r << edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r << " y-y = " - << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl; - */ - } - - //else ;// a Corner with no tangent => nothing to do - } // a tg - else - tg = tg *(lAB/ltg),ltg=lAB; - ltg2[jj] = ltg; - if ( (tg,AB) < 0) - tg = -tg; - //if(edges[i].ref >=4) cout << " tg = " << tg << endl; - edges[i].tg[jj] = tg; - } // for (jj=0;jj<2;jj++) - - if (ltg2[0]!=0) edges[i].SetTgA(); - if (ltg2[1]!=0) edges[i].SetTgB(); - } // for (i=0;i7) - for (i=0;iSetMark(); - nbgem++; - e->CurveNumber=NbOfCurves; - if(curves) { - curves[NbOfCurves].ee=e; - curves[NbOfCurves].ke=k1; - } - - GeometricalVertex *b=(*e)(k1); - if (a == b || b->Required() ) break; - k0 = e->SensAdj[k1];// vertex in next edge - e = e->Adj[k1]; // next edge - - }// for(;;) - if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl; - NbOfCurves++; - if(level) { - if(verbosity>4) - cout << " Warning: Curve "<< NbOfCurves << " without required vertex " - << "so the vertex " << Number(a) << " become required " <SetRequired(); - } - - }} - assert(nbgem && nbe); - - if(step==0) { - curves = new Curve[NbOfCurves]; - } - } - for(int i=0;ilink; - //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; - curves[i].master=true; - if(be->Equi() || be->ReverseEqui() ) - { - assert(eqbe); - int nc = eqbe->CurveNumber; - assert(i!=nc); - curves[i].next=curves[nc].next; - curves[i].master=false; - curves[nc].next=curves+i; - if(be->ReverseEqui()) - curves[i].Reverse(); - } - } - - if(verbosity>3) - cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <4) - for(int i=0;i +#include +#include + +#include "../Mesh2.h" +#include "../QuadTree.h" + +namespace bamg { + +#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))) + +#define IJ(i,j,l) ( ( j & l) ? (( i & l) ? 3 : 2 ) :( ( i & l)? 1 : 0 )) +#define I_IJ(k,l) (( k&1) ? l : 0) +#define J_IJ(k,l) (( k&2) ? l : 0) + + /*Constructors/Destructors*/ + /*FUNCTION QuadTree::QuadTree(Triangles * t,long nbv){{{1*/ + QuadTree::QuadTree(Triangles * t,long nbv) : + lenStorageQuadTreeBox(t->nbvx/8+10), + th(t), + NbQuadTreeBox(0), + NbVertices(0), + NbQuadTreeBoxSearch(0), + NbVerticesSearch(0) + { + if (nbv == -1) nbv = t->nbv; + sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); + root=NewQuadTreeBox(); + assert( MaxISize > MaxICoor); + for (Int4 i=0;ivertices[i]); + } + /*}}}1*/ + /*FUNCTION QuadTree::QuadTree(){{{1*/ + QuadTree::QuadTree() : + lenStorageQuadTreeBox(100), + th(0), + NbQuadTreeBox(0), + NbVertices(0), + NbQuadTreeBoxSearch(0), + NbVerticesSearch(0) + { + sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); + root=NewQuadTreeBox(); + } + /*}}}1*/ + /*FUNCTION QuadTree::~QuadTree(){{{1*/ + QuadTree::~QuadTree() { + delete sb; + root=0; + } + /*}}}1*/ + + /*Methods*/ + /*FUNCTION QuadTree::NearestVertex{{{1*/ + Vertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { + QuadTreeBox * pb[ MaxDeep ]; + int pi[ MaxDeep ]; + Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; + register int l=0; // level + register QuadTreeBox * b; + IntQuad h=MaxISize,h0; + IntQuad hb = MaxISize; + Icoor1 i0=0,j0=0; + Icoor1 iplus( in) + return vn; // empty tree + + while( (n0 = b->n) < 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->n == 0) ) + break; // null box or empty => break + NbQuadTreeBoxSearch++; + 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; + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k];} + NbVerticesSearch++; + } + return vn; + } + // general case ----- + pb[0]= b; + pi[0]=b->n>0 ?(int) b->n : 4 ; + ii[0]=i0; + jj[0]=j0; + h=hb; + do { + b= pb[l]; + while (pi[l]--) + { + register int k = pi[l]; + + if (b->n>0) // Vertex QuadTreeBox none empty + { + NbVerticesSearch++; + I2 i2 = b->v[k]->i; + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k]; + } + } + else // Pointer QuadTreeBox + { + register QuadTreeBox *b0=b; + NbQuadTreeBoxSearch++; + 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->n>0 ?(int) b->n : 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::NearestVertexWithNormal{{{1*/ + Vertex* QuadTree::NearestVertexWithNormal(Icoor1 i,Icoor1 j) { + QuadTreeBox * pb[ MaxDeep ]; + int pi[ MaxDeep ]; + Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; + int l; // level + QuadTreeBox * b; + IntQuad h=MaxISize,h0; + IntQuad hb = MaxISize; + Icoor1 i0=0,j0=0; + Icoor1 iplus( in) + return vn; // empty tree + + while( (n0 = b->n) < 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->n == 0) ) + break; // null box or empty => break + NbQuadTreeBoxSearch++; + 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 sens -- + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k];} + NbVerticesSearch++; + } + if (vn) return vn; + } + // general case ----- + // INITIALISATION OF THE HEAP + l =0; // level + pb[0]= b; + pi[0]=b->n>0 ?(int) b->n : 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->n>0) // Vertex QuadTreeBox none empty + { + NbVerticesSearch++; + I2 i2 = b->v[k]->i; + // if good sens when try -- + + h0 = NORM(iplus,i2.x,jplus,i2.y); + if (h0 v[k]; + } + } + else // Pointer QuadTreeBox + { + register QuadTreeBox *b0=b; + NbQuadTreeBoxSearch++; + 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->n>0 ?(int) b->n : 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::ToClose {{{1*/ + Vertex * QuadTree::ToClose(Vertex & v,Real8 seuil,Icoor1 hx,Icoor1 hy){ + const Icoor1 i=v.i.x; + const Icoor1 j=v.i.y; + const R2 X(v.r); + const Metric Mx(v.m); + + QuadTreeBox * pb[ MaxDeep ]; + int pi[ MaxDeep ]; + Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; + register int l=0; // level + register QuadTreeBox * b; + Icoor1 h=MaxISize; + Icoor1 hb = MaxISize; + Icoor1 i0=0,j0=0; + + // Vertex *vn=0; + + if (!root->n) + return 0; // empty tree + + // general case ----- + pb[0]=root; + pi[0]=root->n>0 ?(int) root->n : 4 ; + ii[0]=i0; + jj[0]=j0; + h=hb; + do { + b= pb[l]; + while (pi[l]--) + { + register int k = pi[l]; + + if (b->n>0) // Vertex QuadTreeBox none empty + { + NbVerticesSearch++; + I2 i2 = b->v[k]->i; + if ( ABS(i-i2.x) v[k]->r); + Real8 dd; + // old code if( Mx(XY) + b->v[k]->m(XY) < seuil ) + if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ) + { + // cout << CurrentTh->Number(v) << "is To Close " + // << CurrentTh->Number( b->v[k]) << " l=" <v[k]; + } + } + } + else // Pointer QuadTreeBox + { + register QuadTreeBox *b0=b; + NbQuadTreeBoxSearch++; + 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->n>0 ?(int) b->n : 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*/ + /*FUNCTION QuadTree::Add{{{1*/ + void QuadTree::Add( Vertex & w){ + QuadTreeBox ** pb , *b; + register long i=w.i.x, j=w.i.y,l=MaxISize; + pb = &root; + // cout << pb << " " << &root << endl; + while( (b=*pb) && (b->n<0)) + { + b->n--; + l >>= 1; + pb = &b->b[IJ(i,j,l)]; + } + if (b) { + if (b->n > 3 && b->v[3] == &w) return; + if (b->n > 2 && b->v[2] == &w) return; + if (b->n > 1 && b->v[1] == &w) return; + if (b->n > 0 && b->v[0] == &w) return; + } + assert(l); + while ((b= *pb) && (b->n == 4)) // the QuadTreeBox is full + { + Vertex *v4[4]; // copy of the QuadTreeBox vertices + + v4[0]= b->v[0]; + v4[1]= b->v[1]; + v4[2]= b->v[2]; + v4[3]= b->v[3]; + b->n = -b->n; // mark is pointer QuadTreeBox + b->b[0]=b->b[1]=b->b[2]=b->b[3]=0; // set empty QuadTreeBox ptr + l >>= 1; // div the size by 2 + for (register int k=0;k<4;k++) // for the 4 vertices find the sub QuadTreeBox ij + { + register int ij; + register QuadTreeBox * bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,l)]; + if (!bb) + bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox + // cout << bb << " " << k << " " << ij << endl; + bb->v[bb->n++] = v4[k]; + } + pb = &b->b[IJ(i,j,l)]; + } + if (!(b = *pb)) + b=*pb= NewQuadTreeBox(); // alloc the QuadTreeBox + // cout << b << " " << b->n << endl; + b->v[b->n++]=&w; // we add the vertex + NbVertices++; + } + /*}}}1*/ + /*FUNCTION QuadTree::StorageQuadTreeBox::StorageQuadTreeBox{{{1*/ + QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) { + len = ll; + n = nn; + b = new QuadTreeBox[ll]; + for (int i = 0; i