![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2850
- Timestamp:
- 01/15/10 09:10:43 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 11 edited
-
Mesh2.h (modified) (1 diff)
-
objects/GeometricalEdge.cpp (modified) (1 diff)
-
objects/Geometry.cpp (modified) (8 diffs)
-
objects/ListofIntersectionTriangles.cpp (modified) (2 diffs)
-
objects/MatVVP2x2.cpp (modified) (1 diff)
-
objects/MetricAnIso.cpp (modified) (1 diff)
-
objects/QuadTree.cpp (modified) (1 diff)
-
objects/SetOfE4.cpp (modified) (3 diffs)
-
objects/Triangle.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (48 diffs)
-
objects/Vertex.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2814 r2850 791 791 catch(...) { this->~Triangles(); throw; } } 792 792 Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR 793 // Triangles(Triangles &){ std::cerr << " BUG call copy opretor of Triangles" << std::endl;MeshError(111);}794 793 Triangles(const Triangles &,const int *flag,const int *bb); // truncature 795 794 -
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
r2841 r2850 8 8 #include "../QuadTree.h" 9 9 #include "../SetOfE4.h" 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "GeometricalEdge" 10 13 11 14 using namespace std; -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2841 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "Geometry" 13 16 14 17 namespace bamg { … … 540 543 const Vertex &v0=e[0],&v1=e[1]; 541 544 V.m = Metric(1.0-s, v0,s, v1); 542 #define MXE__LINE __LINE__+1543 545 const int mxe =100; 544 546 GeometricalEdge *ge[mxe+1]; … … 566 568 << " sg 1= " << vg1 567 569 << "--------------------------------------------" << endl; 568 while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0) 569 { 570 while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0){ 570 571 if (bge<=0) { 571 572 // int kkk; 572 // if (NbTry) cout <<"Read (int) to Show Sioux window", cin >> kkk ; 573 if(NbTry) 574 { 575 cerr << " -- Fatal Error: on the class triangles before call Geometry::ProjectOnCurve" << endl; 576 cerr << " The mesh of the Geometry is to fine: "; 577 cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl; 578 cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl; 579 cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl; 580 cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl; 581 MeshError(222); 573 if(NbTry) { 574 printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); 575 printf("That bug might come from:\n"); 576 printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); 577 printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); 578 printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); 579 throw ErrorException(__FUNCT__,exprintf("see above")); 582 580 } 583 581 NbTry++; … … 603 601 NbTry++; 604 602 if (NbTry<2) goto retry; 605 cerr << " The mesh of the Geometry is to fine:";606 cerr << " 1) a mesh edge contening more than "<< mxe/2 << " geometrical edges." << endl;607 cerr << " 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before " << endl;608 cerr << " To solve the problem do a coarsening of the geometrical mesh " << endl;609 cerr << " or change the constant value of mxe in " << __FILE__ << " line " << MXE__LINE << "( dangerous way )" << endl;610 MeshError(223);603 printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); 604 printf("That bug might come from:\n"); 605 printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); 606 printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); 607 printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); 608 throw ErrorException(__FUNCT__,exprintf("see above")); 611 609 } 612 610 … … 702 700 Int4 * ev = new Int4 [ 2 * nbe ]; 703 701 float * eangle = new float[ nbe ]; 702 double eps = 1e-20; 703 QuadTree quadtree; // to find same vertices 704 Vertex * v0 = vertices; 705 GeometricalVertex * v0g = (GeometricalVertex *) (void *) v0; 706 707 for (i=0;i<nbv;i++) 708 vertices[i].link = vertices +i; 709 for (i=0;i<nbv;i++) 704 710 { 705 double eps = 1e-20; 706 QuadTree quadtree; // to find same vertices 707 Vertex * v0 = vertices; 708 GeometricalVertex * v0g = (GeometricalVertex *) (void *) v0; 709 int k=0; 710 for (i=0;i<nbv;i++) 711 vertices[i].link = vertices +i; 712 for (i=0;i<nbv;i++) 713 { 714 vertices[i].i = toI2(vertices[i].r); // set integer coordinate 715 Vertex *v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 716 if( v && Norme1(v->r - vertices[i]) < eps ) 717 { // link v & vertices[i] 718 // vieille ruse pour recuperer j 719 GeometricalVertex * vg = (GeometricalVertex *) (void *) v; 720 int j = vg-v0g; 721 assert( v == & (Vertex &) vertices[j]); 722 vertices[i].link = vertices + j; 723 k++; 711 vertices[i].i = toI2(vertices[i].r); // set integer coordinate 712 Vertex *v= quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 713 if( v && Norme1(v->r - vertices[i]) < eps ) 714 { // link v & vertices[i] 715 // vieille ruse pour recuperer j 716 GeometricalVertex * vg = (GeometricalVertex *) (void *) v; 717 int j = vg-v0g; 718 assert( v == & (Vertex &) vertices[j]); 719 vertices[i].link = vertices + j; 720 k++; 721 } 722 else quadtree.Add(vertices[i]); 723 } 724 if (k) { 725 printf("number of distinct vertices= %i, over %i\n",nbv - k,nbv); 726 printf("List of duplicate vertices:\n"); 727 for (i=0;i<nbv;i++){ 728 if (!vertices[i].IsThe()) printf(" %i and %i\n",i,Number(vertices[i].The())); 729 } 730 throw ErrorException(__FUNCT__,exprintf("See above")); 731 } 732 733 // verification of cracked edge 734 for (i=0;i<nbe;i++){ 735 if (edges[i].Cracked() ) { 736 // verification of crack 737 GeometricalEdge & e1=edges[i]; 738 GeometricalEdge & e2=*e1.link; 739 cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " << e1[1].The() << " " << e2[1].The() << endl; 740 if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() ) 741 { 724 742 } 725 else quadtree.Add(vertices[i]); 726 } 727 if (k) { 728 cout << " Number of distinte vertices " << nbv - k << " Over " << nbv << endl; 729 //if (verbosity>10) 730 { 731 cout << " The duplicate vertex " << endl; 732 for (i=0;i<nbv;i++) 733 if (!vertices[i].IsThe()) 734 cout << " " << i << " and " << Number(vertices[i].The()) << endl; 735 MeshError(102); 736 //throw(ErrorExec("exit",1)); 737 } 738 } 739 740 // verification of cracked edge 741 int err =0; 742 for (i=0;i<nbe;i++) 743 if (edges[i].Cracked() ) 744 { 745 // verification of crack 746 GeometricalEdge & e1=edges[i]; 747 GeometricalEdge & e2=*e1.link; 748 cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " << e1[1].The() << " " << e2[1].The() << endl; 749 if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() ) 743 else 744 if ( e1[0].The() == e2[1].The() && e1[1].The() == e2[0].The() ) 750 745 { 746 //nothing 751 747 } 752 else 753 if ( e1[0].The() == e2[1].The() && e1[1].The() == e2[0].The() ) 754 { 755 } 756 else 757 { 758 err++; 759 cerr << " Cracked edges with no same vertex " << &e1-edges << " " << &e2 -edges << endl; 760 } 761 } 762 else 763 { 764 // if (!edges[i][0].IsThe()) err++; 765 // if (!edges[i][1].IsThe()) err++; 766 } 767 if (err) 768 { 769 cerr << " Some vertex was not distint and not on cracked edge " << err<< endl; 770 MeshError(222); 771 } 772 } 748 else { 749 throw ErrorException(__FUNCT__,exprintf("Cracked edges with no same vertex")); 750 } 751 } 752 } 753 773 754 if(verbosity>7) 774 755 for (i=0;i<nbv;i++) … … 776 757 cout << " The geo vertices " << i << " is required" << endl; 777 758 778 for (i=0;i<nbv;i++) 779 hv[i]=-1;// empty list 780 781 for (i=0;i<nbe;i++) 782 { 759 for (i=0;i<nbv;i++) hv[i]=-1;// empty list 760 761 for (i=0;i<nbe;i++) { 783 762 R2 v10 = edges[i].v[1]->r - edges[i].v[0]->r; 784 763 Real8 lv10 = Norme2(v10); 785 764 if(lv10 == 0) { 786 cerr << "The length of " <<i<< "th Egde is 0 " << endl;787 MeshError(1);}788 789 790 791 792 793 794 795 796 797 765 throw ErrorException(__FUNCT__,exprintf("Length of edge %i is 0",i)); 766 } 767 eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ] 768 if(verbosity>9) 769 cout << " angle edge " << i <<" " << eangle[i]*180/Pi<< v10<<endl; 770 for (jj=0;jj<2;jj++) 771 { // generation of list 772 Int4 v = Number(edges[i].v[jj]); 773 ev[k] = hv[v]; 774 hv[v] = k++; 775 } 776 } 798 777 // bulle sort on the angle of edge 799 778 for (i=0;i<nbv;i++) { … … 820 799 } // end while (exch) 821 800 822 if (ord >= 1 )823 { /*824 Int4 n = hv[i];825 while ( n >=0)826 { Int4 i1 = n/2,j1 = n%2;827 //float a = 180*(j1 ? OppositeAngle(eangle[i1]): eangle[i1])/Pi;828 n = ev[n];829 }830 */831 }832 801 if(ord == 2) { // angulare test to find a corner 833 802 Int4 n1 = hv[i]; … … 882 851 Int4 n1 = ev[k++]; 883 852 Int4 i1 = n1/2 ,j1=n1%2; 884 if( edges[i1].v[j1] != edges[i].v[jj]) 885 { cerr << " Bug Adj edge " << i << " " << jj << 886 " et " << i1 << " " << j1 << " k=" << k; 887 cerr << Number(edges[i].v[jj]) <<" <> " 888 << Number(edges[i1].v[j1]) <<endl; 889 cerr << "edge " << Number(edges[i].v[0]) << " " 890 << Number(edges[i].v[1]) << endl; 891 // cerr << "in file " <<filename <<endl; 892 MeshError(1); 893 } 853 if( edges[i1].v[j1] != edges[i].v[jj]) { 854 throw ErrorException(__FUNCT__,exprintf("Bug Adj edge")); 855 } 894 856 edges[i1].Adj[j1] = edges + i; 895 857 edges[i1].SensAdj[j1] = jj; -
issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
r2843 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "ListofIntersectionTriangles" 13 16 14 17 namespace bamg { … … 181 184 j= VerticesOfTriangularEdge[ocut][1]; 182 185 } 183 else { 184 cerr << " Bug Split Edge " << endl; 185 cerr << " dt[0]= " << dt[0] 186 << " dt[1]= " << dt[1] 187 << " dt[2]= "<< dt[2] << endl; 188 cerr << i0 << " " << i1 << " " << i2 << endl; 189 cerr << " A = " << A << " B= " << B << endl; 190 cerr << " Triangle t = " << *t << endl; 191 cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl; 192 cerr << " nbt = " << nbt << endl; 193 MeshError(100);}} 186 else { 187 throw ErrorException(__FUNCT__,exprintf("Bug Split Edge")); 188 } 189 } 194 190 195 191 k = OppositeVertex[ocut]; -
issm/trunk/src/c/Bamgx/objects/MatVVP2x2.cpp
r2841 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "MatVVP2x2" 13 16 14 17 namespace bamg { -
issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
r2841 r2850 8 8 #include "../QuadTree.h" 9 9 #include "../SetOfE4.h" 10 11 #undef __FUNCT__ 12 #define __FUNCT__ "MetricAnIso" 10 13 11 14 using namespace std; -
issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
r2841 r2850 5 5 #include "../Mesh2.h" 6 6 #include "../QuadTree.h" 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "QuadTree" 7 10 8 11 namespace bamg { -
issm/trunk/src/c/Bamgx/objects/SetOfE4.cpp
r2842 r2850 1 1 #include <iostream> 2 3 #include "../../shared/shared.h" 4 #include "../../include/macros.h" 5 #include "../../toolkits/toolkits.h" 2 6 #include "../meshtype.h" 3 7 #include "../SetOfE4.h" 8 9 #undef __FUNCT__ 10 #define __FUNCT__ "SetOfEdges4" 4 11 5 12 using namespace std; … … 24 31 Int4 SetOfEdges4::find(Int4 ii,Int4 jj) { 25 32 if (tete == 0 ) { 26 cerr <<"SetOfEdges4::find \nplus de tete de liste\n";27 MeshError(888);}28 33 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::find no more tete de liste (?)")); 34 } 35 Int4 n = tete[ Abs( ii ) % nx ]; 29 36 30 while (n >= 0) 31 if (ii == Edges[n].i && jj == Edges[n].j) 32 return n; 33 else n = Edges[n].next; 34 return -1; // n'existe pas 37 while (n >= 0) 38 if (ii == Edges[n].i && jj == Edges[n].j) return n; 39 else n = Edges[n].next; 40 return -1; //do not exist, return -1 35 41 } 36 42 /*}}}1*/ … … 38 44 Int4 SetOfEdges4::add(Int4 ii,Int4 jj) { 39 45 if (tete == 0 ) { 40 cerr << "SetOfEdges4::add\n plus de tete de liste \n" << endl; 41 MeshError(888);} 46 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add no more tete de liste (?)")); 47 } 48 Int4 h; 49 Int4 n = tete[ h = Abs( ii ) % nx ]; 50 while (n >= 0) 51 if (ii == Edges[n].i && jj == Edges[n].j) 52 return n; 53 else n = Edges[n].next; 54 if (nbax <=NbOfEdges ) { 55 throw ErrorException(__FUNCT__,exprintf("SetOfEdges4::add overflow: NbOfEdges=%i > nbax=%i",NbOfEdges,nbax)); 56 } 42 57 43 Int4 h; 44 Int4 n = tete[ h = Abs( ii ) % nx ]; 45 while (n >= 0) 46 if (ii == Edges[n].i && jj == Edges[n].j) 47 return n; 48 else n = Edges[n].next; 49 if (nbax <=NbOfEdges ) { 50 cerr << " SetOfEdges4::add\noverflow de la pile " << nbax << " " << NbOfEdges << endl; 51 MeshError(888);} 52 53 Edges[NbOfEdges].i=ii; 54 Edges[NbOfEdges].j=jj; 55 Edges[NbOfEdges].next= tete[h]; 56 tete[h] = NbOfEdges; 57 return NbOfEdges ++; 58 Edges[NbOfEdges].i=ii; 59 Edges[NbOfEdges].j=jj; 60 Edges[NbOfEdges].next= tete[h]; 61 tete[h] = NbOfEdges; 62 return NbOfEdges ++; 58 63 } 59 64 /*}}}1*/ -
issm/trunk/src/c/Bamgx/objects/Triangle.cpp
r2841 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "Triangle" 13 16 14 17 namespace bamg { -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2849 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "Triangles" 13 16 14 17 namespace bamg { … … 110 113 assert(kt==nbt); 111 114 if (nbt ==0 && nbv ==0) { 112 cout << "Error all triangles was remove " << endl; 113 MeshError(999,this); 115 throw ErrorException(__FUNCT__,exprintf("All triangles have been removed")); 114 116 } 115 117 delete [] kk; … … 733 735 // -------------------------- 734 736 long int verbosity=0; 735 if (verbosity>1) 736 cout << " -- construction of the geometry from the 2d mesh " << endl; 737 if (nbt<=0 || nbv <=0 ) { MeshError(101);} 737 if (verbosity>1) printf(" construction of the geometry from the 2d mesh\n"); 738 if (nbt<=0 || nbv <=0 ) { 739 throw ErrorException(__FUNCT__,exprintf("nbt or nbv is negative")); 740 } 738 741 739 742 // construction of the edges … … 765 768 edge4->addtrie(Number(edges[i][0]),Number(edges[i][1])); 766 769 } 767 if (nbe != edge4->nb()) 768 { 769 cerr << " Some Double edge in the mesh, the number is " << nbe 770 << " nbe4=" << edge4->nb() << endl; 771 MeshError(1002); 770 if (nbe != edge4->nb()){ 771 throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i, nbe4=%i",nbe,edge4->nb())); 772 772 } 773 for (i=0;i<nbt;i++) 774 for (j=0;j<3;j++) 775 { 776 // Int4 i0,i1; 777 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), 778 Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 779 Int4 invisible = triangles[i].Hidden(j); 780 if(st[k]==-1) 781 st[k]=3*i+j; 782 else if(st[k]>=0) { 783 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); 784 785 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 786 if (invisible) triangles[i].SetHidden(j); 787 if (k<nbe) { 788 triangles[i].SetLocked(j); 789 } 790 st[k]=-2-st[k]; } 791 else { 792 cerr << " The edge (" 793 << Number(triangles[i][VerticesOfTriangularEdge[j][0]]) 794 << " , " 795 << Number(triangles[i][VerticesOfTriangularEdge[j][1]]) 796 << " ) is in more than 2 triangles " <<k <<endl; 797 cerr << " Edge " << j << " Of Triangle " << i << endl; 798 cerr << " Edge " << (-st[k]+2)%3 << " Of Triangle " << (-st[k]+2)/3 << endl; 799 cerr << " Edge " << triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)) 800 << " Of Triangle " << Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3))) 801 << endl; 802 MeshError(9999);} 803 804 773 for (i=0;i<nbt;i++){ 774 for (j=0;j<3;j++) { 775 Int4 k =edge4->addtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 776 Int4 invisible = triangles[i].Hidden(j); 777 if(st[k]==-1) st[k]=3*i+j; 778 else if(st[k]>=0) { 779 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); 780 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 781 if (invisible) triangles[i].SetHidden(j); 782 if (k<nbe) { 783 triangles[i].SetLocked(j); 784 } 785 st[k]=-2-st[k]; 786 } 787 else { 788 printf("The edge (%i,%i) belongs to more than 2 triangles (%i)\n",Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]),k); 789 printf("Edge %i of triangle %i\n",j,i); 790 printf("Edge %i of triangle %i\n",(-st[k]+2)%3,(-st[k]+2)/3); 791 printf("Edge %i of triangle %i\n",triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((int)((-st[k]+2)%3)),Number(triangles[(-st[k]+2)/3].TriangleAdj((int)((-st[k]+2)%3)))); 792 throw ErrorException(__FUNCT__,exprintf("An edge belongs to more than 2 triangles")); 793 } 805 794 } 795 } 806 796 Int4 nbedges = edge4->nb(); // the total number of edges 807 797 delete edge4; … … 809 799 810 800 if(verbosity>5) { 811 if (name) 812 cout << " On Mesh " << name << endl; 813 cout << " - The number of Vertices = " << nbv << endl; 814 cout << " - The number of Triangles = " << nbt << endl; 815 cout << " - The number of given edge = " << nbe << endl; 816 cout << " - The number of all edges = " << nbedges << endl; 817 cout << " - The Euler number = 1-Nb Of Hole = " << nbt-nbedges+nbv << endl; } 818 819 801 printf(" info on Mesh %s:\n",name); 802 printf(" - number of vertices = %i \n",nbv); 803 printf(" - number of triangles = %i \n",nbt); 804 printf(" - number of given edges = %i \n",nbe); 805 printf(" - number of all edges = %i \n" ,edge4->nb()); 806 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 807 } 820 808 // check the consistant of edge[].adj and the geometrical required vertex 821 809 k=0; … … 824 812 825 813 for (i=0;i<nbedges;i++) 826 if (st[i] <-1) // edge internal 827 { 814 if (st[i] <-1) {// edge internal 828 815 it = (-2-st[i])/3; 829 816 j = (int) ((-2-st[i])%3); … … 832 819 if (triangles[it].color != tt.color|| i < nbeold) // Modif FH 06122055 // between 2 sub domai 833 820 k++; 834 821 } 835 822 else if (st[i] >=0) // edge alone 836 823 // if (i >= nbeold) … … 842 829 k += kk; 843 830 kk=0; 844 if (k) 845 { 846 847 // if (nbe) { 848 // cerr << k << " boundary edges are not defined as edges " << endl; 849 // MeshError(9998); 850 // } 831 if (k) { 851 832 // construction of the edges 852 833 nbe = k; … … 1138 1119 } 1139 1120 else 1140 MeshError(103);1121 throw ErrorException(__FUNCT__,exprintf("%i should be >=0")); 1141 1122 } 1142 1123 … … 1157 1138 R2 A=vA,B=vB; 1158 1139 Vertex * pvA=&vA, * pvB=&vB; 1159 if (vA.vint == IsVertexOnVertex) 1160 { 1161 // cout << " Debut vertex = " << BTh.Number(vA.onbv) ; 1140 if (vA.vint == IsVertexOnVertex){ 1162 1141 pA=vA.onbv; 1163 } 1164 else if (vA.vint == IsVertexOnEdge) 1165 { 1142 } 1143 else if (vA.vint == IsVertexOnEdge){ 1166 1144 pA=vA.onbe->be; 1167 1145 tA=vA.onbe->abcisse; 1168 // cout << " Debut edge = " << BTh.Number(vA.onbv) << " " << tA ; 1169 1170 } 1171 else 1172 {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vA) << " " << endl; 1173 cerr << " forget call to SetVertexFieldOnBTh" << endl; 1174 MeshError(-1); 1175 } 1176 1177 if (vB.vint == IsVertexOnVertex) 1178 { 1179 // cout << " Fin vertex = " << BTh.Number(vB.onbv) << endl; 1146 } 1147 else { 1148 throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA))); 1149 } 1150 1151 if (vB.vint == IsVertexOnVertex){ 1180 1152 pB=vB.onbv; 1181 } 1182 else if(vB.vint == IsVertexOnEdge) 1183 { 1153 } 1154 else if(vB.vint == IsVertexOnEdge){ 1184 1155 pB=vB.onbe->be; 1185 1156 tB=vB.onbe->abcisse; 1186 // cout << " Fin edge = " << BTh.Number(vB.onbe->be) << " " << tB ; 1187 1188 } 1189 else 1190 {cerr << "ProjectOnCurve On Vertex " << BTh.Number(vB) << " " << endl; 1191 cerr << " forget call to SetVertexFieldOnBTh" << endl; 1192 MeshError(-1); 1193 } 1157 } 1158 else { 1159 throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB))); 1160 } 1194 1161 Edge * e = &BhAB; 1195 1162 assert( pA && pB && e); … … 1233 1200 Real8 abscisse = -1; 1234 1201 1235 for (int cas=0;cas<2;cas++) 1236 {// 2 times algo:1202 for (int cas=0;cas<2;cas++){ 1203 // 2 times algo: 1237 1204 // 1 for computing the length l 1238 1205 // 2 for find the vertex … … 1289 1256 assert(thetab>=0 && thetab<=1); 1290 1257 BR = VertexOnEdge(&R,eee,thetab); 1291 1292 // cout << kkk << " eee = " << BTh.Number(eee) << " v0= "1293 // << BTh.Number(v0) << " " << te01294 // << " v1 = " << BTh.Number(v1) << " " << tB << endl;1295 1296 //out << Number(R) << " Opt = " << thetab << " on " << BTh.Number(eee)1297 // << " = " << R << endl;1298 1299 1258 return Gh.ProjectOnCurve(*eee,thetab,R,GR); 1300 1259 } … … 1303 1262 1304 1263 } 1305 cerr << " Big Bug" << endl; 1306 MeshError(678); 1264 throw ErrorException(__FUNCT__,exprintf("Big bug...")); 1307 1265 return 0; // just for the compiler 1308 1309 1266 } 1310 1267 /*}}}1*/ … … 1529 1486 // 1530 1487 1531 for (i=0;i<nbt;i++) 1532 { 1533 1488 for (i=0;i<nbt;i++) { 1534 1489 Triangle & t = triangles[i]; 1535 1490 assert(t.link); … … 1553 1508 Vertex &A=vertices[ks]; 1554 1509 Real8 aa,bb,cc,dd; 1555 if ((dd=Area2(v0.r,v1.r,A.r)) >=0) 1556 {// warning PB roundoff error1510 if ((dd=Area2(v0.r,v1.r,A.r)) >=0){ 1511 // warning PB roundoff error 1557 1512 if (t.link && ( (aa=Area2( A.r , t[1].r , t[2].r )) < 0.0 1558 1513 || (bb=Area2( t[0].r , A.r , t[2].r )) < 0.0 1559 || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)) 1560 ferr++, cerr << " Error : " << ke + nbvold << " not in triangle " 1561 << i << " In=" << !!t.link 1562 << " " << aa << " " << bb << " " << cc << " " << dd << endl; 1563 1564 } 1565 1566 else 1567 { 1514 || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)){ 1515 printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd); 1516 throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); 1517 } 1518 } 1519 else { 1568 1520 if (tt.link && ( (aa=Area2( A.r , tt[1].r , tt[2].r )) < 0 1569 1521 || (bb=Area2( tt[0].r , A.r , tt[2].r )) < 0 1570 || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)) 1571 ferr++, cerr << " Warning : " << ke + nbvold << " not in triangle " << ii 1572 << " In=" << !!tt.link 1573 << " " << aa << " " << bb << " " << cc << " " << dd << endl; 1574 1575 } 1576 1522 || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)){ 1523 printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd); 1524 throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); 1525 } 1526 } 1577 1527 } 1578 1528 } 1579 } 1580 if(ferr) 1581 { 1582 cerr << " Number of triangles with P2 interpolation Probleme " << ferr << endl;; 1583 MeshError(9); 1584 } 1585 1586 for (i=0;i<nbt;i++) 1587 { 1529 } 1530 1531 for (i=0;i<nbt;i++){ 1588 1532 ksplit[i]=1; // no split by default 1589 1533 const Triangle & t = triangles[ i]; … … 1966 1910 // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) 1967 1911 OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); 1968 // if(CurrentTh)1969 // cout << "swap " << CurrentTh->Number(sa) << " " << CurrentTh->Number(sb) << " " ;1970 // cout << cosb12 << " " << sinba2 << " " << cosba2 << " " << sinb121971 // << " Onswap = " << OnSwap << endl;1972 1912 break; 1973 1913 } … … 2073 2013 Icoor2 detOld = t->det; 2074 2014 2075 if ( ( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ) 2076 { 2077 cerr << " infv " << infv << " det = " << detOld << endl; 2078 cerr << Number(s) << " "<< Number(s0) << " " 2079 << Number(s1) << " " << Number(s2) << endl; 2080 MeshError(3); 2081 } 2015 if (( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ){ 2016 throw ErrorException(__FUNCT__,exprintf("infv=%g det=%g")); 2017 } 2082 2018 2083 2019 // if det3 do not exist then constuct det3 … … 2114 2050 }} 2115 2051 else { 2116 cerr << " bug " << nbd0 <<endl; 2117 cerr << " Bug double points in " << endl ; 2118 cerr << " s = " << Number(s) << " " << s << endl; 2119 cerr << " s0 = "<< Number(s0) << " " << s0 << endl; 2120 cerr << " s1 = "<< Number(s1) << " " << s1 << endl; 2121 cerr << " s2 = "<< Number(s2) << " " << s2 << endl; 2122 MeshError(5,this);} 2123 2124 // remove de MarkUnSwap edge 2125 t->SetUnMarkUnSwap(0); 2126 t->SetUnMarkUnSwap(1); 2127 t->SetUnMarkUnSwap(2); 2128 2129 tt[0]= t; 2130 tt[1]= &triangles[nbt++]; 2131 tt[2]= &triangles[nbt++]; 2132 2133 if (nbt>nbtx) { 2134 cerr << " No enougth triangles " << endl; 2135 MeshError(999,this); 2136 } 2137 2138 *tt[1]= *tt[2]= *t; 2139 // gestion of the link 2140 tt[0]->link=tt[1]; 2141 tt[1]->link=tt[2]; 2142 2143 (* tt[0])(OppositeVertex[0])=&s; 2144 (* tt[1])(OppositeVertex[1])=&s; 2145 (* tt[2])(OppositeVertex[2])=&s; 2146 2147 tt[0]->det=det3[0]; 2148 tt[1]->det=det3[1]; 2149 tt[2]->det=det3[2]; 2150 2151 // update adj des triangles externe 2152 tt[0]->SetAdjAdj(0); 2153 tt[1]->SetAdjAdj(1); 2154 tt[2]->SetAdjAdj(2); 2155 // update des adj des 3 triangle interne 2156 const int i0 = 0; 2157 const int i1= NextEdge[i0]; 2158 const int i2 = PreviousEdge[i0]; 2159 2160 tt[i0]->SetAdj2(i2,tt[i2],i0); 2161 tt[i1]->SetAdj2(i0,tt[i0],i1); 2162 tt[i2]->SetAdj2(i1,tt[i1],i2); 2163 2164 tt[0]->SetTriangleContainingTheVertex(); 2165 tt[1]->SetTriangleContainingTheVertex(); 2166 tt[2]->SetTriangleContainingTheVertex(); 2167 2168 2169 // swap if the point s is on a edge 2170 if(izerodet>=0) { 2171 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2172 int rswap =tt[izerodet]->swap(iedge); 2173 2174 if (!rswap) 2175 { 2176 cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2177 } 2178 assert(rswap); 2179 } 2052 printf("bug (%i): Bug double points in\n",nbd0); 2053 throw ErrorException(__FUNCT__,exprintf("See above")); 2054 } 2055 2056 // remove de MarkUnSwap edge 2057 t->SetUnMarkUnSwap(0); 2058 t->SetUnMarkUnSwap(1); 2059 t->SetUnMarkUnSwap(2); 2060 2061 tt[0]= t; 2062 tt[1]= &triangles[nbt++]; 2063 tt[2]= &triangles[nbt++]; 2064 2065 if (nbt>nbtx) { 2066 throw ErrorException(__FUNCT__,exprintf("Not ebough triangles")); 2067 } 2068 2069 *tt[1]= *tt[2]= *t; 2070 // gestion of the link 2071 tt[0]->link=tt[1]; 2072 tt[1]->link=tt[2]; 2073 2074 (* tt[0])(OppositeVertex[0])=&s; 2075 (* tt[1])(OppositeVertex[1])=&s; 2076 (* tt[2])(OppositeVertex[2])=&s; 2077 2078 tt[0]->det=det3[0]; 2079 tt[1]->det=det3[1]; 2080 tt[2]->det=det3[2]; 2081 2082 // update adj des triangles externe 2083 tt[0]->SetAdjAdj(0); 2084 tt[1]->SetAdjAdj(1); 2085 tt[2]->SetAdjAdj(2); 2086 // update des adj des 3 triangle interne 2087 const int i0 = 0; 2088 const int i1= NextEdge[i0]; 2089 const int i2 = PreviousEdge[i0]; 2090 2091 tt[i0]->SetAdj2(i2,tt[i2],i0); 2092 tt[i1]->SetAdj2(i0,tt[i0],i1); 2093 tt[i2]->SetAdj2(i1,tt[i1],i2); 2094 2095 tt[0]->SetTriangleContainingTheVertex(); 2096 tt[1]->SetTriangleContainingTheVertex(); 2097 tt[2]->SetTriangleContainingTheVertex(); 2098 2099 2100 // swap if the point s is on a edge 2101 if(izerodet>=0) { 2102 // cout << " the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2103 int rswap =tt[izerodet]->swap(iedge); 2104 2105 if (!rswap) 2106 { 2107 cout << " Pb swap the point s is on a edge =>swap " << iedge << " " << *tt[izerodet] << endl; 2108 } 2109 assert(rswap); 2110 } 2180 2111 } 2181 2112 /*}}}1*/ … … 2187 2118 Int4 nbvold=nbv; 2188 2119 long int verbosity=2; 2189 for (it=0;it<nbt;it++) 2190 { 2120 for (it=0;it<nbt;it++){ 2191 2121 Triangle &t=triangles[it]; 2192 2122 if (t.link) … … 2214 2144 } 2215 2145 } 2216 2146 } 2217 2147 ReMakeTriangleContainingTheVertex(); 2218 if (nbvold!=nbv) 2219 { 2148 if (nbvold!=nbv){ 2220 2149 Int4 iv = nbvold; 2221 2150 Int4 NbSwap = 0; 2222 2151 Icoor2 dete[3]; 2223 for (Int4 i=nbvold;i<nbv;i++) 2224 {// for all the new point 2152 for (Int4 i=nbvold;i<nbv;i++) {// for all the new point 2225 2153 Vertex & vi = vertices[i]; 2226 2154 vi.i = toI2(vi.r); … … 2234 2162 Triangle *tcvi = FindTriangleContening(vi.i,dete); 2235 2163 if (tcvi && !tcvi->link) { 2236 cout << i << " PB insert point " << Number(vi) << vi << Number(vi) 2237 << " tcvi = " << tcvi << " " << tcvi->link << endl; 2238 cout << (*tcvi)[1] << (*tcvi)[2] << endl; 2239 tcvi = FindTriangleContening(vi.i,dete); 2240 cout << (*tcvi)[1] << (*tcvi)[2] << endl; 2241 MeshError(1001,this); 2164 printf("problem inserting point\n"); 2242 2165 } 2243 2244 2166 2245 2167 quadtree->Add(vi); … … 2249 2171 iv++; 2250 2172 // } 2251 } 2252 if (verbosity>3) 2253 { 2254 cout << " Nb Of New Point " << iv ; 2255 cout << " Nb swap = " << NbSwap << " to split internal edges with border vertices" ;} 2256 2173 } 2174 if (verbosity>3) { 2175 printf(" number of points: %i\n",iv); 2176 printf(" number of swap to split internal edges with border vertices: %i\n",NbSwap); 2257 2177 nbv = iv; 2258 } 2259 if (NbSplitEdge > nbv-nbvold) 2260 cout << " Warning not enough vertices to split all internal edges " << endl 2261 << " we lost " << NbSplitEdge - ( nbv-nbvold) << " Edges Sorry " << endl; 2262 if (verbosity>2) 2263 cout << "SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge << endl; 2178 } 2179 } 2180 if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold)); 2181 if (verbosity>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge); 2182 2264 2183 return NbSplitEdge; 2265 2184 } … … 2315 2234 // << " " << Number(vi) << " <--> " << Number(vj) <<endl; 2316 2235 Triangle *tcvj = FindTriangleContening(vj.i,dete); 2317 if (tcvj && !tcvj->link) 2318 { 2319 cerr << i << " PB insert point " << Number(vj) << vj << Number(vi) 2320 << " tcvj = " << tcvj << " " << tcvj->link << endl; 2321 cerr << (*tcvj)[1] << (*tcvj)[2] << endl; 2322 tcvj = FindTriangleContening(vj.i,dete); 2323 cout << (*tcvj)[1] << (*tcvj)[2] << endl; 2324 MeshError(1001,this); 2325 } 2326 2327 2236 if (tcvj && !tcvj->link){ 2237 throw ErrorException(__FUNCT__,exprintf("problem inserting point")); 2238 } 2328 2239 quadtree->Add(vj); 2329 2240 assert (tcvj && tcvj->det >= 0) ;// internal … … 2334 2245 } 2335 2246 if (verbosity>3) { 2336 cout << " Nb Of New Point " << iv << " Nb Of To close Points " << nbv-iv ; 2337 cout << " Nb swap = " << NbSwap << " after " ;} 2338 2247 printf(" number of new points: %i\n",iv); 2248 printf(" number of to close (?) points: %i\n",nbv-iv); 2249 printf(" number of swap after: %i\n",NbSwap); 2250 } 2339 2251 nbv = iv; 2340 2252 } 2341 2253 2342 for (i=nbvold;i<nbv;i++) 2343 NbSwap += vertices[i].Optim(1); 2344 if (verbosity>3) 2345 cout << " NbSwap = " << NbSwap << endl; 2346 2254 for (i=nbvold;i<nbv;i++) NbSwap += vertices[i].Optim(1); 2255 if (verbosity>3) printf(" NbSwap=%i\n",NbSwap); 2347 2256 2348 2257 NbTSwap += NbSwap ; 2349 2258 return nbv-nbvold; 2350 2351 2259 } 2352 2260 /*}}}1*/ … … 2864 2772 SetIntCoor(); 2865 2773 Int4 i; 2866 for (i=0;i<nbv;i++) 2867 ordre[i]= &vertices[i] ; 2774 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ; 2868 2775 2869 2776 // construction d'un ordre aleatoire … … 2873 2780 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 2874 2781 2875 2876 2877 2878 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 2879 if ( ++i >= nbv) { 2880 cerr << " All the vertices are aline " << endl; 2881 MeshError(998,this); } 2882 2883 // echange i et 2 dans ordre afin 2884 // que les 3 premiers ne soit pas aligne 2885 Exchange( ordre[2], ordre[i]); 2886 2887 // on ajoute un point a l'infini pour construire le maillage 2888 // afin d'avoir une definition simple des aretes frontieres 2889 nbt = 2; 2890 2891 2892 // on construit un maillage trivale forme 2893 // d'une arete et de 2 triangles 2894 // construit avec le 2 aretes orientes et 2895 Vertex * v0=ordre[0], *v1=ordre[1]; 2896 2897 triangles[0](0) = 0; // sommet pour infini 2898 triangles[0](1) = v0; 2899 triangles[0](2) = v1; 2900 2901 triangles[1](0) = 0;// sommet pour infini 2902 triangles[1](2) = v0; 2903 triangles[1](1) = v1; 2904 const int e0 = OppositeEdge[0]; 2905 const int e1 = NextEdge[e0]; 2906 const int e2 = PreviousEdge[e0]; 2907 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 2908 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 2909 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 2910 2911 triangles[0].det = -1; // faux triangles 2912 triangles[1].det = -1; // faux triangles 2913 2914 triangles[0].SetTriangleContainingTheVertex(); 2915 triangles[1].SetTriangleContainingTheVertex(); 2916 2917 triangles[0].link=&triangles[1]; 2918 triangles[1].link=&triangles[0]; 2919 2920 // nbtf = 2; 2921 if ( !quadtree ) quadtree = new QuadTree(this,0); 2922 quadtree->Add(*v0); 2923 quadtree->Add(*v1); 2924 2925 // on ajoute les sommets un Ò un 2926 Int4 NbSwap=0; 2927 2928 time1=CPUtime(); 2929 2930 if (verbosity>3) cout << " -- Begin of insertion process " << endl; 2931 2932 for (Int4 icount=2; icount<nbv; icount++) { 2933 Vertex *vi = ordre[icount]; 2934 // cout << " Insert " << Number(vi) << endl; 2935 Icoor2 dete[3]; 2936 Triangle *tcvi = FindTriangleContening(vi->i,dete); 2937 quadtree->Add(*vi); 2938 Add(*vi,tcvi,dete); 2939 NbSwap += vi->Optim(1,0); 2940 2941 }// fin de boucle en icount 2942 time2=CPUtime(); 2943 if (verbosity>3) 2944 cout << " NbSwap of insertion " << NbSwap 2945 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2946 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " 2947 << (float)NbUnSwap /(float) nbv 2948 <<endl; 2949 NbUnSwap = 0; 2950 // construction d'un ordre aleatoire 2951 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; 2782 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){ 2783 if ( ++i >= nbv) { 2784 throw ErrorException(__FUNCT__,exprintf("all the vertices are aligned")); 2785 } 2786 } 2787 // echange i et 2 dans ordre afin 2788 // que les 3 premiers ne soit pas aligne 2789 Exchange( ordre[2], ordre[i]); 2790 2791 // on ajoute un point a l'infini pour construire le maillage 2792 // afin d'avoir une definition simple des aretes frontieres 2793 nbt = 2; 2794 2795 // on construit un maillage trivale forme 2796 // d'une arete et de 2 triangles 2797 // construit avec le 2 aretes orientes et 2798 Vertex * v0=ordre[0], *v1=ordre[1]; 2799 2800 triangles[0](0) = 0; // sommet pour infini 2801 triangles[0](1) = v0; 2802 triangles[0](2) = v1; 2803 2804 triangles[1](0) = 0;// sommet pour infini 2805 triangles[1](2) = v0; 2806 triangles[1](1) = v1; 2807 const int e0 = OppositeEdge[0]; 2808 const int e1 = NextEdge[e0]; 2809 const int e2 = PreviousEdge[e0]; 2810 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 2811 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 2812 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 2813 2814 triangles[0].det = -1; // faux triangles 2815 triangles[1].det = -1; // faux triangles 2816 2817 triangles[0].SetTriangleContainingTheVertex(); 2818 triangles[1].SetTriangleContainingTheVertex(); 2819 2820 triangles[0].link=&triangles[1]; 2821 triangles[1].link=&triangles[0]; 2822 2823 // nbtf = 2; 2824 if ( !quadtree ) quadtree = new QuadTree(this,0); 2825 quadtree->Add(*v0); 2826 quadtree->Add(*v1); 2827 2828 // on ajoute les sommets un Ò un 2829 Int4 NbSwap=0; 2830 2831 time1=CPUtime(); 2832 2833 if (verbosity>3) cout << " -- Begin of insertion process " << endl; 2834 2835 for (Int4 icount=2; icount<nbv; icount++) { 2836 Vertex *vi = ordre[icount]; 2837 // cout << " Insert " << Number(vi) << endl; 2838 Icoor2 dete[3]; 2839 Triangle *tcvi = FindTriangleContening(vi->i,dete); 2840 quadtree->Add(*vi); 2841 Add(*vi,tcvi,dete); 2842 NbSwap += vi->Optim(1,0); 2843 2844 }// fin de boucle en icount 2845 time2=CPUtime(); 2846 if (verbosity>3) 2847 cout << " NbSwap of insertion " << NbSwap 2848 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2849 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " 2850 << (float)NbUnSwap /(float) nbv 2851 <<endl; 2852 NbUnSwap = 0; 2853 // construction d'un ordre aleatoire 2854 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; 2952 2855 #ifdef NBLOOPOPTIM 2953 2856 2954 2955 2956 2957 2958 2959 2960 2961 2962 2963 2964 2965 2966 2967 2968 2969 2970 2971 2972 2973 2974 2975 2976 2857 k3 = rand()%nbv ; 2858 for (int is4=0; is4<nbv; is4++) 2859 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 2860 2861 double timeloop = time2 ; 2862 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 2863 { 2864 double time000 = timeloop; 2865 Int4 NbSwap = 0; 2866 for (int is1=0; is1<nbv; is1++) 2867 NbSwap += ordre[is1]->Optim(0,0); 2868 timeloop = CPUtime(); 2869 if (verbosity>3) 2870 cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap 2871 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2872 << " CPU=" << timeloop - time000 << " s, " 2873 << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv 2874 << endl; 2875 NbUnSwap = 0; 2876 if(!NbSwap) break; 2877 } 2878 ReMakeTriangleContainingTheVertex(); 2879 // because we break the TriangleContainingTheVertex 2977 2880 #endif 2978 2979 2980 2981 2982 2983 2984 2985 2881 time3=CPUtime(); 2882 if (verbosity>4) 2883 cout << " init " << time1 - time0 << " initialisation, " 2884 << time2 - time1 << "s, insert point " 2885 << time3 -time2 << "s, optim " << endl 2886 << " Init Total Cpu Time = " << time3 - time0 << "s " << endl; 2887 2888 CurrentTh=OldCurrentTh; 2986 2889 } 2987 2890 /*}}}1*/ … … 2997 2900 k++,cerr << " det T" << t << " = " << 0 << endl; 2998 2901 if (k!=0) { 2999 cerr << " ther is " << k << " triangles of mes = 0 " << endl;3000 MeshError(11,this);}2902 throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k)); 2903 } 3001 2904 3002 2905 TriangleAdjacent ta(0,0); … … 3026 2929 3027 2930 if (k!=0) { 3028 cerr << " they is " << k << " lost edges " << endl; 3029 cerr << " The boundary is crossing may be!" << endl; 3030 MeshError(10,this); 2931 throw ErrorException(__FUNCT__,exprintf("There are %i lost edges, the boundary might be crossing",k)); 3031 2932 } 3032 2933 for (Int4 j=0;j<nbv;j++) … … 3118 3019 } 3119 3020 it++;} // end while (it<nbt) 3120 if (nbt == NbOutT || !NbSubDomTot) 3121 { 3122 cout << "\n error : " << NbOutT << " " << NbSubDomTot <<" " << nbt << endl; 3123 cerr << "Error: The boundary is not close => All triangles are outside " << endl; 3124 MeshError(888,this); 3125 } 3021 if (nbt == NbOutT || !NbSubDomTot) { 3022 throw ErrorException(__FUNCT__,exprintf("The boundary is not close: all triangles are outside")); 3023 } 3126 3024 3127 3025 delete [] HeapArete; … … 3233 3131 mark[it]=triangles[it].link ? -1 : -2; 3234 3132 Int4 inew =0; 3235 for (Int4 i=0;i<NbSubDomains;i++) 3236 { 3133 for (Int4 i=0;i<NbSubDomains;i++) { 3237 3134 GeometricalEdge &eg = *Gh.subdomains[i].edge; 3238 3135 subdomains[i].ref = Gh.subdomains[i].ref; … … 3258 3155 TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges 3259 3156 3260 while (1) 3261 { 3157 while (1) { 3262 3158 assert( v0 == ta.EdgeVertex(1) ); 3263 3159 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl; … … 3268 3164 subdomains[i].head=t=ta; 3269 3165 //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl; 3270 if(t<triangles || t >= triangles+nbt || t->det < 0 3271 || t->link == 0) // Ajoute aout 200 3166 if(t<triangles || t >= triangles+nbt || t->det < 0 || t->link == 0) { 3167 throw ErrorException(__FUNCT__,exprintf("bad definition of SubSomain %i",i)); 3168 } 3169 Int4 it = Number(t); 3170 if (mark[it] >=0) { 3171 if(verbosity>10){ 3172 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref 3173 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref 3174 << " skip this def " << endl; 3175 } 3176 break; 3177 } 3178 if(i != inew) 3179 Exchange(subdomains[i],subdomains[inew]); 3180 inew++; 3181 Triangle *tt=t; 3182 Int4 kkk=0; 3183 do 3272 3184 { 3273 cerr << " Error in the def of sub domain "<<i 3274 << " form border " << NbSubDomains - i << "/" << NbSubDomains 3275 << ": Bad sens " << Gh.Number(eg) <<" "<< sens << endl; 3276 err++; 3277 break;} 3278 Int4 it = Number(t); 3279 if (mark[it] >=0) { 3280 if(verbosity>10) 3281 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref 3282 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref 3283 << " skip this def " << endl; 3284 break;} 3285 if(i != inew) 3286 Exchange(subdomains[i],subdomains[inew]); 3287 inew++; 3288 Triangle *tt=t; 3289 Int4 kkk=0; 3290 do 3291 { 3292 kkk++; 3293 assert(mark[Number(tt)]<0); 3294 mark[Number(tt)]=i; 3295 tt=tt->link; 3296 } while (tt!=t); 3297 if(verbosity>7) 3298 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl; 3299 break;} 3300 ta = Previous(Adj(ta)); 3301 if(t == (Triangle *) ta) { 3302 err++; 3303 cerr << " Error in the def of sub domain " << i 3304 << " edge=" << Gh.Number(eg) << " " << sens << endl; 3305 break;} 3306 // cout << " NB of remove subdomain " << NbSubDomTot-NbSubDomains<< endl; 3307 3308 } 3309 3310 } 3311 if (err) MeshError(777,this); 3185 kkk++; 3186 assert(mark[Number(tt)]<0); 3187 mark[Number(tt)]=i; 3188 tt=tt->link; 3189 } while (tt!=t); 3190 if(verbosity>7) 3191 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl; 3192 break; 3193 } 3194 ta = Previous(Adj(ta)); 3195 if(t == (Triangle *) ta) { 3196 throw ErrorException(__FUNCT__,exprintf("bad definition of SubSomain %i",i)); 3197 } 3198 } 3199 } 3312 3200 3313 3201 if (inew < NbSubDomains) { … … 3549 3437 VertexOnBThVertex = new VertexOnVertex[NbVerticesOnGeomVertex]; 3550 3438 // 3551 if( NbVerticesOnGeomVertex >= nbvx) 3552 { 3553 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; 3554 MeshError(1,this); 3555 } 3439 if( NbVerticesOnGeomVertex >= nbvx) { 3440 throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); 3441 } 3556 3442 assert(vertices); 3557 3443 for (i=0;i<Gh.nbv;i++) … … 3764 3650 assert (A1-vertices>=0 && A1-vertices <nbv); 3765 3651 break;} 3766 if (!ee.adj[k1]) 3767 {cerr << "Error adj edge " << BTh.Number(ee) << ", nbe = " << nbe 3768 << " Gh.vertices " << Gh.vertices 3769 << " k1 = " << k1 << " on=" << *ee[k1].on << endl; 3770 cerr << ee[k1].on->gv-Gh.vertices << endl; 3771 } 3652 if (!ee.adj[k1]) { 3653 throw ErrorException(__FUNCT__,exprintf(" adj edge %i, nbe=%i, Gh.vertices=%i",BTh.Number(ee),nbe,Gh.vertices)); 3654 } 3772 3655 pe = ee.adj[k1]; // next edge 3773 3656 k0 = pe->Intersection(ee); … … 3823 3706 // end new code 3824 3707 // do the allocation 3825 if(step==0) 3826 { 3708 if(step==0){ 3827 3709 //if(!NbOfNewPoints) break;// nothing ????? bug 3828 if(nbv+NbOfNewPoints > nbvx) 3829 { 3830 cerr << " Too much vertices on geometry " << nbv+NbOfNewPoints << " >= " << nbvx << endl; 3831 MeshError(3,this); 3832 } 3833 //cout << " NbOfNewEdge" << NbOfNewEdge << " NbOfNewPoints " << NbOfNewPoints << endl; 3710 if(nbv+NbOfNewPoints > nbvx) { 3711 throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",nbv+NbOfNewPoints,nbvx)); 3712 } 3834 3713 edges = new Edge[NbOfNewEdge]; 3835 3714 nbex = NbOfNewEdge; … … 3846 3725 3847 3726 delete [] bcurve; 3848 3849 3727 3850 3728 Insert(); … … 3878 3756 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 3879 3757 // 3880 if( NbVerticesOnGeomVertex >= nbvx) 3881 { 3882 cerr << " Too much vertices on geometry " << NbVerticesOnGeomVertex << " >= " << nbvx << endl; 3883 MeshError(1,this); 3884 } 3758 if( NbVerticesOnGeomVertex >= nbvx) { 3759 throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); 3760 } 3885 3761 for (i=0;i<Gh.nbv;i++) 3886 3762 if (Gh[i].Required()&& Gh[i].IsThe() ) {//Gh vertices Required … … 4157 4033 // assert( e[i]); 4158 4034 } 4159 if(kk) MeshError(997,this);4035 if(kk) throw ErrorException(__FUNCT__,exprintf("See above")); 4160 4036 4161 4037 return e; … … 4188 4064 // computation of the det 4189 4065 int Nberr=0; 4190 for (i=0;i<nbt;i++) 4191 { 4066 for (i=0;i<nbt;i++) { 4192 4067 Vertex & v0 = triangles[i][0]; 4193 4068 Vertex & v1 = triangles[i][1]; … … 4196 4071 { 4197 4072 triangles[i].det= det(v0,v1,v2); 4198 if (triangles[i].det <=0 && Nberr++ <10) 4199 { 4073 if (triangles[i].det <=0 && Nberr++ <10){ 4200 4074 if(Nberr==1) 4201 if (strfrom) 4202 cerr << "+++ Fatal Error " << strfrom << "(SetInCoor) Error : area of Triangle < 0 " << endl; 4203 else 4204 cerr << "+++ Fatal Error Triangle (in SetInCoor) area of Triangle < 0" << endl; 4205 cerr << " Triangle " << i << " det (I2) = " << triangles[i].det ; 4206 cerr << " (R2) " << Det(v1.r-v0.r,v2.r-v0.r); 4207 cerr << "; The 3 vertices " << endl; 4208 cerr << Number(v0) << " " << Number(v1) << " " 4209 << Number(v2) << " : " ; 4210 cerr << v0.r << v1.r << v2.r << " ; "; 4211 cerr << v0.i << v1.i << v2.i << endl; 4212 } 4075 if (strfrom){ 4076 throw ErrorException(__FUNCT__,exprintf("Fatal error %s (SetInCoor) : area of Triangle %i < 0",strfrom,i)); 4077 } 4078 else{ 4079 throw ErrorException(__FUNCT__,exprintf("Fatal error (SetInCoor) : area of Triangle %i < 0",i)); 4080 } 4081 } 4213 4082 } 4214 else 4215 triangles[i].det= -1; // Null triangle; 4216 } 4217 if (Nberr) MeshError(899,this); 4218 4083 else triangles[i].det= -1; // Null triangle; 4084 } 4219 4085 } 4220 4086 /*}}}1*/ … … 4338 4204 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;) 4339 4205 if ( ++i >= nbvb) { 4340 cerr << "FillHoleInMesh: All the vertices are aline " << nbvb << endl;4341 MeshError(998,this);}4342 4206 throw ErrorException(__FUNCT__,exprintf("FillHoleInMesh: All the vertices are aligned")); 4207 } 4208 Exchange( ordre[2], ordre[i]); 4343 4209 4344 4210 Vertex * v0=ordre[0], *v1=ordre[1]; … … 4402 4268 } 4403 4269 if(nbloss) { 4404 cerr << " we loss some " << nbloss << " " << " edges other " << knbe << endl; 4405 MeshError(1100,this); 4270 throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe)); 4406 4271 } 4407 4272 … … 4490 4355 */ 4491 4356 if (k) { 4492 cerr << "Error Nb of triangles edge alone = " << k << endl; 4493 MeshError(9997,this); 4357 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); 4494 4358 } 4495 4359 FindSubDomain(); … … 4882 4746 SetVertexFieldOn(); 4883 4747 4884 4885 if (vlast >= vend) 4886 { 4887 cerr << " Not enougth vertices to crack the mesh we need " << nbv << " vertices " << endl; 4888 MeshError(555,this); 4748 if (vlast >= vend) { 4749 throw ErrorException(__FUNCT__,exprintf("Not enougth vertices: to crack the mesh we need %i vertices",nbv)); 4889 4750 } 4890 4751 cout << " NbCrackedVertices " << NbCrackedVertices << endl; … … 4905 4766 4906 4767 if (! a || !a->t ) { 4907 if (a) 4908 {cerr << " Attention PB TriangleConteningTheVertex vertex number=" << Number(a) << endl; 4909 cerr << "We forget a call to ReMakeTriangleContainingTheVertex" << endl;} 4910 cerr << " Pb with " << B << toR2(B) << endl; 4911 MeshError(7777); 4768 if (a) { 4769 printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a)); 4770 } 4771 throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening")); 4912 4772 } 4913 4773 assert(a>= vertices && a < vertices+nbv); … … 6025 5885 } 6026 5886 /*}}}1*/ 6027 5887 /*FUNCTION AGoodNumberPrimeWith{{{1*/ 6028 5888 Int4 AGoodNumberPrimeWith(Int4 n){ 6029 5889 const Int4 BigPrimeNumber[] ={ 567890359L, … … 6042 5902 return pi; 6043 5903 } 6044 5904 /*}}}1*/ 5905 /*FUNCTION MeshError{{{1*/ 6045 5906 void MeshError(int Err,Triangles *Th){ 6046 5907 cerr << " Fatal error in the meshgenerator " << Err << endl ; 6047 5908 exit(1); 6048 5909 } 6049 5910 /*}}}1*/ 5911 /*FUNCTION ostream& operator{{{1*/ 6050 5912 ostream& operator <<(ostream& f, const Triangle & ta) { 6051 5913 if(CurrentTh) … … 6067 5929 << "{" << ta.at[2] << " " << ta.aa[2] << "} " 6068 5930 << "]" ; 6069 return f;} 6070 6071 6072 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , 6073 TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) 6074 { // l'arete ta coupe l'arete pva pvb 6075 // de cas apres le swap sa coupe toujours 6076 // on cherche l'arete suivante 6077 // on suppose que detsa >0 et detsb <0 6078 // attention la routine echange pva et pvb 6079 6080 if(tt1.Locked()) return 0; // frontiere croise 6081 6082 TriangleAdjacent tt2 = Adj(tt1); 6083 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 6084 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 6085 assert ( a1 >= 0 && a1 < 3 ); 6086 6087 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 6088 Vertex & s1= (*t1)[OppositeVertex[a1]]; 6089 Vertex & s2= (*t2)[OppositeVertex[a2]]; 6090 6091 6092 Icoor2 dets2 = det(*pva,*pvb,s2); 6093 Icoor2 det1=t1->det , det2=t2->det ; 6094 Icoor2 detT = det1+det2; 6095 assert((det1>0 ) && (det2 > 0)); 6096 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb 6097 Icoor2 ndet1 = bamg::det(s1,sa,s2); 6098 Icoor2 ndet2 = detT - ndet1; 6099 6100 int ToSwap =0; //pas de swap 6101 if ((ndet1 >0) && (ndet2 >0)) 6102 { // on peut swaper 6103 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 6104 ToSwap =1; 6105 else // swap alleatoire 6106 if (BinaryRand()) 6107 ToSwap =2; 6108 } 6109 if (ToSwap) NbSwap++, 6110 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 6111 6112 int ret=1; 6113 6114 if (dets2 < 0) {// haut 6115 dets1 = ToSwap ? dets1 : detsa ; 6116 detsa = dets2; 6117 tt1 = Previous(tt2) ;} 6118 else if (dets2 > 0){// bas 6119 dets1 = ToSwap ? dets1 : detsb ; 6120 detsb = dets2; 6121 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 6122 if(!ToSwap) tt1 = Next(tt2); 5931 return f; 5932 } 5933 /*}}}1*/ 5934 /*FUNCTION SwapForForcingEdge{{{1*/ 5935 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) { 5936 // l'arete ta coupe l'arete pva pvb 5937 // de cas apres le swap sa coupe toujours 5938 // on cherche l'arete suivante 5939 // on suppose que detsa >0 et detsb <0 5940 // attention la routine echange pva et pvb 5941 5942 if(tt1.Locked()) return 0; // frontiere croise 5943 5944 TriangleAdjacent tt2 = Adj(tt1); 5945 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 5946 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 5947 assert ( a1 >= 0 && a1 < 3 ); 5948 5949 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 5950 Vertex & s1= (*t1)[OppositeVertex[a1]]; 5951 Vertex & s2= (*t2)[OppositeVertex[a2]]; 5952 5953 5954 Icoor2 dets2 = det(*pva,*pvb,s2); 5955 Icoor2 det1=t1->det , det2=t2->det ; 5956 Icoor2 detT = det1+det2; 5957 assert((det1>0 ) && (det2 > 0)); 5958 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb 5959 Icoor2 ndet1 = bamg::det(s1,sa,s2); 5960 Icoor2 ndet2 = detT - ndet1; 5961 5962 int ToSwap =0; //pas de swap 5963 if ((ndet1 >0) && (ndet2 >0)) 5964 { // on peut swaper 5965 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 5966 ToSwap =1; 5967 else // swap alleatoire 5968 if (BinaryRand()) 5969 ToSwap =2; 5970 } 5971 if (ToSwap) NbSwap++, 5972 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 5973 5974 int ret=1; 5975 5976 if (dets2 < 0) {// haut 5977 dets1 = ToSwap ? dets1 : detsa ; 5978 detsa = dets2; 5979 tt1 = Previous(tt2) ;} 5980 else if (dets2 > 0){// bas 5981 dets1 = ToSwap ? dets1 : detsb ; 5982 detsb = dets2; 5983 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 5984 if(!ToSwap) tt1 = Next(tt2); 5985 } 5986 else { // changement de sens 5987 ret = -1; 5988 Exchange(pva,pvb); 5989 Exchange(detsa,detsb); 5990 Exchange(dets1,dets2); 5991 Exchange(tt1,tt2); 5992 dets1=-dets1; 5993 dets2=-dets2; 5994 detsa=-detsa; 5995 detsb=-detsb; 5996 5997 if (ToSwap) 5998 if (dets2 < 0) {// haut 5999 dets1 = (ToSwap ? dets1 : detsa) ; 6000 detsa = dets2; 6001 tt1 = Previous(tt2) ;} 6002 else if (dets2 > 0){// bas 6003 dets1 = (ToSwap ? dets1 : detsb) ; 6004 detsb = dets2; 6005 if(!ToSwap) tt1 = Next(tt2); 6006 } 6007 else {// on a fin ??? 6008 tt1 = Next(tt2); 6009 ret =0;} 6010 6011 } 6012 return ret; 6013 } 6014 /*}}}1*/ 6015 /*FUNCTION ForceEdge{{{1*/ 6016 6017 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) { 6018 int NbSwap =0; 6019 assert(a.t && b.t); // the 2 vertex is in a mesh 6020 int k=0; 6021 taret=TriangleAdjacent(0,0); // erreur 6022 6023 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); 6024 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 6025 // we turn around a in the direct sens 6026 6027 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 6028 if(v2) // normal case 6029 det2 = det(*v2,a,b); 6030 else { // no chance infini vertex try the next 6031 tta= Previous(Adj(tta)); 6032 v2 = tta.EdgeVertex(0); 6033 vbegin =v2; 6034 assert(v2); 6035 det2 = det(*v2,a,b); 6036 // cout << " No Change try the next" << endl; 6037 } 6038 6039 while (v2 != &b) { 6040 TriangleAdjacent tc = Previous(Adj(tta)); 6041 v1 = v2; 6042 v2 = tc.EdgeVertex(0); 6043 det1 = det2; 6044 det2 = v2 ? det(*v2,a,b): det2; 6045 6046 if((det1 < 0) && (det2 >0)) { 6047 // try to force the edge 6048 Vertex * va = &a, *vb = &b; 6049 tc = Previous(tc); 6050 assert ( v1 && v2); 6051 Icoor2 detss = 0,l=0,ks; 6052 // cout << "Real ForcingEdge " << *va << *vb << detss << endl; 6053 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 6054 if(l++ > 10000000) { 6055 throw ErrorException(__FUNCT__,exprintf("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i)too big",NbSwap,l)); 6056 } 6057 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 6058 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { 6059 tc.SetLock(); 6060 a.Optim(1,0); 6061 b.Optim(1,0); 6062 taret = tc; 6063 return NbSwap; 6064 } 6065 else 6066 { 6067 taret = tc; 6068 return -2; // error boundary is crossing 6069 } 6123 6070 } 6124 else { // changement de sens 6125 ret = -1; 6126 Exchange(pva,pvb); 6127 Exchange(detsa,detsb); 6128 Exchange(dets1,dets2); 6129 Exchange(tt1,tt2); 6130 dets1=-dets1; 6131 dets2=-dets2; 6132 detsa=-detsa; 6133 detsb=-detsb; 6134 6135 if (ToSwap) 6136 if (dets2 < 0) {// haut 6137 dets1 = (ToSwap ? dets1 : detsa) ; 6138 detsa = dets2; 6139 tt1 = Previous(tt2) ;} 6140 else if (dets2 > 0){// bas 6141 dets1 = (ToSwap ? dets1 : detsb) ; 6142 detsb = dets2; 6143 if(!ToSwap) tt1 = Next(tt2); 6144 } 6145 else {// on a fin ??? 6146 tt1 = Next(tt2); 6147 ret =0;} 6148 6149 } 6150 return ret; 6151 } 6152 6153 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) 6154 { 6155 int NbSwap =0; 6156 assert(a.t && b.t); // the 2 vertex is in a mesh 6157 int k=0; 6158 taret=TriangleAdjacent(0,0); // erreur 6159 6160 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); 6161 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 6162 // we turn around a in the direct sens 6163 6164 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 6165 if(v2) // normal case 6166 det2 = det(*v2,a,b); 6167 else { // no chance infini vertex try the next 6168 tta= Previous(Adj(tta)); 6169 v2 = tta.EdgeVertex(0); 6170 vbegin =v2; 6171 assert(v2); 6172 det2 = det(*v2,a,b); 6173 // cout << " No Change try the next" << endl; 6174 } 6175 6176 while (v2 != &b) { 6177 TriangleAdjacent tc = Previous(Adj(tta)); 6178 v1 = v2; 6179 v2 = tc.EdgeVertex(0); 6180 det1 = det2; 6181 det2 = v2 ? det(*v2,a,b): det2; 6182 6183 if((det1 < 0) && (det2 >0)) { 6184 // try to force the edge 6185 Vertex * va = &a, *vb = &b; 6186 tc = Previous(tc); 6187 assert ( v1 && v2); 6188 Icoor2 detss = 0,l=0,ks; 6189 // cout << "Real ForcingEdge " << *va << *vb << detss << endl; 6190 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 6191 if(l++ > 10000000) { 6192 cerr << " Loop in forcing Egde AB" 6193 <<"\n vertex A " << a 6194 <<"\n vertex B " << b 6195 <<"\n nb de swap " << NbSwap 6196 <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl; 6197 6198 if ( CurrentTh ) 6199 cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; 6200 MeshError(990); 6201 } 6202 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 6203 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { 6204 tc.SetLock(); 6205 a.Optim(1,0); 6206 b.Optim(1,0); 6207 taret = tc; 6208 return NbSwap; 6209 } 6210 else 6211 { 6212 taret = tc; 6213 return -2; // error boundary is crossing 6214 } 6215 } 6216 tta = tc; 6217 assert(k++<2000); 6218 if ( vbegin == v2 ) return -1;// error 6219 } 6220 6221 tta.SetLock(); 6222 taret=tta; 6223 a.Optim(1,0); 6224 b.Optim(1,0); 6225 return NbSwap; 6226 } 6071 tta = tc; 6072 assert(k++<2000); 6073 if ( vbegin == v2 ) return -1;// error 6074 } 6075 6076 tta.SetLock(); 6077 taret=tta; 6078 a.Optim(1,0); 6079 b.Optim(1,0); 6080 return NbSwap; 6081 } 6082 /*}}}1*/ 6227 6083 6228 6084 } -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r2847 r2850 11 11 #include "../../include/macros.h" 12 12 #include "../../toolkits/toolkits.h" 13 14 #undef __FUNCT__ 15 #define __FUNCT__ "Vertex" 13 16 14 17 namespace bamg {
Note:
See TracChangeset
for help on using the changeset viewer.
