![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
- Timestamp:
- 09/15/17 18:44:44 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp
r21759 r22085 133 133 nbt++; 134 134 } 135 if (nbt==0 && nbv==0) { 135 if (nbt==0 && nbv==0){ 136 delete [] refv; 136 137 _error_("All triangles have been removed"); 137 138 } … … 698 699 for (i=0;i<nbsubdomains;i++){ 699 700 bamgmesh->SubDomains[i*4+0]=3; 700 bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)] ;701 bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing 701 702 bamgmesh->SubDomains[i*4+2]=1; 702 703 bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber; … … 980 981 981 982 /*Methods*/ 982 void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/*{{{*/983 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectGeomMetric)*/984 985 /*Get options*/986 double anisomax = bamgopts->anisomax;987 double errg = bamgopts->errg;988 989 double ss[2]={0.00001,0.99999};990 double errC = 2*sqrt(2*errg);991 double hmax = Gh.MaximalHmax();992 double hmin = Gh.MinimalHmin();993 994 //check that hmax is positive995 _assert_(hmax>0);996 997 //errC cannot be higher than 1998 if(errC>1) errC=1;999 1000 //Set all vertices to "on"1001 SetVertexFieldOn();1002 1003 //loop over all the vertices on edges1004 for(int i=0;i<nbe;i++){1005 for (int j=0;j<2;j++){1006 1007 BamgVertex V;1008 VertexOnGeom GV;1009 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);1010 1011 GeomEdge* eg = GV;1012 double s = GV;1013 R2 tg;1014 double R1= eg->R1tg(s,tg);1015 double ht=hmax;1016 // err relative to the length of the edge1017 if (R1>1.0e-20) {1018 ht = Min(Max(errC/R1,hmin),hmax);1019 }1020 double hn=Min(hmax,ht*anisomax);1021 1022 if (ht<=0 || hn<=0){1023 _error_("ht<=0 || hn<=0");1024 }1025 EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg);1026 Metric MVp(Vp);1027 edges[i][j].m.IntersectWith(MVp);1028 }1029 }1030 // the problem is for the vertex on vertex1031 }1032 /*}}}*/1033 983 void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ 1034 984 // Hessiantype = 0 => H is computed using double L2 projection … … 1049 999 } 1050 1000 /*}}}*/ 1051 void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3){/*{{{*/1001 void Mesh::AddVertex( BamgVertex &s,Triangle* t, long long* det3){/*{{{*/ 1052 1002 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/ 1053 1003 // ------------------------------- … … 1068 1018 /*Intermediaries*/ 1069 1019 Triangle* tt[3]; //the three triangles 1070 Icoor2det3local[3]; //three determinants (integer)1020 long long det3local[3]; //three determinants (integer) 1071 1021 int nbzerodet =0; //number of zeros in det3 1072 1022 int izerodet=-1; //egde containing the vertex s … … 1079 1029 1080 1030 //determinant of t 1081 Icoor2detOld=t->det;1031 long long detOld=t->det; 1082 1032 1083 1033 /* infvertexindex = index of the infinite vertex (NULL) … … 1231 1181 int i,j,k,kk,it,jt; 1232 1182 int verbose=0; 1233 double cutoffradian=10*Pi/180;1234 1183 1235 1184 /*Recover options*/ 1236 if 1185 if(bamgopts){ 1237 1186 verbose=bamgopts->verbose; 1238 cutoffradian=bamgopts->MaxCornerAngle*Pi/180;1239 1187 } 1240 1188 … … 1243 1191 1244 1192 //check that the mesh is not empty 1245 if (nbt<=0 || nbv <=0 ) { 1246 _error_("nbt or nbv is negative (Mesh empty?)"); 1247 } 1248 1249 //Gh is the geometry of the mesh (this), initialize MaxCornerAngle 1250 if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian; 1193 if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); 1251 1194 1252 1195 /*Construction of the edges*/ … … 1447 1390 1448 1391 //check that nbsubdomains is empty 1449 if (nbsubdomains){ 1450 _error_("nbsubdomains should be 0"); 1451 } 1392 if(nbsubdomains) _error_("nbsubdomains should be 0"); 1452 1393 nbsubdomains=0; 1453 1394 … … 1653 1594 long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]); 1654 1595 k = edge4->SortAndFind(i0,i1); 1655 _assert_(k>=0); 1596 if(k<0){ 1597 delete [] colorV; 1598 _error_("This should not happen"); 1599 } 1656 1600 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; 1657 1601 subdomains[i].edge = edges+k; … … 2113 2057 else if (xy==1) dd=dxdy; 2114 2058 else if (xy==2) dd=dydy; 2115 else _error_("not supported yet"); 2059 else{ 2060 delete [] detT; 2061 delete [] Mmass; 2062 delete [] workT; 2063 delete [] workV; 2064 delete [] Mmassxx; 2065 delete [] OnBoundary; 2066 _error_("not supported yet"); 2067 } 2116 2068 // do leat 2 iteration for boundary problem 2117 2069 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ … … 2162 2114 delete [] dxdy; 2163 2115 delete [] dydy; 2164 delete [] 2116 delete [] workT; 2165 2117 delete [] workV; 2166 2118 delete [] Mmassxx; 2167 delete [] 2119 delete [] OnBoundary; 2168 2120 2169 2121 } … … 2453 2405 if (nbt == nbtout || !NbSubDomTot) { 2454 2406 delete [] HeapArete; 2407 delete [] HeapTriangle; 2455 2408 _error_("The boundary is not close: all triangles are outside"); 2456 2409 } … … 2649 2602 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2650 2603 2651 /* initialize random seed: */2652 srand(19999999);2653 2654 2604 /*Initialize fields*/ 2655 NbRef=0; 2656 quadtree=NULL; 2657 nbv=0; 2658 nbt=0; 2659 nbe=0; 2660 edges=NULL; 2661 nbsubdomains=0; 2662 subdomains=NULL; 2663 maxnbv=maxnbv_in; 2664 maxnbt=2 *maxnbv_in-2; 2665 NbVertexOnBThVertex=0; 2666 VertexOnBThVertex=NULL; 2667 NbVertexOnBThEdge=0; 2668 VertexOnBThEdge=NULL; 2669 NbCrackedVertices=0; 2670 CrackedVertices =NULL; 2671 NbCrackedEdges =0; 2672 CrackedEdges =NULL; 2673 NbVerticesOnGeomVertex=0; 2674 VerticesOnGeomVertex=NULL; 2675 NbVerticesOnGeomEdge=0; 2676 VerticesOnGeomEdge=NULL; 2605 this->NbRef = 0; 2606 this->quadtree = NULL; 2607 this->nbv = 0; 2608 this->nbt = 0; 2609 this->nbe = 0; 2610 this->edges = NULL; 2611 this->nbsubdomains = 0; 2612 this->subdomains = NULL; 2613 this->maxnbv = maxnbv_in; 2614 this->maxnbt = 2 *maxnbv_in-2; 2615 this->NbVertexOnBThVertex = 0; 2616 this->VertexOnBThVertex = NULL; 2617 this->NbVertexOnBThEdge = 0; 2618 this->VertexOnBThEdge = NULL; 2619 this->NbCrackedVertices = 0; 2620 this->CrackedVertices = NULL; 2621 this->NbCrackedEdges = 0; 2622 this->CrackedEdges = NULL; 2623 this->NbVerticesOnGeomVertex = 0; 2624 this->VerticesOnGeomVertex = NULL; 2625 this->NbVerticesOnGeomEdge = 0; 2626 this->VerticesOnGeomEdge = NULL; 2627 2628 /*Initialize random seed*/ 2629 this->randomseed = 1; 2677 2630 2678 2631 /*Allocate if maxnbv_in>0*/ 2679 2632 if(maxnbv_in){ 2680 vertices=new BamgVertex[maxnbv];2681 _assert_(vertices);2682 orderedvertices=new BamgVertex* [maxnbv];2683 _assert_(orderedvertices);2684 triangles=new Triangle[maxnbt];2685 _assert_(t riangles);2633 this->vertices=new BamgVertex[this->maxnbv]; 2634 this->orderedvertices=new BamgVertex* [this->maxnbv]; 2635 this->triangles=new Triangle[this->maxnbt]; 2636 _assert_(this->vertices); 2637 _assert_(this->orderedvertices); 2638 _assert_(this->triangles); 2686 2639 } 2687 2640 else{ 2688 vertices=NULL;2689 orderedvertices=NULL;2690 t riangles=NULL;2691 maxnbt=0;2641 this->vertices = NULL; 2642 this->orderedvertices = NULL; 2643 this->triangles = NULL; 2644 this->maxnbt = 0; 2692 2645 } 2693 2646 } 2694 2647 /*}}}*/ 2695 void Mesh::Insert( bool random){/*{{{*/2648 void Mesh::Insert(){/*{{{*/ 2696 2649 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2697 2650 … … 2740 2693 //Get Prime number 2741 2694 const long PrimeNumber= BigPrimeNumber(nbv); 2742 int temp = rand(); if(!random) temp = 756804110; 2743 int k0=temp%nbv; 2695 int k0=this->RandomNumber(nbv); 2744 2696 2745 2697 //Build orderedvertices … … 2805 2757 2806 2758 //Find the triangle in which newvertex is located 2807 Icoor2det3[3];2759 long long det3[3]; 2808 2760 Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates) 2809 2761 … … 2825 2777 } 2826 2778 /*}}}*/ 2827 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap ,bool random) {/*{{{*/2779 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ 2828 2780 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2829 2781 … … 2832 2784 long i; 2833 2785 long NbSwap=0; 2834 Icoor2det3[3];2786 long long det3[3]; 2835 2787 2836 2788 //number of new points … … 2845 2797 /*construction of a random order*/ 2846 2798 const long PrimeNumber= BigPrimeNumber(nbv) ; 2847 //remainder of the division of rand() by nbvnew 2848 int temp = rand(); if(!random) temp = 1098566905; 2849 long k3 = temp%nbvnew; 2799 long k3 = long(this->RandomNumber(nbvnew)); 2850 2800 //loop over the new points 2851 2801 for (int is3=0; is3<nbvnew; is3++){ … … 2945 2895 } 2946 2896 } 2947 if(kk) _error_("See above"); 2897 if(kk){ 2898 delete [] e; 2899 _error_("See above"); 2900 } 2948 2901 2949 2902 return e; … … 3013 2966 } 3014 2967 /*}}}*/ 3015 Metric Mesh::MetricAt(const R2 & A) const{ /*{{{*/2968 Metric Mesh::MetricAt(const R2 & A){ /*{{{*/ 3016 2969 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/ 3017 2970 3018 2971 I2 a = R2ToI2(A); 3019 Icoor2deta[3];2972 long long deta[3]; 3020 2973 Triangle * t =TriangleFindFromCoord(a,deta); 3021 2974 if (t->det <0) { // outside … … 3086 3039 } 3087 3040 } 3088 if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");3041 //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); 3089 3042 Bh.CreateSingleVertexToTriangleConnectivity(); 3090 InsertNewPoints(nbvold,NbTSwap ,bamgopts->random);3043 InsertNewPoints(nbvold,NbTSwap); 3091 3044 } 3092 3045 else Bh.CreateSingleVertexToTriangleConnectivity(); … … 3146 3099 }// for triangle 3147 3100 3148 if (!InsertNewPoints(nbvold,NbTSwap ,bamgopts->random)) break;3101 if (!InsertNewPoints(nbvold,NbTSwap)) break; 3149 3102 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3150 3103 Headt = nbt; // empty list … … 3476 3429 for (int icount=2; icount<nbvb; icount++) { 3477 3430 BamgVertex *vi = orderedvertices[icount]; 3478 Icoor2det3[3];3431 long long det3[3]; 3479 3432 Triangle *tcvi = TriangleFindFromCoord(vi->i,det3); 3480 3433 quadtree->Add(*vi); … … 3893 3846 long iv = nbvold; 3894 3847 long NbSwap = 0; 3895 Icoor2det3[3];3848 long long det3[3]; 3896 3849 for (int i=nbvold;i<nbv;i++) {// for all the new point 3897 3850 BamgVertex & vi = vertices[i]; … … 3934 3887 } 3935 3888 /*}}}*/ 3936 Triangle * Mesh::TriangleFindFromCoord(const I2 & B, Icoor2 det3[3], Triangle *tstart) const{/*{{{*/3889 Triangle * Mesh::TriangleFindFromCoord(const I2 & B,long long det3[3], Triangle *tstart){/*{{{*/ 3937 3890 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/ 3938 3891 … … 3963 3916 } 3964 3917 3965 Icoor2detop ;3918 long long detop ; 3966 3919 3967 3920 /*initialize number of test triangle*/ … … 4009 3962 4010 3963 if (k==0) break; 4011 if (k==2 && BinaryRand()) Exchange(ii[0],ii[1]);3964 if (k==2 && this->RandomNumber(1)) Exchange(ii[0],ii[1]); 4012 3965 _assert_(k<3); 4013 3966 AdjacentTriangle t1 = t->Adj(jj=ii[0]); … … 4095 4048 4096 4049 /*Insert Vertices*/ 4097 Insert( true);4050 Insert(); 4098 4051 } 4099 4052 /*}}}*/ … … 4396 4349 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4397 4350 if (verbose>3) _printf_(" Inserting boundary points\n"); 4398 Insert( bamgopts->random);4351 Insert(); 4399 4352 4400 4353 //Force the boundary … … 4460 4413 printf("\n"); 4461 4414 if(NbVerticesOnGeomVertex >= maxnbv){ 4415 delete [] bcurve; 4462 4416 _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv); 4463 4417 } … … 4722 4676 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4723 4677 if (verbose>3) _printf_(" Inserting boundary points\n"); 4724 Insert( bamgopts->random);4678 Insert(); 4725 4679 4726 4680 //Force the boundary … … 4735 4689 NewPoints(BTh,bamgopts,KeepVertices) ; 4736 4690 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); 4691 } 4692 /*}}}*/ 4693 int Mesh::RandomNumber(int max){/*{{{*/ 4694 /* Generate a random number from 0 to max-1 using linear congruential 4695 * random number generator*/ 4696 4697 this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l; 4698 return this->randomseed / (714025l / max + 1); 4699 } /*}}}*/ 4700 int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4701 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4702 4703 int NbSwap =0; 4704 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4705 _error_("!a.t || !b.t"); 4706 } 4707 int k=0; 4708 taret=AdjacentTriangle(0,0); // erreur 4709 4710 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4711 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4712 // we turn around a in the direct direction 4713 4714 long long det2 = v2 ? det(*v2,a,b): -1 , det1; 4715 if(v2) // normal case 4716 det2 = det(*v2,a,b); 4717 else { // no chance infini vertex try the next 4718 tta= Previous(Adj(tta)); 4719 v2 = tta.EdgeVertex(0); 4720 vbegin =v2; 4721 if (!v2){ 4722 _error_("!v2"); 4723 } 4724 det2 = det(*v2,a,b); 4725 } 4726 4727 while (v2 != &b) { 4728 AdjacentTriangle tc = Previous(Adj(tta)); 4729 v1 = v2; 4730 v2 = tc.EdgeVertex(0); 4731 det1 = det2; 4732 det2 = v2 ? det(*v2,a,b): det2; 4733 4734 if((det1 < 0) && (det2 >0)) { 4735 // try to force the edge 4736 BamgVertex * va = &a, *vb = &b; 4737 tc = Previous(tc); 4738 if (!v1 || !v2){ 4739 _error_("!v1 || !v2"); 4740 } 4741 long long detss = 0,l=0; 4742 while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4743 if(l++ > 10000000) { 4744 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4745 } 4746 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4747 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4748 tc.SetLock(); 4749 a.Optim(1,0); 4750 b.Optim(1,0); 4751 taret = tc; 4752 return NbSwap; 4753 } 4754 else 4755 { 4756 taret = tc; 4757 return -2; // error boundary is crossing 4758 } 4759 } 4760 tta = tc; 4761 k++; 4762 if (k>=2000){ 4763 _error_("k>=2000"); 4764 } 4765 if ( vbegin == v2 ) return -1;// error 4766 } 4767 4768 tta.SetLock(); 4769 taret=tta; 4770 a.Optim(1,0); 4771 b.Optim(1,0); 4772 return NbSwap; 4773 } 4774 /*}}}*/ 4775 int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ 4776 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ 4777 // l'arete ta coupe l'arete pva pvb 4778 // de cas apres le swap sa coupe toujours 4779 // on cherche l'arete suivante 4780 // on suppose que detsa >0 et detsb <0 4781 // attention la routine echange pva et pvb 4782 4783 if(tt1.Locked()) return 0; // frontiere croise 4784 4785 AdjacentTriangle tt2 = Adj(tt1); 4786 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 4787 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 4788 if ( a1<0 || a1>=3 ){ 4789 _error_("a1<0 || a1>=3"); 4790 } 4791 4792 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 4793 BamgVertex & s1= (*t1)[OppositeVertex[a1]]; 4794 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 4795 4796 long long dets2 = det(*pva,*pvb,s2); 4797 long long det1=t1->det , det2=t2->det ; 4798 long long detT = det1+det2; 4799 if ((det1<=0 ) || (det2<=0)){ 4800 _error_("(det1<=0 ) || (det2<=0)"); 4801 } 4802 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 4803 _error_("(detsa>=0) || (detsb<=0)"); 4804 } 4805 long long ndet1 = bamg::det(s1,sa,s2); 4806 long long ndet2 = detT - ndet1; 4807 4808 int ToSwap =0; //pas de swap 4809 if ((ndet1 >0) && (ndet2 >0)) 4810 { // on peut swaper 4811 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 4812 ToSwap =1; 4813 else // swap alleatoire 4814 if (this->RandomNumber(1)) 4815 ToSwap =2; 4816 } 4817 if (ToSwap) NbSwap++, 4818 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 4819 4820 int ret=1; 4821 4822 if (dets2 < 0) {// haut 4823 dets1 = ToSwap ? dets1 : detsa ; 4824 detsa = dets2; 4825 tt1 = Previous(tt2) ;} 4826 else if (dets2 > 0){// bas 4827 dets1 = ToSwap ? dets1 : detsb ; 4828 detsb = dets2; 4829 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 4830 if(!ToSwap) tt1 = Next(tt2); 4831 } 4832 else { // changement de direction 4833 ret = -1; 4834 Exchange(pva,pvb); 4835 Exchange(detsa,detsb); 4836 Exchange(dets1,dets2); 4837 Exchange(tt1,tt2); 4838 dets1=-dets1; 4839 dets2=-dets2; 4840 detsa=-detsa; 4841 detsb=-detsb; 4842 4843 if(ToSwap){ 4844 if (dets2 < 0) {// haut 4845 dets1 = (ToSwap ? dets1 : detsa) ; 4846 detsa = dets2; 4847 tt1 = Previous(tt2) ;} 4848 else if(dets2 > 0){// bas 4849 dets1 = (ToSwap ? dets1 : detsb) ; 4850 detsb = dets2; 4851 if(!ToSwap) tt1 = Next(tt2); 4852 } 4853 else {// on a fin ??? 4854 tt1 = Next(tt2); 4855 ret =0;} 4856 } 4857 4858 } 4859 return ret; 4737 4860 } 4738 4861 /*}}}*/ … … 4748 4871 } 4749 4872 int kkk=0; 4750 Icoor2IJ_IA,IJ_AJ;4873 long long IJ_IA,IJ_AJ; 4751 4874 AdjacentTriangle edge(t,OppositeEdge[k]); 4752 4875 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { … … 4778 4901 } 4779 4902 /*}}}*/ 4780 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4781 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4782 4783 int NbSwap =0; 4784 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4785 _error_("!a.t || !b.t"); 4786 } 4787 int k=0; 4788 taret=AdjacentTriangle(0,0); // erreur 4789 4790 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4791 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4792 // we turn around a in the direct direction 4793 4794 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 4795 if(v2) // normal case 4796 det2 = det(*v2,a,b); 4797 else { // no chance infini vertex try the next 4798 tta= Previous(Adj(tta)); 4799 v2 = tta.EdgeVertex(0); 4800 vbegin =v2; 4801 if (!v2){ 4802 _error_("!v2"); 4803 } 4804 det2 = det(*v2,a,b); 4805 } 4806 4807 while (v2 != &b) { 4808 AdjacentTriangle tc = Previous(Adj(tta)); 4809 v1 = v2; 4810 v2 = tc.EdgeVertex(0); 4811 det1 = det2; 4812 det2 = v2 ? det(*v2,a,b): det2; 4813 4814 if((det1 < 0) && (det2 >0)) { 4815 // try to force the edge 4816 BamgVertex * va = &a, *vb = &b; 4817 tc = Previous(tc); 4818 if (!v1 || !v2){ 4819 _error_("!v1 || !v2"); 4820 } 4821 Icoor2 detss = 0,l=0; 4822 while ((SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4823 if(l++ > 10000000) { 4824 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4825 } 4826 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4827 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4828 tc.SetLock(); 4829 a.Optim(1,0); 4830 b.Optim(1,0); 4831 taret = tc; 4832 return NbSwap; 4833 } 4834 else 4835 { 4836 taret = tc; 4837 return -2; // error boundary is crossing 4838 } 4839 } 4840 tta = tc; 4841 k++; 4842 if (k>=2000){ 4843 _error_("k>=2000"); 4844 } 4845 if ( vbegin == v2 ) return -1;// error 4846 } 4847 4848 tta.SetLock(); 4849 taret=tta; 4850 a.Optim(1,0); 4851 b.Optim(1,0); 4852 return NbSwap; 4853 } 4854 /*}}}*/ 4855 void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){ /*{{{*/ 4903 void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,long long det1,long long det2){ /*{{{*/ 4856 4904 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/ 4857 4905 // -------------------------------------------------------------- … … 4897 4945 } // end swap 4898 4946 /*}}}*/ 4899 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) {/*{{{*/4900 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/4901 // l'arete ta coupe l'arete pva pvb4902 // de cas apres le swap sa coupe toujours4903 // on cherche l'arete suivante4904 // on suppose que detsa >0 et detsb <04905 // attention la routine echange pva et pvb4906 4907 if(tt1.Locked()) return 0; // frontiere croise4908 4909 AdjacentTriangle tt2 = Adj(tt1);4910 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent4911 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles4912 if ( a1<0 || a1>=3 ){4913 _error_("a1<0 || a1>=3");4914 }4915 4916 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];4917 BamgVertex & s1= (*t1)[OppositeVertex[a1]];4918 BamgVertex & s2= (*t2)[OppositeVertex[a2]];4919 4920 Icoor2 dets2 = det(*pva,*pvb,s2);4921 Icoor2 det1=t1->det , det2=t2->det ;4922 Icoor2 detT = det1+det2;4923 if ((det1<=0 ) || (det2<=0)){4924 _error_("(det1<=0 ) || (det2<=0)");4925 }4926 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb4927 _error_("(detsa>=0) || (detsb<=0)");4928 }4929 Icoor2 ndet1 = bamg::det(s1,sa,s2);4930 Icoor2 ndet2 = detT - ndet1;4931 4932 int ToSwap =0; //pas de swap4933 if ((ndet1 >0) && (ndet2 >0))4934 { // on peut swaper4935 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))4936 ToSwap =1;4937 else // swap alleatoire4938 if (BinaryRand())4939 ToSwap =2;4940 }4941 if (ToSwap) NbSwap++,4942 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);4943 4944 int ret=1;4945 4946 if (dets2 < 0) {// haut4947 dets1 = ToSwap ? dets1 : detsa ;4948 detsa = dets2;4949 tt1 = Previous(tt2) ;}4950 else if (dets2 > 0){// bas4951 dets1 = ToSwap ? dets1 : detsb ;4952 detsb = dets2;4953 //xxxx tt1 = ToSwap ? tt1 : Next(tt2);4954 if(!ToSwap) tt1 = Next(tt2);4955 }4956 else { // changement de direction4957 ret = -1;4958 Exchange(pva,pvb);4959 Exchange(detsa,detsb);4960 Exchange(dets1,dets2);4961 Exchange(tt1,tt2);4962 dets1=-dets1;4963 dets2=-dets2;4964 detsa=-detsa;4965 detsb=-detsb;4966 4967 if(ToSwap){4968 if (dets2 < 0) {// haut4969 dets1 = (ToSwap ? dets1 : detsa) ;4970 detsa = dets2;4971 tt1 = Previous(tt2) ;}4972 else if(dets2 > 0){// bas4973 dets1 = (ToSwap ? dets1 : detsb) ;4974 detsb = dets2;4975 if(!ToSwap) tt1 = Next(tt2);4976 }4977 else {// on a fin ???4978 tt1 = Next(tt2);4979 ret =0;}4980 }4981 4982 }4983 return ret;4984 }4985 /*}}}*/4986 4947 }
Note:
See TracChangeset
for help on using the changeset viewer.
