![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2885
- Timestamp:
- 01/21/10 17:44:35 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 4 edited
-
QuadTree.h (modified) (3 diffs)
-
meshtype.h (modified) (2 diffs)
-
objects/Geometry.cpp (modified) (11 diffs)
-
objects/QuadTree.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/QuadTree.h
r2877 r2885 8 8 typedef long IntQuad; 9 9 const IntQuad MaxISize = ( 1L << MaxDeep); 10 11 10 class Triangles; 12 11 class Vertex; 13 14 12 class QuadTree { 15 13 public: 16 14 class QuadTreeBox { 17 15 public: 18 long n; // if n < 4 => Vertex else => QuadTreeBox; 19 union { 20 QuadTreeBox *b[4]; 21 Vertex * v[4]; 16 long n; 17 union{ 18 //contains only one object form the list (either Vertex or QuadTreeBox) 19 // if n < 4 => Vertex else => QuadTreeBox; 20 QuadTreeBox* b[4]; 21 Vertex* v[4]; 22 22 }; 23 23 }; // end class … … 26 26 QuadTreeBox *b,*bc,*be; 27 27 long len; 28 StorageQuadTreeBox *n; // next StorageQuadTreeBox29 StorageQuadTreeBox(long ,StorageQuadTreeBox * =0);28 StorageQuadTreeBox* n; // next StorageQuadTreeBox 29 StorageQuadTreeBox(long ,StorageQuadTreeBox* =NULL); 30 30 ~StorageQuadTreeBox() { 31 31 if(n) delete n; 32 32 delete [] b; 33 33 } 34 long SizeOf() const { 35 return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0); 36 } 34 long SizeOf() const {return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0);} 37 35 }; // end class 38 StorageQuadTreeBox 36 StorageQuadTreeBox* sb; 39 37 long lenStorageQuadTreeBox; 40 38 … … 51 49 QuadTreeBox* NewQuadTreeBox(){ 52 50 if(! (sb->bc<sb->be)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); 53 if (!sb || (sb->bc->n != 0)){ 54 throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)")); 55 } 51 if (!sb || (sb->bc->n != 0)){throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)"));} 56 52 NbQuadTreeBox++; 57 53 return sb->bc++; -
issm/trunk/src/c/Bamgx/meshtype.h
r2851 r2885 1 // -*- Mode : c++ -*-2 //3 // SUMMARY :4 // USAGE :5 // ORG :6 // AUTHOR : Frederic Hecht7 // E-MAIL : hecht@ann.jussieu.fr8 //9 10 /*11 12 This file is part of Freefem++13 14 Freefem++ is free software; you can redistribute it and/or modify15 it under the terms of the GNU Lesser General Public License as published by16 the Free Software Foundation; either version 2.1 of the License, or17 (at your option) any later version.18 19 Freefem++ is distributed in the hope that it will be useful,20 but WITHOUT ANY WARRANTY; without even the implied warranty of21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the22 GNU Lesser General Public License for more details.23 24 You should have received a copy of the GNU Lesser General Public License25 along with Freefem++; if not, write to the Free Software26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA27 */28 1 #ifndef MESHTYPE_H 29 2 #define MESHTYPE_H … … 52 25 const Icoor1 MaxICoor = 1073741823; // 2^30-1 53 26 const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ; 54 55 #elif defined(BAMG_LONG_LONG)56 typedef long Icoor1;57 typedef long long Icoor2;58 const Icoor1 MaxICoor = 1073741823;// 2^30-159 // not a const due to a bug in hp compiler60 #define MaxICoor22 2305843004918726658LL61 //const Icoor2 MaxICoor22 = Icoor2(2)*Icoor2(MaxICoor) * Icoor2(MaxICoor) ;62 27 #else 63 28 typedef int Icoor1; -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2877 r2885 268 268 vertices[i].Set(); 269 269 } 270 //find domain extrema for pmin,pmax 270 //find domain extrema for pmin,pmax that will define the square box 271 //used for by the quadtree 271 272 pmin = vertices[0].r; 272 273 pmax = vertices[0].r; … … 280 281 pmin -= DD05; 281 282 pmax += DD05; 283 //coefIcoor is the coefficient used to have coordinates 284 //in ]0 1[^2: x'=coefIcoor*(x-pmin.x) 282 285 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 283 286 if(coefIcoor <=0){ … … 493 496 float* eangle=new float[nbe]; 494 497 double eps=1e-20; 495 QuadTree quadtree; // to find same vertices498 QuadTree quadtree; // build quadtree to find duplicates 496 499 Vertex* v0=vertices; 497 500 GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0; … … 499 502 k=0; 500 503 //link all vertices to themselves by default 501 for (i=0;i<nbv;i++) vertices[i].link = vertices 504 for (i=0;i<nbv;i++) vertices[i].link = vertices+i; 502 505 503 506 //build quadtree for this geometry (error if we have duplicates (k>0)) 504 507 for (i=0;i<nbv;i++){ 505 // set integer coordinate 508 509 /*build integer coordinates (non unique) 510 these coordinates are used by the quadtree to group 511 the vertices by groups of 5: 512 All the coordinates are transformed to ]0,1[^2 513 then, the integer coordinates are computed using 514 the transformation ]0,1[^2 -> [0 2^(30-1)[^2 for a quadtree of depth 30*/ 506 515 vertices[i].i=toI2(vertices[i].r); 507 Vertex* v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 516 517 //Build the quadtree: 518 Vertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 508 519 if( v && Norme1(v->r - vertices[i]) < eps ){ 509 520 // mama's old trick to get j 510 GeometricalVertex* vg =(GeometricalVertex*) (void*)v;521 GeometricalVertex* vg=(GeometricalVertex*) (void*)v; 511 522 j=vg-v0g; 523 //check that the clostest vertex is not itself... 512 524 if ( v != &(Vertex &) vertices[j]){ 513 525 throw ErrorException(__FUNCT__,exprintf(" v != &(Vertex &) vertices[j]")); … … 516 528 k++; 517 529 } 518 else quadtree.Add(vertices[i]); 530 else quadtree.Add(vertices[i]); //Add vertices[i] to the quadtree 519 531 } 520 532 if (k) { … … 548 560 } 549 561 550 //build hv and ev 551 for (i=0;i<nbv;i++) hv[i]=-1;// empty list 562 //initialize hv as -1 for all vertices 563 for (i=0;i<nbv;i++) hv[i]=-1; 564 k=0; 552 565 for (i=0;i<nbe;i++) { 566 //compute vector of edge i that goes from vertex 0 to vertex 1 553 567 R2 v10=edges[i].v[1]->r - edges[i].v[0]->r; 554 Real8 lv10 = Norme2(v10); 555 if(lv10 == 0) { 568 Real8 lv10=Norme2(v10); 569 //check that its length is not 0 570 if(lv10==0) { 556 571 throw ErrorException(__FUNCT__,exprintf("Length of edge %i is 0",i)); 557 572 } 558 eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ] 559 for (jj=0;jj<2;jj++){ 560 Int4 v=Number(edges[i].v[jj]); 573 //compute angle in [-Pi Pi] 574 eangle[i] = atan2(v10.y,v10.x); 575 //build hv and ev 576 for (j=0;j<2;j++){ 577 Int4 v=Number(edges[i].v[j]); 561 578 ev[k]=hv[v]; 562 579 hv[v]=k++; 563 } 564 } 580 } 581 } 582 565 583 // bulle sort on the angle of edge 566 584 for (i=0;i<nbv;i++) { 567 int exch = 1,ord =0; 568 while (exch) { 569 exch = 0; 570 Int4 *p = hv + i, *po = p; 571 Int4 n = *p; 572 register float angleold = -1000 ; // angle = - infini 585 int exch=1, ord=0; 586 while (exch){ 587 exch =0; 588 Int4* p=hv+i; 589 Int4* po=p; 590 Int4 n=*p; 591 register float angleold=-1000 ; // angle = - infinity 573 592 ord = 0; 574 while (n >=0) 575 { 593 while (n >=0){ 576 594 ord++; 577 register Int4 i1= n/2;578 register Int4 j1 = n %2;579 register Int4 *pn = ev +n;595 register Int4 i1=n/2; 596 register Int4 j1=n%2; 597 register Int4 *pn=ev+n; 580 598 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; 581 599 n = *pn; … … 584 602 else // to have : po -> p -> pn 585 603 angleold = angle, po = p,p = pn; 586 } 587 } // end while (exch) 588 589 if(ord == 2) { // angulare test to find a corner 604 } 605 } 606 607 // angulare test to find a corner 608 if(ord == 2) { 590 609 Int4 n1 = hv[i]; 591 610 Int4 n2 = ev[n1]; 592 611 Int4 i1 = n1 /2, i2 = n2/2; // edge number 593 612 Int4 j1 = n1 %2, j2 = n2%2; // vertex in the edge 594 float angle1= j1 ? OppositeAngle(eangle[i1]) :eangle[i1];595 float angle2= !j2 ? OppositeAngle(eangle[i2]) : 613 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; 614 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; 596 615 float da12 = Abs(angle2-angle1); 597 616 if (( da12 >= MaximalAngleOfCorner ) && (da12 <= 2*Pi -MaximalAngleOfCorner)) { … … 606 625 } 607 626 } 608 609 627 if(ord != 2) { 610 628 vertices[i].SetCorner(); 611 629 } 630 612 631 // close the liste around the vertex 613 { Int4 no=-1, ne = hv[i]; 614 while ( ne >=0) 615 ne = ev[no=ne]; 616 if(no>=0) 617 ev[no] = hv[i]; 618 } // now the list around the vertex is circular 619 620 } // end for (i=0;i<nbv;i++) 632 Int4 no=-1, ne = hv[i]; 633 while (ne >=0) ne = ev[no=ne]; 634 if(no>=0) ev[no] = hv[i]; 635 // now the list around the vertex is circular 636 } 621 637 622 638 k =0; 623 for (i=0;i<nbe;i++) 624 for (jj=0;jj<2;jj++){ 625 Int4 n1 = ev[k++]; 626 Int4 i1 = n1/2 ,j1=n1%2; 627 if( edges[i1].v[j1] != edges[i].v[jj]) { 628 throw ErrorException(__FUNCT__,exprintf("Bug Adj edge")); 629 } 630 edges[i1].Adj[j1] = edges + i; 631 edges[i1].DirAdj[j1] = jj; 632 } 639 for (i=0;i<nbe;i++){ 640 for (jj=0;jj<2;jj++){ 641 Int4 n1 = ev[k++]; 642 Int4 i1 = n1/2 ,j1=n1%2; 643 if( edges[i1].v[j1] != edges[i].v[jj]) { 644 throw ErrorException(__FUNCT__,exprintf("Bug Adj edge")); 645 } 646 edges[i1].Adj[j1] = edges + i; 647 edges[i1].DirAdj[j1] = jj; 648 } 649 } 633 650 634 651 // generation of all the tangente … … 648 665 tg = tg *(lAB/ltg),ltg=lAB; 649 666 } 650 651 //else ;// a Corner with no tangent => nothing to do652 } // a tg653 else654 tg = tg *(lAB/ltg),ltg=lAB;667 //else: a Corner with no tangent => nothing to do 668 } 669 else{ 670 tg = tg *(lAB/ltg),ltg=lAB; 671 } 655 672 ltg2[jj] = ltg; 656 if ( (tg,AB) < 0) 657 tg = -tg; 673 if ((tg,AB)<0) tg = -tg; 658 674 edges[i].tg[jj] = tg; 659 } // for (jj=0;jj<2;jj++) 660 675 } 661 676 if (ltg2[0]!=0) edges[i].SetTgA(); 662 677 if (ltg2[1]!=0) edges[i].SetTgB(); 663 } // for (i=0;i<nbe;i++)678 } 664 679 665 680 for (int step=0;step<2;step++){ 666 681 for (i=0;i<nbe;i++) edges[i].SetUnMark(); 667 668 682 NbOfCurves = 0; 669 683 Int4 nbgem=0; … … 671 685 for (i=0;i<nbe;i++) { 672 686 GeometricalEdge & ei = edges[i]; 673 for(jj=0;jj<2;jj++) 674 if (!ei.Mark() && (level || ei[jj].Required())) { 675 // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++) 676 int k0=jj,k1; 677 GeometricalEdge *e = & ei; 678 GeometricalVertex *a=(*e)(k0); // begin 679 if(curves) { 680 curves[NbOfCurves].be=e; 681 curves[NbOfCurves].kb=k0; 682 } 683 int nee=0; 684 for(;;) { 685 nee++; 686 k1 = 1-k0; // next vertex of the edge 687 e->SetMark(); 688 nbgem++; 689 e->CurveNumber=NbOfCurves; 690 if(curves) { 691 curves[NbOfCurves].ee=e; 692 curves[NbOfCurves].ke=k1; 693 } 694 695 GeometricalVertex *b=(*e)(k1); 696 if (a == b || b->Required() ) break; 697 k0 = e->DirAdj[k1];// vertex in next edge 698 e = e->Adj[k1]; // next edge 699 700 }// for(;;) 701 NbOfCurves++; 702 if(level) { 703 a->SetRequired(); 704 } 705 706 } 687 for(jj=0;jj<2;jj++){ 688 if (!ei.Mark() && (level || ei[jj].Required())) { 689 // warning ei.Mark() can be change in loop for(jj=0;jj<2;jj++) 690 int k0=jj,k1; 691 GeometricalEdge *e = & ei; 692 GeometricalVertex *a=(*e)(k0); // begin 693 if(curves) { 694 curves[NbOfCurves].be=e; 695 curves[NbOfCurves].kb=k0; 696 } 697 int nee=0; 698 for(;;) { 699 nee++; 700 k1 = 1-k0; // next vertex of the edge 701 e->SetMark(); 702 nbgem++; 703 e->CurveNumber=NbOfCurves; 704 if(curves) { 705 curves[NbOfCurves].ee=e; 706 curves[NbOfCurves].ke=k1; 707 } 708 GeometricalVertex *b=(*e)(k1); 709 if (a == b || b->Required() ) break; 710 k0 = e->DirAdj[k1];// vertex in next edge 711 e = e->Adj[k1]; // next edge 712 } 713 NbOfCurves++; 714 if(level) a->SetRequired(); 715 } 716 } 707 717 } 708 718 if (nbgem==0 || nbe==0){ 709 719 throw ErrorException(__FUNCT__,exprintf("nbgem==0 || nbe==0")); 710 720 } 711 712 721 if(step==0) { 713 722 curves = new Curve[NbOfCurves]; … … 734 743 } 735 744 745 /*clean up*/ 736 746 delete []ev; 737 747 delete []hv; -
issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
r2877 r2885 65 65 /*Methods*/ 66 66 /*FUNCTION QuadTree::Add{{{1*/ 67 void QuadTree::Add( Vertex & w){ 68 QuadTreeBox ** pb , *b; 67 void QuadTree::Add(Vertex &w){ 68 QuadTreeBox** pb; 69 QuadTreeBox* b; 69 70 register long i=w.i.x, j=w.i.y,l=MaxISize; 70 71 pb = &root; 71 while( 72 while((b=*pb) && (b->n<0)){ 72 73 b->n--; 73 74 l >>= 1; … … 109 110 /*FUNCTION QuadTree::NearestVertex{{{1*/ 110 111 Vertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { 112 113 /*Build QuadTree*/ 111 114 QuadTreeBox* pb[MaxDeep]; 112 int pi[MaxDeep]; 113 Icoor1 ii[MaxDeep], jj [MaxDeep]; 115 int pi[MaxDeep]; 116 Icoor1 ii[MaxDeep]; 117 Icoor1 jj[MaxDeep]; 114 118 register int l=0; // level 115 register QuadTreeBox 119 register QuadTreeBox* b; 116 120 IntQuad h=MaxISize,h0; 117 IntQuad hb =MaxISize;118 Icoor1 i0=0,j0=0;119 Icoor1 iplus( i<MaxISize?(i<0?0:i):MaxISize-1);120 Icoor1 jplus( j<MaxISize?(j<0?0:j):MaxISize-1);121 IntQuad hb=MaxISize; 122 Icoor1 i0=0,j0=0; 123 Icoor1 iplus( i<MaxISize?(i<0?0:i):MaxISize-1); 124 Icoor1 jplus( j<MaxISize?(j<0?0:j):MaxISize-1); 121 125 122 126 Vertex *vn=0; 123 127 124 // initfor optimization128 //initialization for optimization 125 129 b = root; 126 130 register Int4 n0; … … 188 192 ii[l]= iii; 189 193 jj[l]= jjj; 190 191 194 } 192 195 else{
Note:
See TracChangeset
for help on using the changeset viewer.
