![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2965
- Timestamp:
- 02/05/10 11:05:19 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 5 edited
-
Mesh2.h (modified) (3 diffs)
-
objects/Geometry.cpp (modified) (6 diffs)
-
objects/ListofIntersectionTriangles.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (32 diffs)
-
objects/Vertex.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.h
r2958 r2965 114 114 Direction DirOfSearch; 115 115 union { 116 Triangle 116 Triangle* t; // one triangle which contained the vertex 117 117 Int4 color; 118 118 Vertex * to;// use in geometry Vertex to now the Mesh Vertex associed … … 366 366 } 367 367 368 void SetTriangleContainingTheVertex() 369 { 368 void SetTriangleContainingTheVertex() { 370 369 if (ns[0]) (ns[0]->t=this,ns[0]->vint=0); 371 370 if (ns[1]) (ns[1]->t=this,ns[1]->vint=1); 372 371 if (ns[2]) (ns[2]->t=this,ns[2]->vint=2); 373 372 } 374 373 375 374 int swap(Int2 a1,int=0); … … 782 781 783 782 Vertex * NearestVertex(Icoor1 i,Icoor1 j) ; 784 Triangle * FindTriangleCont ening(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;783 Triangle * FindTriangleContaining(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; 785 784 786 785 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2957 r2965 234 234 Real8 Hmin = HUGE_VAL;// the infinie value 235 235 Int4 hvertices=0; 236 int i,j,k,n ;236 int i,j,k,n,i1,i2; 237 237 238 238 //initialize some variables … … 290 290 //Edges 291 291 if (bamggeom->Edges){ 292 int i1,i2; 293 R2 zero2(0,0); 294 Real4 *len =0; 292 R2 zero2(0,0); 293 Real4* len=NULL; 295 294 296 295 if(verbose>3) printf(" processing Edges\n"); 297 296 edges = new GeometricalEdge[nbe]; 298 297 298 //if hvertices==0, initialize len (length of each edge) 299 299 if (!hvertices) { 300 300 len = new Real4[nbv]; 301 for(i=0;i<nbv;i++) 302 len[i]=0; 301 for(i=0;i<nbv;i++) len[i]=0; 303 302 } 304 303 305 304 for (i=0;i<nbe;i++){ 305 306 306 i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing 307 307 i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing 308 308 edges[i].ref=(Int4)bamggeom->Edges[i*3+2]; 309 edges[i].v[0]= vertices + i1; 310 edges[i].v[1]= vertices + i2; 311 R2 x12 = vertices[i2].r-vertices[i1].r; 309 310 edges[i].v[0]= vertices + i1; 311 edges[i].v[1]= vertices + i2; 312 313 //get length of edge 314 R2 x12=vertices[i2].r-vertices[i1].r; 312 315 Real8 l12=sqrt((x12,x12)); 316 Hmin=Min(Hmin,l12); 317 318 //initialize other fields 313 319 edges[i].tg[0]=zero2; 314 320 edges[i].tg[1]=zero2; … … 322 328 len[i2] += l12; 323 329 } 324 325 Hmin = Min(Hmin,l12);326 330 } 327 331 … … 330 334 for (i=0;i<nbv;i++) 331 335 if (vertices[i].color > 0) 332 vertices[i].m= 336 vertices[i].m=Metric(len[i] /(Real4) vertices[i].color); 333 337 else 334 vertices[i].m= 338 vertices[i].m=Metric(Hmin); 335 339 delete [] len; 336 340 } … … 460 464 j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing 461 465 if (j>nbe-1 || j<0){ 462 throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbv)); 463 } 464 edges[j].SetRequired(); } 466 throw ErrorException(__FUNCT__,exprintf("Bad RequiredEdges definition: should in [0 %i]",nbe)); 467 } 468 edges[j].SetRequired(); 469 } 465 470 } 466 471 else{ … … 741 746 } 742 747 743 // generation of all the tangent e748 // generation of all the tangents 744 749 for (i=0;i<nbe;i++) { 745 R2 AB = edges[i].v[1]->r -edges[i].v[0]->r; 746 Real8 lAB = Norme2(AB); // length of current edge AB 750 751 //Get AB = vertex1->vertex2 752 R2 AB =edges[i].v[1]->r -edges[i].v[0]->r; 753 //Get length of AB 754 Real8 lAB=Norme2(AB); 755 //initialize tangent 747 756 Real8 ltg2[2]; 748 757 ltg2[0]=0;ltg2[1]=0; 758 759 //loop over the 2 vertices of the edge 749 760 for (jj=0;jj<2;jj++) { 750 R2 tg = edges[i].tg[jj]; 751 Real8 ltg = Norme2(tg); // length of tg 752 if(ltg == 0) {// no tg 753 if( ! edges[i].v[jj]->Corner()) { // not a Corner 754 tg = edges[i].v[1-jj]->r 755 - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; 756 ltg = Norme2(tg); 761 R2 tg =edges[i].tg[jj]; 762 Real8 ltg=Norme2(tg); 763 764 //by default, tangent=[0 0] 765 if(ltg == 0) { 766 //if the current vertex of the edge is not a corner 767 if(!edges[i].v[jj]->Corner()){ 768 tg = edges[i].v[1-jj]->r - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; 769 ltg= Norme2(tg); 757 770 tg = tg *(lAB/ltg),ltg=lAB; 758 771 } -
issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
r2865 r2965 41 41 else {// not optimisation 42 42 init(); 43 t=tbegin = Bh.FindTriangleCont ening(a,deta);43 t=tbegin = Bh.FindTriangleContaining(a,deta); 44 44 if( t->det>=0) 45 45 ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det); -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2960 r2965 730 730 /*Methods*/ 731 731 /*FUNCTION Triangles::Add{{{1*/ 732 void Triangles::Add( Vertex & s,Triangle * t, Icoor2* det3) {732 void Triangles::Add( Vertex &s,Triangle* t, Icoor2* det3) { 733 733 // ------------------------------------------- 734 734 // s2 … … 746 746 //-------------------------------------------- 747 747 748 Triangle * tt[3]; // the 3 new Triangles 749 Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2]; 750 Icoor2 det3local[3]; 751 int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; 752 // infv = ordre of the infini vertex (null) 753 register int nbd0 =0; // number of zero det3 754 register int izerodet=-1,iedge; // izerodet = egde contening the vertex s 748 //the three triangles tt 749 Triangle* tt[3]; 750 //three vertices of t 751 Vertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2]; 752 //three determinants 753 Icoor2 det3local[3]; 754 // number of zero det3 755 register int nbzerodet =0; 756 // izerodet = egde contening the vertex s 757 register int izerodet=-1,iedge; 758 //determinant of t 755 759 Icoor2 detOld = t->det; 756 760 757 if (( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ){ 758 throw ErrorException(__FUNCT__,exprintf("infv=%g det=%g")); 759 } 760 761 // if det3 do not exist then constuct det3 761 // infinitevertexpos = order of the infinite vertex (NULL) 762 // if no infinite vertex (NULL) infinitevertexpos=-1 763 // else if v_i is infinite, infinitevertexpos=i 764 int infinitevertexpos = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; 765 766 //some checks 767 if (( infinitevertexpos <0 ) && (detOld <0) || ( infinitevertexpos >=0 ) && (detOld >0) ){ 768 throw ErrorException(__FUNCT__,exprintf("bug in Triangles::Add, bad configuration")); 769 } 770 771 // if det3 does not exist, build it 762 772 if (!det3) { 763 det3 = det3local; // alloc 764 if ( infv<0 ) { 773 //allocate 774 det3 = det3local; 775 //if no infinite vertex 776 if (infinitevertexpos<0 ) { 765 777 det3[0]=bamg::det(s ,s1,s2); 766 778 det3[1]=bamg::det(s0,s ,s2); … … 768 780 else { 769 781 // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk 770 det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; 771 det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; 772 det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}} 773 774 775 if (!det3[0]) izerodet=0,nbd0++; 776 if (!det3[1]) izerodet=1,nbd0++; 777 if (!det3[2]) izerodet=2,nbd0++; 778 779 if (nbd0 >0 ) // point s on a egde or on a vertex 780 if (nbd0 == 1) { 781 iedge = OppositeEdge[izerodet]; 782 TriangleAdjacent ta = t->Adj(iedge); 783 784 // the point is on the edge 785 // if the point is one the boundary 786 // add the point in outside part 787 if ( t->det >=0) { // inside triangle 788 if ((( Triangle *) ta)->det < 0 ) { 789 // add in outside triangle 790 Add(s,( Triangle *) ta); 791 return;} 792 }} 793 else { 794 printf("bug (%i): Bug double points in\n",nbd0); 795 throw ErrorException(__FUNCT__,exprintf("See above")); 796 } 797 798 // remove de MarkUnSwap edge 799 t->SetUnMarkUnSwap(0); 800 t->SetUnMarkUnSwap(1); 801 t->SetUnMarkUnSwap(2); 802 803 tt[0]= t; 804 tt[1]= &triangles[nbt++]; 805 tt[2]= &triangles[nbt++]; 806 807 if (nbt>nbtx) { 808 throw ErrorException(__FUNCT__,exprintf("Not ebough triangles")); 809 } 810 811 *tt[1]= *tt[2]= *t; 812 // gestion of the link 813 tt[0]->link=tt[1]; 814 tt[1]->link=tt[2]; 815 816 (* tt[0])(OppositeVertex[0])=&s; 817 (* tt[1])(OppositeVertex[1])=&s; 818 (* tt[2])(OppositeVertex[2])=&s; 819 820 tt[0]->det=det3[0]; 821 tt[1]->det=det3[1]; 822 tt[2]->det=det3[2]; 823 824 // update adj des triangles externe 825 tt[0]->SetAdjAdj(0); 826 tt[1]->SetAdjAdj(1); 827 tt[2]->SetAdjAdj(2); 828 // update des adj des 3 triangle interne 829 const int i0 = 0; 830 const int i1= NextEdge[i0]; 831 const int i2 = PreviousEdge[i0]; 832 833 tt[i0]->SetAdj2(i2,tt[i2],i0); 834 tt[i1]->SetAdj2(i0,tt[i0],i1); 835 tt[i2]->SetAdj2(i1,tt[i1],i2); 836 837 tt[0]->SetTriangleContainingTheVertex(); 838 tt[1]->SetTriangleContainingTheVertex(); 839 tt[2]->SetTriangleContainingTheVertex(); 840 841 842 // swap if the point s is on a edge 843 if(izerodet>=0) { 844 int rswap =tt[izerodet]->swap(iedge); 845 846 if (!rswap) { 847 throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge")); 782 det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; 783 det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; 784 det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ; 785 } 786 } 787 788 if (!det3[0]) izerodet=0,nbzerodet++; 789 if (!det3[1]) izerodet=1,nbzerodet++; 790 if (!det3[2]) izerodet=2,nbzerodet++; 791 792 //if nbzerodet>0, point s is on an egde or on a vertex 793 if (nbzerodet >0 ){ 794 if (nbzerodet == 1) { 795 iedge = OppositeEdge[izerodet]; 796 TriangleAdjacent ta = t->Adj(iedge); 797 798 // the point is on the edge 799 // if the point is one the boundary 800 // add the point in outside part 801 if ( t->det >=0) { // inside triangle 802 if ((( Triangle *) ta)->det < 0 ) { 803 // add in outside triangle 804 Add(s,( Triangle *) ta); 805 return; 848 806 } 849 807 } 808 } 809 else { 810 t->Echo(); 811 printf("\nproblem while trying to add:\n"); 812 s.Echo(); 813 throw ErrorException(__FUNCT__,exprintf("Bug in Triangles::Add points duplicated %i times",nbzerodet)); 814 } 815 } 816 817 // remove de MarkUnSwap edge 818 t->SetUnMarkUnSwap(0); 819 t->SetUnMarkUnSwap(1); 820 t->SetUnMarkUnSwap(2); 821 822 tt[0]= t; 823 tt[1]= &triangles[nbt++]; 824 tt[2]= &triangles[nbt++]; 825 826 if (nbt>nbtx) { 827 throw ErrorException(__FUNCT__,exprintf("Not ebough triangles")); 828 } 829 830 *tt[1]= *tt[2]= *t; 831 // gestion of the link 832 tt[0]->link=tt[1]; 833 tt[1]->link=tt[2]; 834 835 (* tt[0])(OppositeVertex[0])=&s; 836 (* tt[1])(OppositeVertex[1])=&s; 837 (* tt[2])(OppositeVertex[2])=&s; 838 839 tt[0]->det=det3[0]; 840 tt[1]->det=det3[1]; 841 tt[2]->det=det3[2]; 842 843 // update adj des triangles externe 844 tt[0]->SetAdjAdj(0); 845 tt[1]->SetAdjAdj(1); 846 tt[2]->SetAdjAdj(2); 847 // update des adj des 3 triangle interne 848 const int i0 = 0; 849 const int i1= NextEdge[i0]; 850 const int i2 = PreviousEdge[i0]; 851 852 tt[i0]->SetAdj2(i2,tt[i2],i0); 853 tt[i1]->SetAdj2(i0,tt[i0],i1); 854 tt[i2]->SetAdj2(i1,tt[i1],i2); 855 856 tt[0]->SetTriangleContainingTheVertex(); 857 tt[1]->SetTriangleContainingTheVertex(); 858 tt[2]->SetTriangleContainingTheVertex(); 859 860 861 // swap if the point s is on a edge 862 if(izerodet>=0) { 863 int rswap =tt[izerodet]->swap(iedge); 864 865 if (!rswap) { 866 throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge")); 867 } 868 } 850 869 } 851 870 /*}}}1*/ … … 2011 2030 /*}}}1*/ 2012 2031 /*FUNCTION Triangles::ForceBoundary{{{1*/ 2013 void Triangles::ForceBoundary() { 2014 long int verbosity=2; 2015 if (verbosity > 2) printf(" ForceBoundary nb of edge: %i\n",nbe); 2016 int k=0; 2017 Int4 nbfe=0,nbswp=0,Nbswap=0; 2018 for (Int4 t = 0; t < nbt; t++){ 2019 if (!triangles[t].det) k++; 2020 } 2021 if (k!=0) { 2022 throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k)); 2023 } 2024 2032 void Triangles::ForceBoundary() { 2033 long int verbosity=2; 2034 int k=0; 2035 int nbfe=0,nbswp=0,Nbswap=0; 2036 2037 //display 2038 if (verbosity > 2) printf(" ForceBoundary nb of edge: %i\n",nbe); 2039 2040 //check that there is no triangle with 0 determinant 2041 for (Int4 t = 0; t < nbt; t++){ 2042 if (!triangles[t].det) k++; 2043 } 2044 if (k!=0) { 2045 throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k)); 2046 } 2047 2048 //Force Edges 2025 2049 TriangleAdjacent ta(0,0); 2026 for (Int4 i = 0; i < nbe; i++) 2027 { 2050 for (Int4 i = 0; i < nbe; i++){ 2028 2051 nbswp = ForceEdge(edges[i][0],edges[i][1],ta); 2029 2052 … … 2033 2056 if ( nbswp < 0 && k < 5){ 2034 2057 throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1]))); 2035 2058 } 2036 2059 if ( nbswp >=0 && edges[i].onGeometry->Cracked()) 2037 2060 ta.SetCracked(); 2038 } 2039 2061 } 2040 2062 2041 2063 if (k!=0) { 2042 2064 throw ErrorException(__FUNCT__,exprintf("There are %i lost edges, the boundary might be crossing",k)); 2043 2065 } 2044 for (Int4 j=0;j<nbv;j++) 2045 Nbswap += vertices[j].Optim(1,0); 2066 for (Int4 j=0;j<nbv;j++){ 2067 Nbswap += vertices[j].Optim(1,0); 2068 } 2046 2069 if (verbosity > 3) printf(" number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 2047 }2070 } 2048 2071 /*}}}1*/ 2049 2072 /*FUNCTION Triangles::FillHoleInMesh{{{1*/ … … 2209 2232 Vertex *vi = ordre[icount]; 2210 2233 Icoor2 dete[3]; 2211 Triangle *tcvi = FindTriangleCont ening(vi->i,dete);2234 Triangle *tcvi = FindTriangleContaining(vi->i,dete); 2212 2235 quadtree->Add(*vi); 2213 2236 Add(*vi,tcvi,dete); … … 2584 2607 } 2585 2608 /*}}}1*/ 2586 /*FUNCTION Triangles::FindTriangleCont ening{{{1*/2587 Triangle * Triangles::FindTriangleCont ening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const {2609 /*FUNCTION Triangles::FindTriangleContaining{{{1*/ 2610 Triangle * Triangles::FindTriangleContaining(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { 2588 2611 Triangle * t=0; 2589 2612 int j,jp,jn,jj; … … 2599 2622 printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a)); 2600 2623 } 2601 throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleCont ening"));2624 throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContaining")); 2602 2625 } 2603 2626 if (a<vertices || a>=vertices+nbv){ … … 2682 2705 Gh.NbRef++;// add a ref to GH 2683 2706 2684 Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve; 2707 int i,j,k; 2708 int NbOfCurves=0,NbNewPoints,NbEdgeCurve; 2685 2709 Real8 lcurve,lstep,s; 2710 const int MaxSubEdge = 10; 2686 2711 2687 2712 R2 AB; … … 2710 2735 } 2711 2736 2712 // compute nbv (number of vertices)2737 //Add all the geometrical vertices to the mesh 2713 2738 nbv=0; 2714 2739 for (i=0;i<Gh.nbv;i++){ … … 2730 2755 Gh[i].to=vertices+nbv; 2731 2756 2732 2733 2757 //Build VerticesOnGeomVertex for current point 2734 2758 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); … … 2753 2777 nbe=0; 2754 2778 Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 2755 2756 2757 2779 Gh.UnMarkEdges(); 2758 2780 NbOfCurves = 0; … … 2761 2783 for (i=0;i<Gh.nbe;i++) { 2762 2784 2785 //ei = current Geometrical edge 2763 2786 GeometricalEdge &ei=Gh.edges[i]; 2764 2787 … … 2773 2796 Edge* PreviousNewEdge=NULL; 2774 2797 2775 lstep = -1;//do not create points 2798 //do not create points if required 2799 lstep = -1; 2776 2800 if(ei.Required()){ 2777 2801 if (j==0){ … … 2781 2805 a=ei(0)->The(); 2782 2806 b=ei(1)->The(); 2807 2808 //check that edges has been allocated 2783 2809 if (!edges){ 2784 throw ErrorException(__FUNCT__,exprintf(" !edges"));2810 throw ErrorException(__FUNCT__,exprintf("edges has not been allocated...")); 2785 2811 } 2786 2812 edges[nbe].v[0]=a->to; … … 2794 2820 } 2795 2821 } 2796 else { // on curve ------ 2822 2823 //if not required: on curve 2824 else { 2797 2825 for ( int kstep=0;kstep<= step;kstep++){ 2798 // if 2nd step where 2 step 2826 //step=0, nohing 2827 //step=1, 2828 //step=2 2799 2829 // 1 compute the length of the curve 2800 2830 // 2 create the points and edge … … 2807 2837 lcurve =0; 2808 2838 s = lstep; 2809 int k=j; 2810 e = & ei; 2811 a=ei(k)->The(); 2812 va = a->to; 2813 e->SetMark(); 2814 2815 // if SameGeo We have go in the background geometry 2816 // to find the discretisation of the curve 2839 2840 // i = edge number, j=[0;1] vertex number in edge 2841 2842 k=j; // k = vertex number in edge 2843 e=&ei; // e = address of current edge 2844 a=ei(k)->The(); // a = pointer toward the kth vertex of the current edge 2845 va = a->to; // va= pointer toward vertex of ??? 2846 e->SetMark(); // Mark e 2847 2848 //if SameGeo We have go to the background geometry 2849 //to find the discretisation of the curve 2817 2850 for(;;){ 2818 2851 k = 1-k; 2819 b= (*e)(k)->The(); 2820 AB = b->r - a->r; 2821 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; 2822 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; 2823 Real8 ledge = (MA(AB) + MB(AB))/2; 2824 2825 const int MaxSubEdge = 10; 2852 b = (*e)(k)->The(); // b = pointer toward the other vertex of the current edge 2853 AB = b->r - a->r; // AB = vector of the current edge 2854 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 2855 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 2856 Real8 ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 2857 2858 /* We are now creating the edges of the mesh from the 2859 * geometrical edge selected above. 2860 * The edge will be divided according to the metric 2861 * previously computed and cannot be divided more 2862 * than 10 times (MaxSubEdge). */ 2863 2864 //By default, there is only one subedge that is the geometrical edge itself 2826 2865 int NbSubEdge = 1; 2866 2867 //initialize lSubEdge, holding the length of each subedge (cannot be higher than 10) 2827 2868 Real8 lSubEdge[MaxSubEdge]; 2828 R2 A,B; 2869 2870 //Build Subedges according to the edge length 2871 //if ledge < 1.5 (between one and 2), take the edge as is 2829 2872 if (ledge < 1.5) lSubEdge[0] = ledge; 2873 //else, divide the edge 2830 2874 else { 2875 //compute number of subedges (division of the edge) 2831 2876 NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); 2832 A= a->r; 2833 Metric MAs =MA,MBs; 2834 ledge = 0; 2835 Real8 x =0, xstep= 1. / NbSubEdge; 2877 //A and B are the position of points on the edge 2878 R2 A,B; 2879 A=a->r; 2880 Metric MAs=MA,MBs; 2881 ledge=0; 2882 Real8 x =0, xstep= 1./NbSubEdge; 2836 2883 for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { 2837 2884 x += xstep; … … 2844 2891 2845 2892 Real8 lcurveb = lcurve+ ledge ; 2893 2894 //Now, create corresponding points 2895 //s = lstep 2846 2896 while (lcurve<=s && s <= lcurveb && nbv < nbvend){ 2847 // New points2848 2849 // Real8 aa=(lcurveb-s)/ledge;2850 // Real8 bb=(s-lcurve)/ledge;2851 2897 2852 2898 Real8 ss = s-lcurve; 2853 // 1) find the SubEdge containing ss by dichotomie 2899 2900 /*Find the SubEdge containing ss using Dichotomy*/ 2901 2854 2902 int kk0=-1,kk1=NbSubEdge-1,kkk; 2855 2903 Real8 ll0=0,ll1=ledge,llk; 2856 while (kk1-kk0>1) 2857 { 2904 while (kk1-kk0>1){ 2858 2905 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) 2859 2906 kk1=kkk,ll1=llk; 2860 2907 else 2861 kk0=kkk,ll0=llk; }2862 if (kk1 == kk0){2863 throw ErrorException(__FUNCT__,exprintf("kk1 == kk0"));2864 }2865 2866 Real8 sbb = (ss-ll0 )/(ll1-ll0); 2867 Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb;2868 2869 // new vertex on edge 2870 vb = &vertices[nbv++];2871 vb->m = Metric(aa,a->m,bb,b->m);2872 vb->ReferenceNumber = e->ref;2873 vb->DirOfSearch =NoDirOfSearch;2874 Real8 abcisse = k ? bb : aa;2875 vb->r = e->F( abcisse );2876 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse);2877 2878 // to take into account the direction of the edge 2879 s += lstep;2880 edges[nbe].v[0]=va;2881 edges[nbe].v[1]=vb;2882 edges[nbe].ref = e->ref;2883 edges[nbe].onGeometry = e;2884 edges[nbe].adj[0] = PreviousNewEdge;2885 if(PreviousNewEdge)2886 PreviousNewEdge->adj[1] =&edges[nbe];2887 PreviousNewEdge = edges +nbe;2888 2889 2908 kk0=kkk,ll0=llk; 2909 } 2910 if (kk1 == kk0){ 2911 throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); 2912 } 2913 2914 Real8 sbb = (ss-ll0 )/(ll1-ll0); 2915 Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb; 2916 2917 // new vertex on edge 2918 vb = &vertices[nbv++]; 2919 vb->m = Metric(aa,a->m,bb,b->m); 2920 vb->ReferenceNumber = e->ref; 2921 vb->DirOfSearch =NoDirOfSearch; 2922 Real8 abcisse = k ? bb : aa; 2923 vb->r = e->F( abcisse ); 2924 VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); 2925 2926 // to take into account the direction of the edge 2927 s += lstep; 2928 edges[nbe].v[0]=va; 2929 edges[nbe].v[1]=vb; 2930 edges[nbe].ref = e->ref; 2931 edges[nbe].onGeometry = e; 2932 edges[nbe].adj[0] = PreviousNewEdge; 2933 if(PreviousNewEdge) PreviousNewEdge->adj[1]=&edges[nbe]; 2934 PreviousNewEdge=edges+nbe; 2935 nbe++; 2936 va = vb; 2890 2937 } 2891 2938 lcurve = lcurveb; … … 2947 2994 } 2948 2995 2996 //Insert points inside existing triangles 2949 2997 Insert(); 2998 2999 //Force the boundary 2950 3000 ForceBoundary(); 3001 3002 //Extract SubDomains 2951 3003 FindSubDomain(); 2952 3004 … … 3369 3421 /*FUNCTION Triangles::Insert{{{1*/ 3370 3422 void Triangles::Insert() { 3423 /*Insert points in the existing Geometry*/ 3424 3425 //Intermediary 3426 int i; 3427 3428 /*Get options*/ 3371 3429 long int verbosity=2; 3430 3431 //Display info 3372 3432 if (verbosity>2) printf(" Insert initial %i vertices\n",nbv); 3433 3434 //Compute integer coordinates and determinants for the existing vertices (from Geometry) 3373 3435 SetIntCoor(); 3374 Int4 i; 3375 3376 //fill out ordre with the adresses of each vertex 3436 3437 /*Now we want to build a list (ordre) of the vertices in a random 3438 * order. To do so, we use the following method: 3439 * 3440 * From an initial k0 in [0 nbv[ random (vertex number) 3441 * the next k (vertex number) is computed using a big 3442 * prime number (PN>>nbv) following: 3443 * 3444 * k_{i+1} = k_i + PN [nbv] 3445 * 3446 * let's show that: 3447 * 3448 * for all j in [0 nbv[, ∃! i in [0 nbv[ such that k_i=j 3449 * 3450 * Let's assume that there are 2 distinct j1 and j2 such that 3451 * k_j1 = k_j2 3452 * 3453 * This means that 3454 * 3455 * k0+j1*PN = k0+j2*PN [nbv] 3456 * (j1-j2)*PN =0 [nbv] 3457 * since PN is a prime number larger than nbv, and nbv!=1 3458 * j1-j2=0 [nbv] 3459 * BUT 3460 * j1 and j2 are in [0 nbv[ which is impossible. 3461 * 3462 * We hence have built a random list of nbv elements of 3463 * [0 nbv[ all distincts*/ 3377 3464 for (i=0;i<nbv;i++) ordre[i]= &vertices[i] ; 3378 3379 //Build a random order3380 3465 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; 3381 Int4 k3 = rand()%nbv ; 3382 for (int is3=0; is3<nbv; is3++) 3383 ordre[is3]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 3384 3466 int k0=rand()%nbv; 3467 for (int is3=0; is3<nbv; is3++){ 3468 ordre[is3]= &vertices[k0=(k0+PrimeNumber)%nbv]; 3469 } 3470 3471 /*Modify ordre such that the first 3 vertices form a triangle*/ 3472 3473 //get first vertex i such that [0,1,i] are not aligned 3385 3474 for (i=2 ; det( ordre[0]->i, ordre[1]->i, ordre[i]->i ) == 0;){ 3475 //if i is higher than nbv, it means that all the determinants are 0, 3476 //all vertices are aligned! 3386 3477 if ( ++i >= nbv) { 3387 3478 throw ErrorException(__FUNCT__,exprintf("all the vertices are aligned")); … … 3389 3480 } 3390 3481 // exchange i et 2 in "ordre" so that 3391 // the first 3 vertices are not aligned 3482 // the first 3 vertices are not aligned (real triangle) 3392 3483 Exchange(ordre[2], ordre[i]); 3393 3484 3394 // We add a vertex at the infinity to build the mesh 3395 // so that we have a simple definition of the boundary edges 3485 /*Take the first edge formed by the first two vertices and build 3486 * the initial simple mesh from this edge and 2 boundary triangles*/ 3487 3488 Vertex * v0=ordre[0], *v1=ordre[1]; 3489 3396 3490 nbt = 2; 3397 3398 // We build a trivial mesh from one edge3399 // and 2 triangles built with the 2 oriented3400 // et3401 Vertex * v0=ordre[0], *v1=ordre[1];3402 3403 3491 triangles[0](0) = 0; //infinite vertex 3404 3492 triangles[0](1) = v0; 3405 3493 triangles[0](2) = v1; 3406 3407 3494 triangles[1](0) = 0;//infinite vertex 3408 3495 triangles[1](2) = v0; 3409 3496 triangles[1](1) = v1; 3497 3498 //Build adjacence 3410 3499 const int e0 = OppositeEdge[0]; 3411 3500 const int e1 = NextEdge[e0]; … … 3429 3518 quadtree->Add(*v1); 3430 3519 3431 //the vertices are added one by one 3520 /*Now, add the vertices One by One*/ 3521 3432 3522 Int4 NbSwap=0; 3433 3434 3435 3523 if (verbosity>3) printf(" Begining of insertion process...\n"); 3436 3524 3437 3525 for (Int4 icount=2; icount<nbv; icount++) { 3438 Vertex *vi=ordre[icount]; 3526 3527 //Get new vertex 3528 Vertex *newvertex=ordre[icount]; 3529 3530 //Find the triangle in which newvertex is located 3439 3531 Icoor2 dete[3]; 3440 Triangle *tcvi = FindTriangleContening(vi->i,dete); 3441 quadtree->Add(*vi); 3442 Add(*vi,tcvi,dete); 3443 NbSwap += vi->Optim(1,0); 3444 } 3532 Triangle* tcvi = FindTriangleContaining(newvertex->i,dete); //(newvertex->i = integer coordinates) 3533 3534 //Add newvertex to the quadtree 3535 quadtree->Add(*newvertex); 3536 3537 //Add newvertex to the existing mesh 3538 Add(*newvertex,tcvi,dete); 3539 3540 //Make the mesh Delaunay around newvertex by swaping the edges 3541 NbSwap += newvertex->Optim(1,0); 3542 } 3543 3544 //Display info 3445 3545 if (verbosity>3) { 3446 3546 printf(" NbSwap of insertion: %i\n",NbSwap); … … 3450 3550 } 3451 3551 NbUnSwap = 0; 3452 // construction d'un ordre aleatoire 3453 // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; 3552 3454 3553 #ifdef NBLOOPOPTIM 3455 3554 3456 k 3= rand()%nbv ;3555 k0 = rand()%nbv ; 3457 3556 for (int is4=0; is4<nbv; is4++) 3458 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 3459 3460 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 3461 { 3557 ordre[is4]= &vertices[k0 = (k0 + PrimeNumber)% nbv]; 3558 3559 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++){ 3462 3560 Int4 NbSwap = 0; 3463 3561 for (int is1=0; is1<nbv; is1++) … … 3471 3569 NbUnSwap = 0; 3472 3570 if(!NbSwap) break; 3473 3571 } 3474 3572 ReMakeTriangleContainingTheVertex(); 3475 3573 // because we break the TriangleContainingTheVertex … … 3528 3626 } 3529 3627 vj.ReferenceNumber=0; 3530 Triangle *tcvj=FindTriangleCont ening(vj.i,dete);3628 Triangle *tcvj=FindTriangleContaining(vj.i,dete); 3531 3629 if (tcvj && !tcvj->link){ 3532 3630 tcvj->Echo(); … … 3703 3801 I2 a = toI2(A); 3704 3802 Icoor2 deta[3]; 3705 Triangle * t =FindTriangleCont ening(a,deta);3803 Triangle * t =FindTriangleContaining(a,deta); 3706 3804 if (t->det <0) { // outside 3707 3805 double ba,bb; … … 4221 4319 /*FUNCTION Triangles::SetIntCoor{{{1*/ 4222 4320 void Triangles::SetIntCoor(const char * strfrom) { 4321 /*Set integer coordinate for existing vertices*/ 4322 4323 //Get extrema coordinates of the existing vertices 4223 4324 pmin = vertices[0].r; 4224 4325 pmax = vertices[0].r; 4225 4226 // recherche des extrema des vertices pmin,pmax4227 4326 Int4 i; 4228 4327 for (i=0;i<nbv;i++) { … … 4235 4334 pmin = pmin-DD; 4236 4335 pmax = pmax+DD; 4336 4337 //Compute coefIcoor 4237 4338 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 4238 4339 if (coefIcoor<=0){ 4239 throw ErrorException(__FUNCT__,exprintf("coefIcoor <=0"));4340 throw ErrorException(__FUNCT__,exprintf("coefIcoor should be positive, a problem in the geometry is likely")); 4240 4341 } 4241 4342 … … 4251 4352 Vertex & v1 = triangles[i][1]; 4252 4353 Vertex & v2 = triangles[i][2]; 4253 if ( &v0 && &v1 && &v2 ) // a good triangles; 4254 { 4354 4355 //If this is not a boundary triangle 4356 if ( &v0 && &v1 && &v2 ){ 4255 4357 triangles[i].det= det(v0,v1,v2); 4256 4358 if (triangles[i].det <=0 && Nberr++ <10){ … … 4263 4365 } 4264 4366 } 4265 } 4266 else triangles[i].det= -1; // Null triangle; 4367 } 4368 4369 //else, set as -1 4370 else triangles[i].det=-1; 4267 4371 } 4268 4372 } … … 5043 5147 vi.ReferenceNumber=0; 5044 5148 vi.DirOfSearch =NoDirOfSearch; 5045 Triangle *tcvi = FindTriangleCont ening(vi.i,dete);5149 Triangle *tcvi = FindTriangleContaining(vi.i,dete); 5046 5150 if (tcvi && !tcvi->link) { 5047 5151 printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n"); … … 5244 5348 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 5245 5349 Icoor2 dete[3]; 5246 Triangle & tb = *Th.FindTriangleCont ening(I,dete);5350 Triangle & tb = *Th.FindTriangleContaining(I,dete); 5247 5351 5248 5352 if (tb.link) -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r2944 r2965 55 55 I2 IBTh = BTh.toI2(PNew); 56 56 57 tstart=BTh.FindTriangleCont ening(IBTh,deta,tstart);57 tstart=BTh.FindTriangleContaining(IBTh,deta,tstart); 58 58 59 59 if (tstart->det <0){ // outside
Note:
See TracChangeset
for help on using the changeset viewer.
