![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2848
- Timestamp:
- 01/14/10 17:42:33 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 2 edited
-
Mesh2.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.cpp
r2847 r2848 18 18 int ForDebugging = 0; 19 19 const Direction NoDirOfSearch = Direction(); 20 21 Int4 AGoodNumberPrimeWith(Int4 n){22 const Int4 BigPrimeNumber[] ={ 567890359L,23 567890431L, 567890437L, 567890461L, 567890471L,24 567890483L, 567890489L, 567890497L, 567890507L,25 567890591L, 567890599L, 567890621L, 567890629L , 0};26 27 Int4 o = 0;28 Int4 pi = BigPrimeNumber[1];29 for (int i=0; BigPrimeNumber[i]; i++) {30 Int4 r = BigPrimeNumber[i] % n;31 Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r)));32 if ( o < oo)33 o=oo,pi=BigPrimeNumber[i];}34 // cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl;35 return pi;36 }37 38 void MeshError(int Err,Triangles *Th){39 cerr << " Fatal error in the meshgenerator " << Err << endl ;40 exit(1);41 }42 43 ostream& operator <<(ostream& f, const Triangle & ta) {44 if(CurrentTh)45 f << "[" << CurrentTh->Number(ta) << "::"46 << CurrentTh->Number(ta.ns[0]) << ","47 << CurrentTh->Number(ta.ns[1]) << ","48 << CurrentTh->Number(ta.ns[2]) << ","49 << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} "50 << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} "51 << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} "52 << "]" ;53 else54 f << "["55 << ta.ns[0] << ","56 << ta.ns[1] << ","57 << ta.ns[2] << ","58 << "{" << ta.at[0] << " " << ta.aa[0] << "} "59 << "{" << ta.at[1] << " " << ta.aa[1] << "} "60 << "{" << ta.at[2] << " " << ta.aa[2] << "} "61 << "]" ;62 return f;}63 64 65 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,66 TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap)67 { // l'arete ta coupe l'arete pva pvb68 // de cas apres le swap sa coupe toujours69 // on cherche l'arete suivante70 // on suppose que detsa >0 et detsb <071 // attention la routine echange pva et pvb72 73 if(tt1.Locked()) return 0; // frontiere croise74 75 TriangleAdjacent tt2 = Adj(tt1);76 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent77 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles78 assert ( a1 >= 0 && a1 < 3 );79 80 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];81 Vertex & s1= (*t1)[OppositeVertex[a1]];82 Vertex & s2= (*t2)[OppositeVertex[a2]];83 84 85 Icoor2 dets2 = det(*pva,*pvb,s2);86 Icoor2 det1=t1->det , det2=t2->det ;87 Icoor2 detT = det1+det2;88 assert((det1>0 ) && (det2 > 0));89 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb90 Icoor2 ndet1 = bamg::det(s1,sa,s2);91 Icoor2 ndet2 = detT - ndet1;92 93 int ToSwap =0; //pas de swap94 if ((ndet1 >0) && (ndet2 >0))95 { // on peut swaper96 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))97 ToSwap =1;98 else // swap alleatoire99 if (BinaryRand())100 ToSwap =2;101 }102 if (ToSwap) NbSwap++,103 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);104 105 int ret=1;106 107 if (dets2 < 0) {// haut108 dets1 = ToSwap ? dets1 : detsa ;109 detsa = dets2;110 tt1 = Previous(tt2) ;}111 else if (dets2 > 0){// bas112 dets1 = ToSwap ? dets1 : detsb ;113 detsb = dets2;114 //xxxx tt1 = ToSwap ? tt1 : Next(tt2);115 if(!ToSwap) tt1 = Next(tt2);116 }117 else { // changement de sens118 if (ForDebugging) cout << "changement de sens" << endl;119 ret = -1;120 Exchange(pva,pvb);121 Exchange(detsa,detsb);122 Exchange(dets1,dets2);123 Exchange(tt1,tt2);124 dets1=-dets1;125 dets2=-dets2;126 detsa=-detsa;127 detsb=-detsb;128 129 if (ToSwap)130 if (dets2 < 0) {// haut131 dets1 = (ToSwap ? dets1 : detsa) ;132 detsa = dets2;133 tt1 = Previous(tt2) ;}134 else if (dets2 > 0){// bas135 dets1 = (ToSwap ? dets1 : detsb) ;136 detsb = dets2;137 if(!ToSwap) tt1 = Next(tt2);138 }139 else {// on a fin ???140 tt1 = Next(tt2);141 ret =0;}142 143 }144 return ret;145 }146 147 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret)148 {149 int NbSwap =0;150 assert(a.t && b.t); // the 2 vertex is in a mesh151 int k=0;152 taret=TriangleAdjacent(0,0); // erreur153 154 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]);155 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;156 // we turn around a in the direct sens157 158 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1;159 if(v2) // normal case160 det2 = det(*v2,a,b);161 else { // no chance infini vertex try the next162 tta= Previous(Adj(tta));163 v2 = tta.EdgeVertex(0);164 vbegin =v2;165 assert(v2);166 det2 = det(*v2,a,b);167 // cout << " No Change try the next" << endl;168 }169 170 while (v2 != &b) {171 TriangleAdjacent tc = Previous(Adj(tta));172 v1 = v2;173 v2 = tc.EdgeVertex(0);174 det1 = det2;175 det2 = v2 ? det(*v2,a,b): det2;176 177 if((det1 < 0) && (det2 >0)) {178 // try to force the edge179 Vertex * va = &a, *vb = &b;180 tc = Previous(tc);181 assert ( v1 && v2);182 Icoor2 detss = 0,l=0,ks;183 // cout << "Real ForcingEdge " << *va << *vb << detss << endl;184 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))185 if(l++ > 10000000) {186 cerr << " Loop in forcing Egde AB"187 <<"\n vertex A " << a188 <<"\n vertex B " << b189 <<"\n nb de swap " << NbSwap190 <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl;191 192 if ( CurrentTh )193 cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl;194 MeshError(990);195 }196 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);197 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) {198 tc.SetLock();199 a.Optim(1,0);200 b.Optim(1,0);201 taret = tc;202 return NbSwap;203 }204 else205 {206 taret = tc;207 return -2; // error boundary is crossing208 }209 }210 tta = tc;211 assert(k++<2000);212 if ( vbegin == v2 ) return -1;// error213 }214 215 tta.SetLock();216 taret=tta;217 a.Optim(1,0);218 b.Optim(1,0);219 return NbSwap;220 }221 222 20 } 223 -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2843 r2848 6024 6024 /*}}}1*/ 6025 6025 6026 Int4 AGoodNumberPrimeWith(Int4 n){ 6027 const Int4 BigPrimeNumber[] ={ 567890359L, 6028 567890431L, 567890437L, 567890461L, 567890471L, 6029 567890483L, 567890489L, 567890497L, 567890507L, 6030 567890591L, 567890599L, 567890621L, 567890629L , 0}; 6031 6032 Int4 o = 0; 6033 Int4 pi = BigPrimeNumber[1]; 6034 for (int i=0; BigPrimeNumber[i]; i++) { 6035 Int4 r = BigPrimeNumber[i] % n; 6036 Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); 6037 if ( o < oo) 6038 o=oo,pi=BigPrimeNumber[i];} 6039 // cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl; 6040 return pi; 6041 } 6042 6043 void MeshError(int Err,Triangles *Th){ 6044 cerr << " Fatal error in the meshgenerator " << Err << endl ; 6045 exit(1); 6046 } 6047 6048 ostream& operator <<(ostream& f, const Triangle & ta) { 6049 if(CurrentTh) 6050 f << "[" << CurrentTh->Number(ta) << "::" 6051 << CurrentTh->Number(ta.ns[0]) << "," 6052 << CurrentTh->Number(ta.ns[1]) << "," 6053 << CurrentTh->Number(ta.ns[2]) << "," 6054 << "{" << CurrentTh->Number(ta.at[0]) << " " << ta.aa[0] << "} " 6055 << "{" << CurrentTh->Number(ta.at[1]) << " " << ta.aa[1] << "} " 6056 << "{" << CurrentTh->Number(ta.at[2]) << " " << ta.aa[2] << "} " 6057 << "]" ; 6058 else 6059 f << "[" 6060 << ta.ns[0] << "," 6061 << ta.ns[1] << "," 6062 << ta.ns[2] << "," 6063 << "{" << ta.at[0] << " " << ta.aa[0] << "} " 6064 << "{" << ta.at[1] << " " << ta.aa[1] << "} " 6065 << "{" << ta.at[2] << " " << ta.aa[2] << "} " 6066 << "]" ; 6067 return f;} 6068 6069 6070 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , 6071 TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) 6072 { // l'arete ta coupe l'arete pva pvb 6073 // de cas apres le swap sa coupe toujours 6074 // on cherche l'arete suivante 6075 // on suppose que detsa >0 et detsb <0 6076 // attention la routine echange pva et pvb 6077 6078 if(tt1.Locked()) return 0; // frontiere croise 6079 6080 TriangleAdjacent tt2 = Adj(tt1); 6081 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 6082 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 6083 assert ( a1 >= 0 && a1 < 3 ); 6084 6085 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 6086 Vertex & s1= (*t1)[OppositeVertex[a1]]; 6087 Vertex & s2= (*t2)[OppositeVertex[a2]]; 6088 6089 6090 Icoor2 dets2 = det(*pva,*pvb,s2); 6091 Icoor2 det1=t1->det , det2=t2->det ; 6092 Icoor2 detT = det1+det2; 6093 assert((det1>0 ) && (det2 > 0)); 6094 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb 6095 Icoor2 ndet1 = bamg::det(s1,sa,s2); 6096 Icoor2 ndet2 = detT - ndet1; 6097 6098 int ToSwap =0; //pas de swap 6099 if ((ndet1 >0) && (ndet2 >0)) 6100 { // on peut swaper 6101 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 6102 ToSwap =1; 6103 else // swap alleatoire 6104 if (BinaryRand()) 6105 ToSwap =2; 6106 } 6107 if (ToSwap) NbSwap++, 6108 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 6109 6110 int ret=1; 6111 6112 if (dets2 < 0) {// haut 6113 dets1 = ToSwap ? dets1 : detsa ; 6114 detsa = dets2; 6115 tt1 = Previous(tt2) ;} 6116 else if (dets2 > 0){// bas 6117 dets1 = ToSwap ? dets1 : detsb ; 6118 detsb = dets2; 6119 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 6120 if(!ToSwap) tt1 = Next(tt2); 6121 } 6122 else { // changement de sens 6123 ret = -1; 6124 Exchange(pva,pvb); 6125 Exchange(detsa,detsb); 6126 Exchange(dets1,dets2); 6127 Exchange(tt1,tt2); 6128 dets1=-dets1; 6129 dets2=-dets2; 6130 detsa=-detsa; 6131 detsb=-detsb; 6132 6133 if (ToSwap) 6134 if (dets2 < 0) {// haut 6135 dets1 = (ToSwap ? dets1 : detsa) ; 6136 detsa = dets2; 6137 tt1 = Previous(tt2) ;} 6138 else if (dets2 > 0){// bas 6139 dets1 = (ToSwap ? dets1 : detsb) ; 6140 detsb = dets2; 6141 if(!ToSwap) tt1 = Next(tt2); 6142 } 6143 else {// on a fin ??? 6144 tt1 = Next(tt2); 6145 ret =0;} 6146 6147 } 6148 return ret; 6149 } 6150 6151 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) 6152 { 6153 int NbSwap =0; 6154 assert(a.t && b.t); // the 2 vertex is in a mesh 6155 int k=0; 6156 taret=TriangleAdjacent(0,0); // erreur 6157 6158 TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); 6159 Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 6160 // we turn around a in the direct sens 6161 6162 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 6163 if(v2) // normal case 6164 det2 = det(*v2,a,b); 6165 else { // no chance infini vertex try the next 6166 tta= Previous(Adj(tta)); 6167 v2 = tta.EdgeVertex(0); 6168 vbegin =v2; 6169 assert(v2); 6170 det2 = det(*v2,a,b); 6171 // cout << " No Change try the next" << endl; 6172 } 6173 6174 while (v2 != &b) { 6175 TriangleAdjacent tc = Previous(Adj(tta)); 6176 v1 = v2; 6177 v2 = tc.EdgeVertex(0); 6178 det1 = det2; 6179 det2 = v2 ? det(*v2,a,b): det2; 6180 6181 if((det1 < 0) && (det2 >0)) { 6182 // try to force the edge 6183 Vertex * va = &a, *vb = &b; 6184 tc = Previous(tc); 6185 assert ( v1 && v2); 6186 Icoor2 detss = 0,l=0,ks; 6187 // cout << "Real ForcingEdge " << *va << *vb << detss << endl; 6188 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 6189 if(l++ > 10000000) { 6190 cerr << " Loop in forcing Egde AB" 6191 <<"\n vertex A " << a 6192 <<"\n vertex B " << b 6193 <<"\n nb de swap " << NbSwap 6194 <<"\n nb of try swap too big = " << l << " gearter than " << 1000000 << endl; 6195 6196 if ( CurrentTh ) 6197 cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; 6198 MeshError(990); 6199 } 6200 Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 6201 if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { 6202 tc.SetLock(); 6203 a.Optim(1,0); 6204 b.Optim(1,0); 6205 taret = tc; 6206 return NbSwap; 6207 } 6208 else 6209 { 6210 taret = tc; 6211 return -2; // error boundary is crossing 6212 } 6213 } 6214 tta = tc; 6215 assert(k++<2000); 6216 if ( vbegin == v2 ) return -1;// error 6217 } 6218 6219 tta.SetLock(); 6220 taret=tta; 6221 a.Optim(1,0); 6222 b.Optim(1,0); 6223 return NbSwap; 6224 } 6225 6026 6226 }
Note:
See TracChangeset
for help on using the changeset viewer.
