Index: /issm/trunk/src/c/Bamgx/QuadTree.h =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2884) +++ /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2885) @@ -8,16 +8,16 @@ typedef long IntQuad; const IntQuad MaxISize = ( 1L << MaxDeep); - class Triangles; class Vertex; - class QuadTree { public: class QuadTreeBox { public: - long n; // if n < 4 => Vertex else => QuadTreeBox; - union { - QuadTreeBox *b[4]; - Vertex * v[4]; + long n; + union{ + //contains only one object form the list (either Vertex or QuadTreeBox) + // if n < 4 => Vertex else => QuadTreeBox; + QuadTreeBox* b[4]; + Vertex* v[4]; }; }; // end class @@ -26,15 +26,13 @@ QuadTreeBox *b,*bc,*be; long len; - StorageQuadTreeBox *n; // next StorageQuadTreeBox - StorageQuadTreeBox(long ,StorageQuadTreeBox * =0); + StorageQuadTreeBox* n; // next StorageQuadTreeBox + StorageQuadTreeBox(long ,StorageQuadTreeBox* =NULL); ~StorageQuadTreeBox() { if(n) delete n; delete [] b; } - long SizeOf() const { - return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0); - } + long SizeOf() const {return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);} }; // end class - StorageQuadTreeBox * sb; + StorageQuadTreeBox* sb; long lenStorageQuadTreeBox; @@ -51,7 +49,5 @@ QuadTreeBox* NewQuadTreeBox(){ if(! (sb->bcbe)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); - if (!sb || (sb->bc->n != 0)){ - throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)")); - } + if (!sb || (sb->bc->n != 0)){throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));} NbQuadTreeBox++; return sb->bc++; Index: /issm/trunk/src/c/Bamgx/meshtype.h =================================================================== --- /issm/trunk/src/c/Bamgx/meshtype.h (revision 2884) +++ /issm/trunk/src/c/Bamgx/meshtype.h (revision 2885) @@ -1,29 +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 - */ #ifndef MESHTYPE_H #define MESHTYPE_H @@ -52,12 +25,4 @@ const Icoor1 MaxICoor = 1073741823; // 2^30-1 const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; - -#elif defined(BAMG_LONG_LONG) - typedef long Icoor1; - typedef long long Icoor2; - const Icoor1 MaxICoor = 1073741823;// 2^30-1 -// not a const due to a bug in hp compiler -#define MaxICoor22 2305843004918726658LL -//const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; #else typedef int Icoor1; Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2884) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2885) @@ -268,5 +268,6 @@ vertices[i].Set(); } - //find domain extrema for pmin,pmax + //find domain extrema for pmin,pmax that will define the square box + //used for by the quadtree pmin = vertices[0].r; pmax = vertices[0].r; @@ -280,4 +281,6 @@ pmin -= DD05; pmax += DD05; + //coefIcoor is the coefficient used to have coordinates + //in ]0 1[^2: x'=coefIcoor*(x-pmin.x) coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); if(coefIcoor <=0){ @@ -493,5 +496,5 @@ float* eangle=new float[nbe]; double eps=1e-20; - QuadTree quadtree; // to find same vertices + QuadTree quadtree; // build quadtree to find duplicates Vertex* v0=vertices; GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0; @@ -499,15 +502,24 @@ k=0; //link all vertices to themselves by default - for (i=0;i0)) for (i=0;i [0 2^(30-1)[^2 for a quadtree of depth 30*/ vertices[i].i=toI2(vertices[i].r); - Vertex* v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); + + //Build the quadtree: + Vertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); if( v && Norme1(v->r - vertices[i]) < eps ){ // mama's old trick to get j - GeometricalVertex* vg = (GeometricalVertex*) (void*)v; + GeometricalVertex* vg=(GeometricalVertex*) (void*)v; j=vg-v0g; + //check that the clostest vertex is not itself... if ( v != &(Vertex &) vertices[j]){ throw ErrorException(__FUNCT__,exprintf(" v != &(Vertex &) vertices[j]")); @@ -516,5 +528,5 @@ k++; } - else quadtree.Add(vertices[i]); + else quadtree.Add(vertices[i]); //Add vertices[i] to the quadtree } if (k) { @@ -548,34 +560,40 @@ } - //build hv and ev - for (i=0;ir - edges[i].v[0]->r; - Real8 lv10 = Norme2(v10); - if(lv10 == 0) { + Real8 lv10=Norme2(v10); + //check that its length is not 0 + if(lv10==0) { throw ErrorException(__FUNCT__,exprintf("Length of edge %i is 0",i)); } - eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ] - for (jj=0;jj<2;jj++){ - Int4 v=Number(edges[i].v[jj]); + //compute angle in [-Pi Pi] + eangle[i] = atan2(v10.y,v10.x); + //build hv and ev + for (j=0;j<2;j++){ + Int4 v=Number(edges[i].v[j]); ev[k]=hv[v]; hv[v]=k++; - } - } + } + } + // bulle sort on the angle of edge for (i=0;i=0) - { + while (n >=0){ ord++; - register Int4 i1= n /2; - register Int4 j1 = n % 2; - register Int4 *pn = ev + n; + 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; @@ -584,14 +602,15 @@ else // to have : po -> p -> pn angleold = angle, po = p,p = pn; - } - } // end while (exch) - - if(ord == 2) { // angulare test to find a corner + } + } + + // angulare test to find a corner + if(ord == 2) { 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 angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; + float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; float da12 = Abs(angle2-angle1); if (( da12 >= MaximalAngleOfCorner ) && (da12 <= 2*Pi -MaximalAngleOfCorner)) { @@ -606,29 +625,27 @@ } } - if(ord != 2) { vertices[i].SetCorner(); } + // 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=0) ne = ev[no=ne]; + if(no>=0) ev[no] = hv[i]; + // now the list around the vertex is circular + } k =0; - for (i=0;i nothing to do - } // a tg - else - tg = tg *(lAB/ltg),ltg=lAB; + //else: a Corner with no tangent => nothing to do + } + else{ + tg = tg *(lAB/ltg),ltg=lAB; + } ltg2[jj] = ltg; - if ( (tg,AB) < 0) - tg = -tg; + if ((tg,AB)<0) tg = -tg; 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;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->DirAdj[k1];// vertex in next edge - e = e->Adj[k1]; // next edge - - }// for(;;) - NbOfCurves++; - if(level) { - a->SetRequired(); - } - - } + for(jj=0;jj<2;jj++){ + if (!ei.Mark() && (level || ei[jj].Required())) { + // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++) + int k0=jj,k1; + GeometricalEdge *e = & ei; + GeometricalVertex *a=(*e)(k0); // begin + if(curves) { + curves[NbOfCurves].be=e; + curves[NbOfCurves].kb=k0; + } + int nee=0; + for(;;) { + nee++; + k1 = 1-k0; // next vertex of the edge + e->SetMark(); + 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->DirAdj[k1];// vertex in next edge + e = e->Adj[k1]; // next edge + } + NbOfCurves++; + if(level) a->SetRequired(); + } + } } if (nbgem==0 || nbe==0){ throw ErrorException(__FUNCT__,exprintf("nbgem==0 || nbe==0")); } - if(step==0) { curves = new Curve[NbOfCurves]; @@ -734,4 +743,5 @@ } + /*clean up*/ delete []ev; delete []hv; Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2884) +++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2885) @@ -65,9 +65,10 @@ /*Methods*/ /*FUNCTION QuadTree::Add{{{1*/ - void QuadTree::Add( Vertex & w){ - QuadTreeBox ** pb , *b; + void QuadTree::Add(Vertex &w){ + QuadTreeBox** pb; + QuadTreeBox* b; register long i=w.i.x, j=w.i.y,l=MaxISize; pb = &root; - while( (b=*pb) && (b->n<0)){ + while((b=*pb) && (b->n<0)){ b->n--; l >>= 1; @@ -109,18 +110,21 @@ /*FUNCTION QuadTree::NearestVertex{{{1*/ Vertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { + + /*Build QuadTree*/ QuadTreeBox* pb[MaxDeep]; - int pi[MaxDeep]; - Icoor1 ii[MaxDeep], jj [MaxDeep]; + int pi[MaxDeep]; + Icoor1 ii[MaxDeep]; + Icoor1 jj[MaxDeep]; register int l=0; // level - register QuadTreeBox * b; + register QuadTreeBox* b; IntQuad h=MaxISize,h0; - IntQuad hb = MaxISize; - Icoor1 i0=0,j0=0; - Icoor1 iplus( i