![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2856
- Timestamp:
- 01/15/10 12:26:00 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 8 edited
-
Mesh2.h (modified) (3 diffs)
-
QuadTree.h (modified) (1 diff)
-
SetOfE4.h (modified) (3 diffs)
-
objects/GeometricalEdge.cpp (modified) (1 diff)
-
objects/Geometry.cpp (modified) (22 diffs)
-
objects/MetricAnIso.cpp (modified) (3 diffs)
-
objects/QuadTree.cpp (modified) (4 diffs)
-
objects/Triangles.cpp (modified) (32 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2855 r2856 500 500 IntersectionTriangles & operator[](int i) {return lIntTria[i];} 501 501 operator int&() {return Size;} 502 ListofIntersectionTriangles(int n=256,int m=16) 503 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , 504 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]) 505 { 506 long int verbosity=0; 507 508 if (verbosity>9)509 std::cout << " construct ListofIntersectionTriangles" 510 << MaxSize << " " << MaxNbSeg<< std::endl;}; 502 503 ListofIntersectionTriangles(int n=256,int m=16) 504 : MaxSize(n), Size(0), len(-1),state(-1),lIntTria(new IntersectionTriangles[n]) , 505 NbSeg(0), MaxNbSeg(m), lSegsI(new SegInterpolation[m]) 506 { 507 long int verbosity=0; 508 if (verbosity>9) printf(" construct ListofIntersectionTriangles %i %i\n",MaxSize,MaxNbSeg); 509 }; 510 511 511 ~ListofIntersectionTriangles(){ 512 512 if (lIntTria) delete [] lIntTria,lIntTria=0; … … 555 555 nw[i] = lIntTria[i]; 556 556 long int verbosity=0; 557 if(verbosity>3) 558 std::cout << " ListofIntersectionTriangles ReShape MaxSize " 559 << MaxSize << " -> " 560 << newsize << std::endl; 557 if(verbosity>3) printf(" ListofIntersectionTriangles ReShape Maxsize %i -> %i\n",MaxSize,MaxNbSeg); 561 558 MaxSize = newsize; 562 559 delete [] lIntTria;// remove old … … 603 600 VertexOnGeom(): mv(0),abscisse(0){gv=0;} 604 601 VertexOnGeom(Vertex & m,GeometricalVertex &g) : mv(&m),abscisse(-1){gv=&g;} 605 // std::cout << " mv = " <<mv << " gv = " << gv << std::endl;}606 602 VertexOnGeom(Vertex & m,GeometricalEdge &g,Real8 s) : mv(&m),abscisse(s){ge=&g;} 607 //std::cout << &g << " " << ge << std::endl;608 603 operator Vertex * () const {return mv;} 609 604 operator GeometricalVertex * () const {return gv;} 610 605 operator GeometricalEdge * () const {return ge;} 611 // operator Real8 & () {return abscisse;}612 606 operator const Real8 & () const {return abscisse;} 613 607 int IsRequiredVertex(){ return this? (( abscisse<0 ? (gv?gv->Required():0):0 )) : 0;} -
issm/trunk/src/c/Bamgx/QuadTree.h
r2854 r2856 34 34 StorageQuadTreeBox *n; // next StorageQuadTreeBox 35 35 StorageQuadTreeBox(long ,StorageQuadTreeBox * =0); 36 ~StorageQuadTreeBox() 37 { //cout << "~StorageQuadTreeBox " << this << " n " << n << " b " << b << endl; 36 ~StorageQuadTreeBox() { 38 37 if(n) delete n; 39 38 delete [] b; -
issm/trunk/src/c/Bamgx/SetOfE4.h
r2740 r2856 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 29 1 #ifndef _SetOfEdge4_h 30 2 #define _SetOfEdge4_h … … 47 19 public: 48 20 SetOfEdges4(Int4 ,Int4);// nb Edges mx , nb de sommet 49 ~SetOfEdges4() { // cout << " delete SetofArete " << endl ;21 ~SetOfEdges4() { 50 22 delete [] tete; delete [] Edges;} 51 23 Int4 add (Int4 ii,Int4 jj); … … 54 26 Int4 find (Int4 ii,Int4 jj); 55 27 Int4 findtrie (Int4 ii,Int4 jj) {return ii <=jj ? find (ii,jj) : find (jj,ii) ;} 56 // inline void close();57 28 Int4 i(Int4 k){return Edges[k].i;} 58 29 Int4 j(Int4 k){return Edges[k].j;} -
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
r2854 r2856 112 112 cta = (1-theta)*(1-theta)*theta; 113 113 ctb = (theta-1)*theta*theta ; 114 // if(ref==4 || ref==5)115 // cout << " FFF " << tg[0] << tg[1] << A << B << " => " << A*ca + B*cb + tg[0]* cta + tg[1] * ctb << endl;116 114 } 117 115 else { // 1-t*t, t-t*t, t*t -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2854 r2856 54 54 throw ErrorException(__FUNCT__,exprintf("NbRef>0")); 55 55 } 56 if(verbosity>9)57 cout << "DELETE ~Geometry "<< this << endl;58 56 if(vertices) delete [] vertices;vertices=0; 59 57 if(edges) delete [] edges;edges=0; … … 511 509 GeometricalEdge* on =start,* pon=0; 512 510 // walk with the cos on geometry 513 // cout << P ;514 511 int k=0; 515 512 while(pon != on){ … … 524 521 R2 AP = P-A; 525 522 R2 BP = P-B; 526 // cout << ":: " << on - edges << " " << AB*AP << " " << AB*BP << " " << A << B << endl;527 523 if ( (AB,AP) < 0) 528 524 on = on->Adj[0]; … … 560 556 R2 V0 = v0,V1=v1,V01=V1-V0; 561 557 VertexOnGeom vg0= *v0.on, vg1=*v1.on; 562 if(NbTry) cout << "bug: s==== " << s << " e=" << V0 << " " << V1 << endl;563 558 564 559 // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; 565 560 GeometricalEdge * eg0=on, *eg1=on; 566 561 R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; 567 if(NbTry) cout <<" G edge= " << Ag << Bg << endl << " v edge" << V01 << " v geom " << AB << (V01,AB) <<endl;568 562 int OppositeSens = (V01,AB) < 0; 569 563 int sens0=0,sens1=1; 570 564 if (OppositeSens) 571 565 s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); 572 if(NbTry) cout << "bug: edge = " << v0.r << " -> " << v1.r << endl573 << "sg 0 = " << vg0574 << " on = " << Number(on) << ":" << Ag << Bg << "; "575 << " sg 1= " << vg1576 << "--------------------------------------------" << endl;577 566 while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0){ 578 567 if (bge<=0) { … … 589 578 goto retry;} 590 579 GeometricalEdge* tmpge = eg0; 591 if(NbTry)592 cout << "bug: --Edge @" << Number(tmpge) << " = "<< Number(eg0) << ":" <<Number(eg0->Adj[0]) << "," <<593 Number(eg0->Adj[1]) <<"," ;594 580 ge[--bge] =eg0 = eg0->Adj[sens0]; 595 581 if (bge<0 || bge>mxe){ … … 597 583 } 598 584 sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); 599 if(NbTry)600 cout << "bug: Edge " << Number(eg0) << " "<< 1-sens0 << " S "601 << Number((*eg0)[1-sens0]) <<":" << Number(eg0->Adj[0]) << ","602 << Number(eg0->Adj[1]) <<"," << endl603 <<Number(eg0)<< (*eg0)[sens0].r << "v = " << Number((*eg1)(sens0)) << " e = " << eg0 << endl;604 585 } 605 if(NbTry) cout << Number((GeometricalEdge*) vg1) << " " << Number((GeometricalVertex*) vg1) << endl; 606 while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) 607 { 586 while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) { 608 587 if(tge>=mxe ) { 609 cerr << " --Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl;588 printf("WARNING: on the class triangles before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); 610 589 NbTry++; 611 590 if (NbTry<2) goto retry; … … 619 598 620 599 GeometricalEdge* tmpge = eg1; 621 if(NbTry)622 cout << "++Edge @" << tmpge << " = " << Number(eg1) <<"%" << Number(eg1->Adj[0]) << ","623 << Number(eg1->Adj[1]) <<"," ;624 600 ge[++tge] =eg1 = eg1->Adj[sens1]; 625 601 sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; … … 627 603 throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)")); 628 604 } 629 if(NbTry)630 cout << " Edge " << Number(eg1) << " " << sens1 << " S "631 <<Number((*eg1)[sens1]) <<"%"<< Number(eg1->Adj[0]) << "," << Number(eg1->Adj[1]) <<","632 <<Number(eg1)<< (*eg1)[sens1].r << "v = " << Number((*eg1)(sens1)) << " e = " << Number(eg1) << endl;633 605 } 634 635 if(NbTry) cout << endl;636 606 637 607 … … 643 613 644 614 Real8 sg; 645 // cout << " " << Number(on) << " " << Number(eg0) << " " << Number(eg1) << " " ;646 615 if (eg0 == eg1) { 647 616 register Real8 s0= vg0,s1=vg1; 648 617 sg = s0 * (1.0-s) + s * s1; 649 // cout <<" s0=" << s0 << " s1=" << s1650 // << " s = " << s << " sens= " << OppositeSens << "\t\t sg = " << sg << endl ;651 618 on=eg0;} 652 619 else { 653 620 R2 AA=V0,BB; 654 621 Real8 s0,s1; 655 656 //cout << endl << "s= " << s << Number(eg0) << " " << (Real8) vg0 << " "657 // << Number(eg1) << " " << (Real8) vg1 << V0 << V1 << " Interpol = "658 // << V0*(1-s)+V1*s << ";;; " << endl;659 622 int i=bge; 660 623 Real8 ll=0; … … 665 628 BB = (*ge[i])[sensge[i]]; 666 629 lge[i]=ll += Norme2(AA-BB); 667 // cout << " ll " << i << BB << ll << " " <<sensge[i] <<" on = " <<668 // Number(ge[i]) << " sens= " << sensge[i] ;669 630 AA=BB ;} 670 631 lge[tge]=ll+=Norme2(AA-V1); 671 // cout << " ll " << tge << " " << ll << sensge[tge]672 // <<" on = " << Number(ge[tge]) << " sens= " << sensge[tge] << endl;673 632 // search the geometrical edge 674 633 if (s>1.0){ … … 691 650 692 651 s=(ls-l0)/(l1-l0); 693 // cout << "on =" << Number(on) << sens0 << sens1 << "s0 " << s0 << " s1 ="694 // << s1 << " l0 =" << l0 << " ls= " << ls << " l1= " << l1 << " s= " << s;695 652 sg = s0 * (1.0-s) + s * s1; 696 653 } … … 699 656 } 700 657 V.r= on->F(sg); 701 // if (eg0 != eg1)702 // cout << "----- sg = "<< sg << " Sens =" << OppositeSens << " Edge = "703 // << Number(on) <<" V=" << V << endl;704 658 GV=VertexOnGeom(V,*on,sg); 705 659 return on; … … 710 664 long int verbosity=0; 711 665 712 if (verbosity>20)713 cout << "Geometry::AfterRead()" << nbv << " " << nbe << endl;714 666 Int4 i,k=0; ; 715 667 int jj; // jj in [0,1] … … 771 723 } 772 724 773 if(verbosity>7)774 for (i=0;i<nbv;i++)775 if (vertices[i].Required())776 cout << " The geo vertices " << i << " is required" << endl;777 778 725 for (i=0;i<nbv;i++) hv[i]=-1;// empty list 779 726 … … 785 732 } 786 733 eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ] 787 if(verbosity>9)788 cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<<endl;789 734 for (jj=0;jj<2;jj++) 790 735 { // generation of list … … 826 771 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; 827 772 float da12 = Abs(angle2-angle1); 828 if(verbosity>9) 829 cout <<" check angle " << i << " " << i1 << " " << i2 << " " << 180*(da12)/Pi 830 << " " << 180*MaximalAngleOfCorner/Pi << vertices[i] << endl; 831 832 if (( da12 >= MaximalAngleOfCorner ) 833 && (da12 <= 2*Pi -MaximalAngleOfCorner)) { 773 if (( da12 >= MaximalAngleOfCorner ) && (da12 <= 2*Pi -MaximalAngleOfCorner)) { 834 774 vertices[i].SetCorner() ; 835 if(verbosity>7) 836 cout << " The vertex " << i << " is a corner (angle) " 837 << 180*(da12)/ Pi<< " " << 180*MaximalAngleOfCorner/Pi << endl;} 838 // if the ref a changing then is SetRequired(); 839 840 if (edges[i1].flag != edges[i2].flag || edges[i1].Required()) 841 { 842 vertices[i].SetRequired(); 843 if(verbosity>7) 844 cout << " The vertex " << i << " is Required the flag change (crack or equi, or require)" << endl;} 845 846 if (edges[i1].ref != edges[i2].ref) { 847 vertices[i].SetRequired(); 848 if(verbosity>7) 849 cout << " The vertex " << i << " is Required ref" << endl;} 850 } ; 775 } 776 // if the ref a changing then is SetRequired(); 777 if (edges[i1].flag != edges[i2].flag || edges[i1].Required()){ 778 vertices[i].SetRequired(); 779 } 780 if (edges[i1].ref != edges[i2].ref) { 781 vertices[i].SetRequired(); 782 } 783 } 851 784 852 785 if(ord != 2) { 853 786 vertices[i].SetCorner(); 854 if(verbosity>7)855 cout << " the vertex " << i << " is a corner ordre = " << ord << endl;856 787 } 857 788 // close the liste around the vertex … … 875 806 edges[i1].Adj[j1] = edges + i; 876 807 edges[i1].SensAdj[j1] = jj; 877 if (verbosity>10)878 cout << " edges. Adj " << i1 << " " << j1 << " <--- " << i << " " << jj << endl;879 808 } 880 809 … … 894 823 ltg = Norme2(tg); 895 824 tg = tg *(lAB/ltg),ltg=lAB; 896 /*897 if(edges[i].ref >=4)898 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 = "899 << edges[i].v[1-jj]->r.y -edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r.y << endl;900 */901 825 } 902 826 … … 908 832 if ( (tg,AB) < 0) 909 833 tg = -tg; 910 //if(edges[i].ref >=4) cout << " tg = " << tg << endl;911 834 edges[i].tg[jj] = tg; 912 835 } // for (jj=0;jj<2;jj++) … … 915 838 if (ltg2[1]!=0) edges[i].SetTgB(); 916 839 } // for (i=0;i<nbe;i++) 917 918 if(verbosity>7)919 for (i=0;i<nbv;i++)920 if (vertices[i].Required())921 cout << " The geo vertices " << i << " is required " << endl;922 840 923 841 for (int step=0;step<2;step++) … … 959 877 960 878 }// for(;;) 961 if(verbosity>10 && curves==0) cout << NbOfCurves <<" curve : nb edges= "<< nee<< endl;962 879 NbOfCurves++; 963 880 if(level) { 964 if(verbosity>4)965 cout << " Warning: Curve "<< NbOfCurves << " without required vertex "966 << "so the vertex " << Number(a) << " become required " <<endl;967 881 a->SetRequired(); 968 882 } … … 999 913 } 1000 914 1001 if(verbosity>3)1002 cout << " End ReadGeometry: Number of curves in geometry is " << NbOfCurves <<endl;1003 if(verbosity>4)1004 for(int i=0;i<NbOfCurves ;i++)1005 {1006 cout << " Curve " << i << " begin e=" << Number(curves[i].be) << " k=" << curves[i].kb1007 << " end e= " << Number(curves[i].ee) << " k=" << curves[i].ke << endl;1008 }1009 915 delete []ev; 1010 916 delete []hv; -
issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
r2854 r2856 207 207 } 208 208 else 209 L[i]= l += s,S[i]=sss+=k,i++; //cout << i << " l = " << l << " sss = " << sss << endl;209 L[i]= l += s,S[i]=sss+=k,i++; 210 210 } 211 211 // warning for optimisation S is in [0:0.5] not in [0:1] … … 215 215 LastMetricInterpole.lab=l; 216 216 LastMetricInterpole.opt=i; 217 if (i>200 && kkk++<10) 218 cout << "Warning LengthInterpole: ( i = " << i << " l = " << l << " sss " << sss << " ) " << sstop <<endl; 217 if (i>200 && kkk++<10) printf("WARNING: LengthInterpole: ( i=%i l=%i sss=%i ) %g\n",i,l,sss,sstop); 219 218 return l; 220 219 } … … 243 242 i=k; // L[i]<l 244 243 }; 245 // cout << i << " " << j <<" " << L[i] << " " << L[j] << " " << S[i] << " " << S[j] << " l=" << l << endl;246 244 if (i==j) 247 245 r = 2*S[i]; -
issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
r2854 r2856 321 321 R2 XY(X,b->v[k]->r); 322 322 Real8 dd; 323 // old code if( Mx(XY) + b->v[k]->m(XY) < seuil ) 324 if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ) 325 { 326 // cout << CurrentTh->Number(v) << "is To Close " 327 // << CurrentTh->Number( b->v[k]) << " l=" <<dd<<endl; 323 if( (dd= LengthInterpole(Mx(XY), b->v[k]->m(XY))) < seuil ){ 328 324 return b->v[k]; 329 325 } … … 366 362 register long i=w.i.x, j=w.i.y,l=MaxISize; 367 363 pb = &root; 368 // cout << pb << " " << &root << endl;369 364 while( (b=*pb) && (b->n<0)) 370 365 { … … 399 394 if (!bb) 400 395 bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox 401 // cout << bb << " " << k << " " << ij << endl;402 396 bb->v[bb->n++] = v4[k]; 403 397 } … … 406 400 if (!(b = *pb)) 407 401 b=*pb= NewQuadTreeBox(); // alloc the QuadTreeBox 408 // cout << b << " " << b->n << endl;409 402 b->v[b->n++]=&w; // we add the vertex 410 403 NbVertices++; -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2854 r2856 79 79 } 80 80 k=0; 81 for (i=0;i<Tho.nbv;i++) 82 if (kk[i]>=0)83 kk[i]=k++;84 cout << " number of vertices " << k << " remove = " << Tho.nbv - k << endl;85 cout << " number of triangles " << kt << " remove = " << nbInT-kt << endl;86 cout << " number of New boundary edge " << nbNewBedge << endl;81 for (i=0;i<Tho.nbv;i++){ 82 if (kk[i]>=0) kk[i]=k++; 83 } 84 printf(" number of vertices %i, remove = %i\n",k,Tho.nbv - k); 85 printf(" number of triangles %i, remove = \n",kt,nbInT-kt); 86 printf(" number of New boundary edge %i\n",nbNewBedge); 87 87 Int4 inbvx =k; 88 88 PreInit(inbvx,cname); … … 465 465 /*Recreate geometry if needed*/ 466 466 if(1!=0) { 467 printf("W arning: mesh present but no geometry found. Reconstructing...\n");467 printf("WARNING: mesh present but no geometry found. Reconstructing...\n"); 468 468 /*Recreate geometry: */ 469 469 ConsGeometry(bamgopts->MaximalAngleOfCorner*Pi/180); … … 767 767 768 768 Int4 nbeold = nbe; 769 for (i=0;i<nbe;i++) 770 { 771 // cout << i << " " << Number(edges[i][0]) << " " << Number(edges[i][1]) << endl; 769 for (i=0;i<nbe;i++){ 772 770 edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])); 773 771 } 774 772 if (nbe != edge4->nb()){ 775 773 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i, nbe4=%i",nbe,edge4->nb())); … … 822 820 j = (int) ((-2-st[i])%3); 823 821 Triangle & tt = * triangles[it].TriangleAdj(j); 824 //cout << it << " c=" << triangles[it].color << " " << Number(tt) << " c=" << tt.color << endl; 825 if (triangles[it].color != tt.color|| i < nbeold) // Modif FH 06122055 // between 2 sub domai 826 k++; 822 if (triangles[it].color != tt.color|| i < nbeold) k++; 827 823 } 828 824 else if (st[i] >=0) // edge alone 829 // if (i >= nbeold)830 825 kk++; 831 826 832 if(verbosity>4 && (k+kk) )833 cout << " Nb of ref edge " << kk+k << " (internal " << k << ")"834 << " in file " << nbe << endl;835 827 k += kk; 836 828 kk=0; … … 842 834 k =0; 843 835 // construction of the edges 844 if(verbosity>4) 845 cout << " Construction of the edges " << nbe << endl; 846 847 for (i=0;i<nbedges;i++) 848 { 836 if(verbosity>4) printf(" Construction of the edges %i\n",nbe); 837 838 for (i=0;i<nbedges;i++){ 849 839 Int4 add= -1; 850 840 … … 952 942 else 953 943 level-=2; 954 if (verbosity>5)955 cout << " Nb of triangles " << k << " of Subdomain "956 << NbSubDomains << " " << kolor << endl;957 944 NbSubDomains++; 958 945 } 959 946 } 960 if (verbosity> 3) 961 cout << " The Number of sub domain = " << NbSubDomains << endl; 947 if (verbosity> 3) printf(" The Number of sub domain = %i\n",NbSubDomains); 962 948 963 949 Int4 isd; … … 1009 995 Gh.NbSubDomains = NbSubDomains; 1010 996 Gh.subdomains = new GeometricalSubDomain[NbSubDomains]; 1011 if (verbosity>3) 1012 cout << " Nb of vertices = " << Gh.nbv << " Nb of edges = " << Gh.nbe << endl; 997 if (verbosity>3) printf(" number of vertices = %i\n number of edges = %i\n",Gh.nbv,Gh.nbe); 1013 998 NbVerticesOnGeomVertex = Gh.nbv; 1014 999 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; … … 1074 1059 int sens=equiedges[i]%2; 1075 1060 if(i!=j && equiedges[i]>=0) { 1076 if(verbosity>9)1077 cout << " Edges Equi " << i << " <=> " << j << " sens = " << sens << endl;1078 1061 if( sens==0) 1079 1062 Gh.edges[i].SetEqui(); … … 1200 1183 e = vB.onbe->be; 1201 1184 1202 // cout << " EXCHANGE A et B) " << endl;1203 1185 } 1204 1186 else{ // do the search by walking … … 1226 1208 Real8 te0; 1227 1209 // we suppose take the curve's abcisse 1228 // cout << kkk << " e = " << BTh.Number(e) << " v0= "1229 // << BTh.Number(v0) << " v1 = " << BTh.Number((*e)[sens]) << endl;1230 1210 for ( eee=e,iii=sens,te0=tA; 1231 1211 eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; 1232 1212 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { 1233 // cout << kkk << " eee = " << BTh.Number(eee) << " v0= "1234 // << BTh.Number(v0) << " v1 = " << BTh.Number(v1) << endl;1235 1213 1236 1214 kkk=kkk+1; … … 1289 1267 long int verbosity=0; 1290 1268 1291 if (verbosity>2) 1292 cout << " -- MakeQuadrangles costheta = " << costheta << endl; 1293 if (verbosity>5) 1294 cout << " (in) Nb of Quadrilaterals = " << NbOfQuad 1295 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 1296 << " Nb of outside triangles = " << NbOutT << endl; 1269 if (verbosity>2) printf("MakeQuadrangles costheta = %g\n",costheta); 1297 1270 1298 1271 if (costheta >1) { 1299 if (verbosity>5) 1300 cout << " do nothing costheta >1 "<< endl; 1301 return;} 1272 if (verbosity>5) printf(" do nothing: costheta > 1\n"); 1273 } 1302 1274 1303 1275 … … 1316 1288 1317 1289 Int4 kk=0; 1318 for (ij=0;ij<k;ij++) 1319 { 1320 // cout << qq[ij].q << " " << endl; 1290 for (ij=0;ij<k;ij++) { 1321 1291 i=qq[ij].i3j/3; 1322 1292 j=(int) (qq[ij].i3j%3); … … 1326 1296 } 1327 1297 NbOfQuad = kk; 1328 if (verbosity>2) 1329 { 1330 cout << " (out) Nb of Quadrilaterals = " << NbOfQuad 1331 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 1332 << " Nb of outside triangles = " << NbOutT << endl; 1333 } 1298 if (verbosity>2){ 1299 printf(" number of quadrilaterals = %i\n",NbOfQuad); 1300 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); 1301 printf(" number of outside triangles = %i\n",NbOutT); 1302 } 1334 1303 delete [] qq; 1335 1304 } … … 1341 1310 Direction NoDirOfSearch; 1342 1311 const int withBackground = &BTh != this && &BTh; 1343 if (verbosity>2)1344 cout << " -- SplitElement " << (choice? " Q->4Q and T->4T " : " Q->4Q or T->3Q " ) << endl;;1345 if (verbosity>5)1346 cout << endl << " (in) Nb of Quadrilaterals = " << NbOfQuad1347 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*21348 << " Nb of outside triangles = " << NbOutT << endl;1349 1312 1350 1313 ReNumberingTheTriangleBySubDomain(); … … 1417 1380 newedges[ie].adj[1]=newedges + ie +1; 1418 1381 R2 A = edges[i][0],B = edges[i][1]; 1419 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl;1420 1382 1421 1383 … … 1431 1393 throw ErrorException(__FUNCT__,exprintf("!edgesGtoB")); 1432 1394 } 1433 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl;1434 1395 ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].on)], 1435 1396 edges[i][0],edges[i][1],0.5,vertices[k], … … 1446 1407 vertices[k].m = Metric(1-s,bv0,s,bv1); 1447 1408 kvb++; 1448 // cout << " ie = " << ie <<" v0 = " << Number(newedges[ie][0]) << endl;1449 1409 } 1450 1410 else … … 1479 1439 newedges[ie].on = Gh.Contening(BB,ong); 1480 1440 newedges[ie++].v[0]=vertices+k; 1481 // cout << " ie = " << ie-2 << " vm " << k << " v0 = " << Number(newedges[ie-2][0])1482 // << " v1 = " << Number(newedges[ie-1][1])1483 // << " ong =" << ong-Gh.edges1484 // << " on 0 =" << newedges[ie-2].on -Gh.edges << AA1485 // << " on 1 =" << newedges[ie-1].on -Gh.edges << BB1486 // << endl;1487 1441 k++; 1488 1442 } … … 1554 1508 ksplit[i]=1; // no split by default 1555 1509 const Triangle & t = triangles[ i]; 1556 // cout << " Triangle " << i << " " << t << !!t.link << ":: " ;1557 1510 int nbsplitedge =0; 1558 1511 int nbinvisible =0; … … 1569 1522 const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]]; 1570 1523 const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]]; 1571 // cout << " ke = " << kedge[3*i+j] << " " << Number(v0) << " " << Number(v1) << "/ ";1572 1524 if ( kedge[3*i+j] < 0) 1573 1525 { 1574 1526 Int4 ke =edge4->findtrie(Number(v0),Number(v1)); 1575 // cout << ":" << ke << "," << !!t.link << " " << &tt ;1576 1527 if (ke<0) // new 1577 1528 { … … 1611 1562 throw ErrorException(__FUNCT__,exprintf("nbinvisible>=2")); 1612 1563 } 1613 // cout << " " << nbinvisible << " " << nbsplitedge << endl;1614 1564 switch (nbsplitedge) { 1615 1565 case 0: ksplit[i]=10; newnbt++; break; // nosplit … … 1628 1578 newNbOfQuad = 4*NbOfQuad; 1629 1579 nbv = k; 1630 // cout << " Nbv = " << nbv << endl;1631 1580 kkk = nbt; 1632 1581 ksplit[-1] = nbt; 1633 1582 // look on old true triangles 1634 1583 1635 for (i=0;i<nbtsave;i++) 1636 { 1637 // cout << "Triangle " << i << " " << ksplit[i] << ":" << triangles[i] 1638 // << " ----------------------------------------------- " <<endl; 1639 // Triangle * tc=0; 1584 for (i=0;i<nbtsave;i++){ 1640 1585 int nbmkadj=0; 1641 1586 Int4 mkadj [100]; … … 1727 1672 Vertex * v02 = vertices + kedge[3*i+k1]; 1728 1673 Vertex * v01 = vertices + kedge[3*i+k2]; 1729 // cout << Number(t0(i0)) << " " << Number(t0(i1))1730 // << " " << Number(t0(i2))1731 // << " " << kedge[3*i+k0]1732 // << " " << kedge[3*i+k1]1733 // << " " << kedge[3*i+k2] << endl;1734 1674 t0(i1) = v01; 1735 1675 t0(i2) = v02; … … 1790 1730 } 1791 1731 1792 // cout << " -- " << i << " " << nbmkadj << " " << kkk << " " << tc << endl;1793 // t0.SetDetf();1794 1732 // save all the new triangles 1795 1733 mkadj[nbmkadj++]=i; … … 1801 1739 t0.link= triangles+jj; 1802 1740 mkadj[nbmkadj++]=jj; 1803 // triangles[jj].SetDet();1804 1741 } 1805 // cout << " -- " << i << " " << nbmkadj << endl;1806 1742 if (nbmkadj>13){// 13 = 6 + 4 + 1807 1743 throw ErrorException(__FUNCT__,exprintf("nbmkadj>13")); … … 1809 1745 1810 1746 if (kk==6) newNbOfQuad+=3; 1811 // triangles[i].Draw(); 1812 1813 for (jj=ksplit[i-1];jj<kkk;jj++) 1814 // triangles[jj].SetDet(); 1815 // triangles[jj].Draw(); 1816 1817 1818 1819 nbt = kkk; 1747 for (jj=ksplit[i-1];jj<kkk;jj++) nbt = kkk; 1820 1748 ksplit[i]= nbt; // save last adresse of the new triangles 1821 1749 kkk = nbt; 1822 1823 1750 } 1824 1751 1825 // cout << " nv = " << nbv << " nbt = " << nbt << endl; 1826 for (i=0;i<nbv;i++) 1827 vertices[i].m = vertices[i].m*2.; 1828 // 1752 for (i=0;i<nbv;i++) vertices[i].m = vertices[i].m*2.; 1753 1829 1754 if(withBackground) 1830 1755 for (i=0;i<BTh.nbv;i++) … … 1867 1792 FillHoleInMesh(); 1868 1793 1869 if (verbosity>2) 1870 cout << " (out) Nb of Quadrilaterals = " << NbOfQuad 1871 << " Nb Of Triangles = " << nbt-NbOutT- NbOfQuad*2 1872 << " Nb of outside triangles = " << NbOutT << endl; 1794 if (verbosity>2){ 1795 printf(" number of quadrilaterals = %i\n",NbOfQuad); 1796 printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); 1797 printf(" number of outside triangles = %i\n",NbOutT); 1798 } 1873 1799 1874 1800 CurrentTh=OCurrentTh; … … 2139 2065 // swap if the point s is on a edge 2140 2066 if(izerodet>=0) { 2141 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl;2142 2067 int rswap =tt[izerodet]->swap(iedge); 2143 2068 … … 2165 2090 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 2166 2091 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; 2167 if (v0.on && v1.on) 2168 { 2092 if (v0.on && v1.on){ 2169 2093 R2 P= ((R2) v0 + (R2) v1)*0.5; 2170 2094 if ( nbv<nbvx) { … … 2175 2099 } 2176 2100 NbSplitEdge++; 2177 if (verbosity>7)2178 cout <<" Internal edge with two vertices on boundary"2179 << Number(v0) << " " << Number(v1) << " by " << endl;2180 2101 } 2181 2102 } … … 2191 2112 vi.i = toI2(vi.r); 2192 2113 vi.r = toR2(vi.i); 2193 // if (!quadtree->ToClose(vi,seuil,hi,hj)) { 2114 2194 2115 // a good new point 2195 2116 vi.ReferenceNumber=0; 2196 2117 vi.DirOfSearch =NoDirOfSearch; 2197 // cout << " Add " << Number(vi) << " " << vi2198 // << " " << Number(vi) << " <--> " << Number(vi) <<endl;2199 2118 Triangle *tcvi = FindTriangleContening(vi.i,dete); 2200 2119 if (tcvi && !tcvi->link) { … … 2231 2150 2232 2151 const Int4 nbvnew = nbv-nbvold; 2233 if (verbosity>5) 2234 cout << " Try to Insert the " <<nbvnew<< " new points " << endl; 2152 if (verbosity>5) printf(" Try to Insert the %i new points\n",nbvnew); 2235 2153 Int4 NbSwap=0; 2236 2154 Icoor2 dete[3];
Note:
See TracChangeset
for help on using the changeset viewer.
