Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2847) +++ /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2848) @@ -18,206 +18,3 @@ int ForDebugging = 0; const Direction NoDirOfSearch = Direction(); - - Int4 AGoodNumberPrimeWith(Int4 n){ - const Int4 BigPrimeNumber[] ={ 567890359L, - 567890431L, 567890437L, 567890461L, 567890471L, - 567890483L, 567890489L, 567890497L, 567890507L, - 567890591L, 567890599L, 567890621L, 567890629L , 0}; - - Int4 o = 0; - Int4 pi = BigPrimeNumber[1]; - for (int i=0; BigPrimeNumber[i]; i++) { - Int4 r = BigPrimeNumber[i] % n; - Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); - if ( o < oo) - o=oo,pi=BigPrimeNumber[i];} - // cout << " AGoodNumberPrimeWith " << n << " " <Number(ta) << "::" - << CurrentTh->Number(ta.ns[0]) << "," - << CurrentTh->Number(ta.ns[1]) << "," - << CurrentTh->Number(ta.ns[2]) << "," - << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} " - << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} " - << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} " - << "]" ; - else - f << "[" - << ta.ns[0] << "," - << ta.ns[1] << "," - << ta.ns[2] << "," - << "{" << ta.at[0] << " " << ta.aa[0] << "} " - << "{" << ta.at[1] << " " << ta.aa[1] << "} " - << "{" << ta.at[2] << " " << ta.aa[2] << "} " - << "]" ; - return f;} - - - int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , - TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) - { // l'arete ta coupe l'arete pva pvb - // de cas apres le swap sa coupe toujours - // on cherche l'arete suivante - // on suppose que detsa >0 et detsb <0 - // attention la routine echange pva et pvb - - if(tt1.Locked()) return 0; // frontiere croise - - TriangleAdjacent tt2 = Adj(tt1); - Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent - Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles - assert ( a1 >= 0 && a1 < 3 ); - - Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; - Vertex & s1= (*t1)[OppositeVertex[a1]]; - Vertex & s2= (*t2)[OppositeVertex[a2]]; - - - Icoor2 dets2 = det(*pva,*pvb,s2); - Icoor2 det1=t1->det , det2=t2->det ; - Icoor2 detT = det1+det2; - assert((det1>0 ) && (det2 > 0)); - assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb - Icoor2 ndet1 = bamg::det(s1,sa,s2); - Icoor2 ndet2 = detT - ndet1; - - int ToSwap =0; //pas de swap - if ((ndet1 >0) && (ndet2 >0)) - { // on peut swaper - if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) - ToSwap =1; - else // swap alleatoire - if (BinaryRand()) - ToSwap =2; - } - if (ToSwap) NbSwap++, - bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); - - int ret=1; - - if (dets2 < 0) {// haut - dets1 = ToSwap ? dets1 : detsa ; - detsa = dets2; - tt1 = Previous(tt2) ;} - else if (dets2 > 0){// bas - dets1 = ToSwap ? dets1 : detsb ; - detsb = dets2; - //xxxx tt1 = ToSwap ? tt1 : Next(tt2); - if(!ToSwap) tt1 = Next(tt2); - } - else { // changement de sens - if (ForDebugging) cout << "changement de sens" << endl; - ret = -1; - Exchange(pva,pvb); - Exchange(detsa,detsb); - Exchange(dets1,dets2); - Exchange(tt1,tt2); - dets1=-dets1; - dets2=-dets2; - detsa=-detsa; - detsb=-detsb; - - if (ToSwap) - if (dets2 < 0) {// haut - dets1 = (ToSwap ? dets1 : detsa) ; - detsa = dets2; - tt1 = Previous(tt2) ;} - else if (dets2 > 0){// bas - dets1 = (ToSwap ? dets1 : detsb) ; - detsb = dets2; - if(!ToSwap) tt1 = Next(tt2); - } - else {// on a fin ??? - tt1 = Next(tt2); - ret =0;} - - } - return ret; - } - - int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) - { - int NbSwap =0; - assert(a.t && b.t); // the 2 vertex is in a mesh - int k=0; - taret=TriangleAdjacent(0,0); // erreur - - TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); - Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; - // we turn around a in the direct sens - - Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; - if(v2) // normal case - det2 = det(*v2,a,b); - else { // no chance infini vertex try the next - tta= Previous(Adj(tta)); - v2 = tta.EdgeVertex(0); - vbegin =v2; - assert(v2); - det2 = det(*v2,a,b); - // cout << " No Change try the next" << endl; - } - - while (v2 != &b) { - TriangleAdjacent tc = Previous(Adj(tta)); - v1 = v2; - v2 = tc.EdgeVertex(0); - det1 = det2; - det2 = v2 ? det(*v2,a,b): det2; - - if((det1 < 0) && (det2 >0)) { - // try to force the edge - Vertex * va = &a, *vb = &b; - tc = Previous(tc); - assert ( v1 && v2); - Icoor2 detss = 0,l=0,ks; - // cout << "Real ForcingEdge " << *va << *vb << detss << endl; - while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) - if(l++ > 10000000) { - cerr << " Loop in forcing Egde AB" - <<"\n vertex A " << a - <<"\n vertex B " << b - <<"\n nb de swap " << NbSwap - <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl; - - if ( CurrentTh ) - cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; - MeshError(990); - } - Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); - if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { - tc.SetLock(); - a.Optim(1,0); - b.Optim(1,0); - taret = tc; - return NbSwap; - } - else - { - taret = tc; - return -2; // error boundary is crossing - } - } - tta = tc; - assert(k++<2000); - if ( vbegin == v2 ) return -1;// error - } - - tta.SetLock(); - taret=tta; - a.Optim(1,0); - b.Optim(1,0); - return NbSwap; - } - } - Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2847) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 2848) @@ -6024,3 +6024,203 @@ /*}}}1*/ + Int4 AGoodNumberPrimeWith(Int4 n){ + const Int4 BigPrimeNumber[] ={ 567890359L, + 567890431L, 567890437L, 567890461L, 567890471L, + 567890483L, 567890489L, 567890497L, 567890507L, + 567890591L, 567890599L, 567890621L, 567890629L , 0}; + + Int4 o = 0; + Int4 pi = BigPrimeNumber[1]; + for (int i=0; BigPrimeNumber[i]; i++) { + Int4 r = BigPrimeNumber[i] % n; + Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); + if ( o < oo) + o=oo,pi=BigPrimeNumber[i];} + // cout << " AGoodNumberPrimeWith " << n << " " <Number(ta) << "::" + << CurrentTh->Number(ta.ns[0]) << "," + << CurrentTh->Number(ta.ns[1]) << "," + << CurrentTh->Number(ta.ns[2]) << "," + << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} " + << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} " + << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} " + << "]" ; + else + f << "[" + << ta.ns[0] << "," + << ta.ns[1] << "," + << ta.ns[2] << "," + << "{" << ta.at[0] << " " << ta.aa[0] << "} " + << "{" << ta.at[1] << " " << ta.aa[1] << "} " + << "{" << ta.at[2] << " " << ta.aa[2] << "} " + << "]" ; + return f;} + + + int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , + TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) + { // l'arete ta coupe l'arete pva pvb + // de cas apres le swap sa coupe toujours + // on cherche l'arete suivante + // on suppose que detsa >0 et detsb <0 + // attention la routine echange pva et pvb + + if(tt1.Locked()) return 0; // frontiere croise + + TriangleAdjacent tt2 = Adj(tt1); + Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent + Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles + assert ( a1 >= 0 && a1 < 3 ); + + Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; + Vertex & s1= (*t1)[OppositeVertex[a1]]; + Vertex & s2= (*t2)[OppositeVertex[a2]]; + + + Icoor2 dets2 = det(*pva,*pvb,s2); + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + assert((det1>0 ) && (det2 > 0)); + assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb + Icoor2 ndet1 = bamg::det(s1,sa,s2); + Icoor2 ndet2 = detT - ndet1; + + int ToSwap =0; //pas de swap + if ((ndet1 >0) && (ndet2 >0)) + { // on peut swaper + if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) + ToSwap =1; + else // swap alleatoire + if (BinaryRand()) + ToSwap =2; + } + if (ToSwap) NbSwap++, + bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); + + int ret=1; + + if (dets2 < 0) {// haut + dets1 = ToSwap ? dets1 : detsa ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = ToSwap ? dets1 : detsb ; + detsb = dets2; + //xxxx tt1 = ToSwap ? tt1 : Next(tt2); + if(!ToSwap) tt1 = Next(tt2); + } + else { // changement de sens + ret = -1; + Exchange(pva,pvb); + Exchange(detsa,detsb); + Exchange(dets1,dets2); + Exchange(tt1,tt2); + dets1=-dets1; + dets2=-dets2; + detsa=-detsa; + detsb=-detsb; + + if (ToSwap) + if (dets2 < 0) {// haut + dets1 = (ToSwap ? dets1 : detsa) ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = (ToSwap ? dets1 : detsb) ; + detsb = dets2; + if(!ToSwap) tt1 = Next(tt2); + } + else {// on a fin ??? + tt1 = Next(tt2); + ret =0;} + + } + return ret; + } + + int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) + { + int NbSwap =0; + assert(a.t && b.t); // the 2 vertex is in a mesh + int k=0; + taret=TriangleAdjacent(0,0); // erreur + + TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); + Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; + // we turn around a in the direct sens + + Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; + if(v2) // normal case + det2 = det(*v2,a,b); + else { // no chance infini vertex try the next + tta= Previous(Adj(tta)); + v2 = tta.EdgeVertex(0); + vbegin =v2; + assert(v2); + det2 = det(*v2,a,b); + // cout << " No Change try the next" << endl; + } + + while (v2 != &b) { + TriangleAdjacent tc = Previous(Adj(tta)); + v1 = v2; + v2 = tc.EdgeVertex(0); + det1 = det2; + det2 = v2 ? det(*v2,a,b): det2; + + if((det1 < 0) && (det2 >0)) { + // try to force the edge + Vertex * va = &a, *vb = &b; + tc = Previous(tc); + assert ( v1 && v2); + Icoor2 detss = 0,l=0,ks; + // cout << "Real ForcingEdge " << *va << *vb << detss << endl; + while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) + if(l++ > 10000000) { + cerr << " Loop in forcing Egde AB" + <<"\n vertex A " << a + <<"\n vertex B " << b + <<"\n nb de swap " << NbSwap + <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl; + + if ( CurrentTh ) + cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; + MeshError(990); + } + Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); + if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { + tc.SetLock(); + a.Optim(1,0); + b.Optim(1,0); + taret = tc; + return NbSwap; + } + else + { + taret = tc; + return -2; // error boundary is crossing + } + } + tta = tc; + assert(k++<2000); + if ( vbegin == v2 ) return -1;// error + } + + tta.SetLock(); + taret=tta; + a.Optim(1,0); + b.Optim(1,0); + return NbSwap; + } + }