![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2802
- Timestamp:
- 01/12/10 14:13:29 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 5 edited
-
Mesh2.cpp (modified) (37 diffs)
-
Mesh2.h (modified) (2 diffs)
-
MeshGeom.cpp (modified) (1 diff)
-
MeshQuad.cpp (modified) (4 diffs)
-
Metric.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.cpp
r2801 r2802 42 42 43 43 namespace bamg { 44 45 46 #ifdef DEBUG147 extern int SHOW ; // for debugging48 int SHOW = 0; // for debugging49 50 #endif51 44 52 45 int Triangles::counter = 0; … … 146 139 int ap1 = PreviousEdge[a1]; 147 140 int ap2 = PreviousEdge[a2]; 148 #ifdef DEBUG1149 t1->check();150 t2->check();151 #endif152 141 (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb 153 142 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa … … 157 146 tas1(t1,as1), tas2(t2,as2), 158 147 ta1(t1,a1),ta2(t2,a2); 159 #ifdef DEBUG160 assert( ! ta1.Locked());161 assert( ! ta2.Locked());162 #endif163 148 // externe haut gauche 164 149 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); … … 176 161 t1->SetTriangleContainingTheVertex(); 177 162 t2->SetTriangleContainingTheVertex(); 178 #ifdef DEBUG1179 t1->check();180 t2->check();181 #endif182 163 } // end swap 183 164 … … 263 244 { 264 245 // turn around the vertex ns[i] also call s 265 #ifdef DEBUG266 register Vertex * s = ns[i];267 #endif268 246 Triangle *t = (Triangle *) this , *ttc; 269 247 int k=0,j = EdgesVertexTriangle[i][0],jc; … … 274 252 int exterieurp = exterieur; 275 253 k++; 276 #ifdef DEBUG277 assert( s == & (*t)[VerticesOfTriangularEdge[j][1]] );278 #endif279 254 ttc = t->at[j]; 280 255 exterieur = !ttc->link; … … 673 648 dt[1]=bamg::det((*t)[2],(*t)[0],b); 674 649 dt[2]=bamg::det((*t)[0],(*t)[1],b); 675 #ifdef DEBUG676 assert(dt[0] >= 0);677 assert(dt[1] >= 0);678 assert(dt[2] >= 0);679 #endif680 650 Real8 dd = t->det; 681 651 NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd); … … 728 698 if ((v=(*tt)(2))) m2 = v->m; 729 699 lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2); 730 #ifdef DEBUG1731 if(SHOW) { cout << "SHOW ++ NewItem =" << Size << x ;732 cout << " " << d0 << " " << d1 << " " << d2 <<endl;}733 #endif734 700 n=Size++;} 735 701 else n=Size-1; … … 743 709 lIntTria[Size].x=A; 744 710 lIntTria[Size].m=mm; 745 #ifdef DEBUG1746 if (SHOW) cout << "SHOW ++ NewItem A" << Size << A << endl;747 #endif748 711 n=Size++; 749 712 } … … 836 799 sx =lIntTria[ii1].s; 837 800 Metric Mx=lIntTria[ii1].m; 838 #ifdef DEBUG839 double lx = lIntTria[ii-1].sn;840 #endif841 801 y =lIntTria[ii].x; 842 802 sy =lIntTria[ii].s; 843 803 Metric My=lIntTria[ii].m; 844 #ifdef DEBUG845 double ly =lIntTria[ii].sp;846 assert( sx <= si);847 assert( si <= sy);848 assert( sy != sx);849 #endif850 851 804 Real8 lxy = sy-sx; 852 805 Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy); … … 863 816 cout << " -- Add point " << nbv-1 << " " << vertices[nbv-1] << " " << vertices[nbv-1].m << endl; 864 817 865 #ifdef DEBUG866 if(k>1) {867 R2 AB = vertices[nbv-2].r - vertices[nbv-1].r ;868 Real8 dp = LengthInterpole(vertices[nbv-2].m,vertices[nbv-1].m,AB);869 if (dp > 1.6) {870 cerr << "PB calcul new Int. points trop loin l=" << dp << " v=" << nbv-1 << " " << nbv-2 <<Mx<<My<<y-x << endl;871 }872 }873 #endif874 818 } 875 819 else return nbv-nbvold; … … 899 843 900 844 Icoor2 dets2 = det(*pva,*pvb,s2); 901 902 #ifdef DEBUG903 Vertex & sb= (*t1)[VerticesOfTriangularEdge[a1][1]];904 Icoor2 wdets1 = det(*pva,*pvb,s1);905 Icoor2 wdetsa = det(*pva,*pvb,sa);906 Icoor2 wdetsb = det(*pva,*pvb,sb);907 assert(wdets1 == dets1);908 assert(wdetsa == detsa);909 assert(wdetsb == detsb);910 #endif911 912 845 Icoor2 det1=t1->det , det2=t2->det ; 913 #ifdef DEBUG914 assert(det1>0 && det2 >0);915 Icoor2 ddet1 = det((*t1)[0],(*t1)[1],(*t1)[2]);916 Icoor2 ddet2 = det((*t2)[0],(*t2)[1],(*t2)[2]);917 if ((det1 != ddet1) || (det2 != ddet2) )918 {919 assert(det1 == ddet1);920 assert(det2 == ddet2);921 }922 Icoor2 detvasasb = det(*pva,sa,sb);923 Icoor2 detvbsasb = det(*pvb,sa,sb);924 if ( CurrentTh && ! ( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0))))925 {926 cout << " detvasasb =" << detvasasb << "detvbsasb = " << detvbsasb927 << " " << pva << " " << pvb << " " <<CurrentTh <<endl;928 }929 assert( ( (detvasasb <= 0) && (detvbsasb >= 0)) || ( (detvasasb >= 0) && (detvbsasb <= 0)));930 #endif931 932 846 Icoor2 detT = det1+det2; 933 847 assert((det1>0 ) && (det2 > 0)); … … 945 859 ToSwap =2; 946 860 } 947 #ifdef DEBUG948 if (ForDebugging) {949 cerr << "swap = " << ToSwap << " ndet1 " << ndet1 << ", ndet2 " << ndet2 << "det1 " << det1 << " det2 " << det2950 << " if1 = " << ((ndet1 >0) && (ndet2 >0))951 << " if2 = " << ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) << endl;952 }953 #endif954 861 if (ToSwap) NbSwap++, 955 862 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 956 863 957 #ifdef DEBUG958 if (ForDebugging) {959 }960 #endif961 864 int ret=1; 962 865 … … 1003 906 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) 1004 907 { 1005 #ifdef DEBUG1006 restart: // for debug1007 #endif1008 908 int NbSwap =0; 1009 909 assert(a.t && b.t); // the 2 vertex is in a mesh … … 1032 932 v2 = tc.EdgeVertex(0); 1033 933 det1 = det2; 1034 #ifdef DEBUG1035 assert( v1 == tta.EdgeVertex(0));1036 assert( &a == tc.EdgeVertex(1) );1037 #endif1038 934 det2 = v2 ? det(*v2,a,b): det2; 1039 935 … … 1045 941 Icoor2 detss = 0,l=0,ks; 1046 942 // cout << "Real ForcingEdge " << *va << *vb << detss << endl; 1047 #ifdef DEBUG1048 Icoor2 dettt1 = det(*v1,a,b);1049 Icoor2 dettt2 = det(*v2,a,b);1050 1051 if (!(dettt1==det1 && dettt2==det2))1052 {1053 assert(ForDebugging==0);1054 ForDebugging=1;1055 goto restart;1056 }1057 1058 #endif1059 943 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 1060 944 if(l++ > 10000000) { … … 1067 951 if ( CurrentTh ) 1068 952 cerr << " vertex number " << CurrentTh->Number(a) << " " << CurrentTh->Number(b) << endl; 1069 #ifdef DEBUG1070 ForDebugging = 1;1071 #endif1072 953 MeshError(990); 1073 954 } … … 1108 989 1109 990 int Triangle::swap(Int2 a,int koption){ 1110 #ifdef DEBUG1111 if(a &4 ) return 0;// arete lock1112 int munswap1 = a/4;1113 a &=3;1114 #else1115 991 if(a/4 !=0) return 0;// arete lock or MarkUnSwap 1116 #endif1117 992 1118 993 register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent 1119 994 register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles 1120 #ifdef DEBUG1121 if(a2 & 4) return 0; // arete lock1122 int munswap2 = a2/4;1123 a2 &= 3;1124 #else1125 995 if(a2/4 !=0) return 0; // arete lock or MarkUnSwap 1126 #endif1127 996 1128 997 register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; … … 1130 999 register Vertex *s1=t1->ns[OppositeVertex[a1]]; 1131 1000 register Vertex *s2=t2->ns[OppositeVertex[a2]]; 1132 1133 #ifdef DEBUG1134 assert ( a >= 0 && a < 3 );1135 #endif1136 1001 1137 1002 Icoor2 det1=t1->det , det2=t2->det ; … … 1232 1097 } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) 1233 1098 } 1234 #ifdef DEBUG11235 if (OnSwap && ( munswap1 || munswap2)) {1236 cout << " erreur Mark unswap T " << CurrentTh->Number(t1) << " " << CurrentTh->Number(t2) << endl1237 << *t1 << endl1238 << *t2 << endl;1239 return 0;1240 }1241 #endif1242 1099 if( OnSwap ) 1243 1100 bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); … … 1251 1108 Real8 Vertex::Smoothing(Triangles & Th,const Triangles & BTh,Triangle * & tstart ,Real8 omega) 1252 1109 { 1253 #ifdef DEBUG1254 register Int4 NbSwap =0;1255 #endif1256 1110 register Vertex * s = this; 1257 1111 Vertex &vP = *s,vPsave=vP; … … 1267 1121 k++; 1268 1122 1269 #ifdef DEBUG1270 assert( s == & (*tria)[VerticesOfTriangularEdge[j][1]] );1271 assert( tria->det >0);1272 #endif1273 1123 if (!tria->Hidden(j)) 1274 1124 { … … 1423 1273 TriangleAdjacent ta = t->Adj(iedge); 1424 1274 1425 #ifdef DEBUG11426 cout << " the point " << Number(s) << " is the edge " << izerodet1427 << " of " << Number(t) << " det3 = "1428 << det3[0] << " " << det3[1] << " " << det3[2] << " " << endl;1429 cout << " ta = " << ta << "ta->det =" << ((Triangle*) ta)->det1430 << " "<< t->det<< endl;1431 #endif1432 1433 1275 // the point is on the edge 1434 1276 // if the point is one the boundary … … 1505 1347 assert(rswap); 1506 1348 } 1507 1508 #ifdef DEBUG1509 tt[0]->check();1510 tt[1]->check();1511 tt[2]->check();1512 #endif1513 1349 } 1514 1350 … … 1680 1516 1681 1517 NbTSwap += NbSwap ; 1682 #ifdef DEBUG1683 {1684 Int4 NbErr=0;1685 Int4 i;1686 for (i=0;i<nbt;i++)1687 if (triangles[i].link)1688 {1689 double dd =Det(triangles[i][1].r-triangles[i][0].r,triangles[i][2].r-triangles[i][0].r);1690 if(dd <=0)1691 {1692 NbErr++;1693 cerr << " det triangle i " << i << " = " << triangles[i].det ;1694 cerr << " det triangle " << dd ;1695 cerr << " Les trois sommets " ;1696 cerr << Number(triangles[i][0]) << " " << Number(triangles[i][1]) << " "1697 << Number(triangles[i][2]) << endl;1698 cerr << "I2 " <<triangles[i][0].r << triangles[i][1].r << triangles[i][2].r << endl;1699 cerr << "R2 " << triangles[i][0].i << triangles[i][1].i << triangles[i][2].i << endl;1700 cerr << "I2-R2 =" <<toR2(triangles[i][0].i)-triangles[i][0].r1701 << toR2(triangles[i][1].i)-triangles[i][1].r1702 << toR2(triangles[i][2].i)-triangles[i][2].r << endl;1703 }1704 }1705 if(NbErr) {1706 }1707 1708 }1709 #endif1710 1518 return nbv-nbvold; 1711 1519 … … 2068 1876 if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1; 2069 1877 } 2070 #ifdef DEBUG22071 cout << "ii1 = " << ii12072 << " ii0 = " << ii0 << endl;2073 cout << " cccc = " << cii0 << " " << ciii2074 << " " << cii1 << " sss=" << sss << endl;2075 #endif2076 1878 if (kstack >5) // bug ? 2077 1879 cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl; … … 2216 2018 } while (nbv!=nbvold); 2217 2019 delete [] first_np_or_next_t; 2218 #ifdef DEBUG2219 int nberr=0;2220 for (int it=0;it<nbt;it++)2221 if(triangles[it].det==0)2222 if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;2223 if(nberr) MeshError(992,this);2224 #endif2225 2020 cout << " end : Triangles::NewPoints old nbv=" << nbv << endl; 2226 2021 … … 2291 2086 triangles[1].link=&triangles[0]; 2292 2087 2293 #ifdef DEBUG2294 triangles[0].check();2295 triangles[1].check();2296 #endif2297 2088 // nbtf = 2; 2298 2089 if ( !quadtree ) quadtree = new QuadTree(this,0); … … 2352 2143 ReMakeTriangleContainingTheVertex(); 2353 2144 // because we break the TriangleContainingTheVertex 2354 #endif2355 #ifdef DEBUG2356 int nberr=0;2357 for (int it=0;it<nbt;it++)2358 if(triangles[it].det==0)2359 if(nberr++<10) cerr << "Bug Triangle nul" << it << triangles[it] << endl;2360 if(nberr) MeshError(991,this);2361 2145 #endif 2362 2146 time3=CPUtime(); … … 2847 2631 nt = nbt - NbOutT; 2848 2632 2849 #ifdef DEBUG2850 // verif2851 for ( it=0;it<nbt;it++)2852 triangles[it].check();2853 #endif2854 2633 } 2855 2634 Int4 Triangles::ConsRefTriangle(Int4 *reft) const … … 3257 3036 Real8 se= (sNew-L0)/LAB; 3258 3037 assert(se>=0 && se < 1.000000001); 3259 #ifdef DEBUG3260 se = abscisseInterpole(v0.m,v1.m,AB,se); // because code \ref(xxx)3261 #else3262 3038 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 3263 #endif3264 3039 assert(se>=0 && se <= 1); 3265 3040 //((k1==1) != (k1==k1equi)) … … 3279 3054 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl; 3280 3055 3281 #ifdef DEBUG3282 // code \label(xxx)3283 R2 A1A0 = A1->r - A0->r;3284 Real8 dp = LengthInterpole(A0->m,A1->m,A1A0);3285 if (dp > 1.4) {3286 cerr << " PB new Points "<< nbv-1 ;3287 cerr << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= " ;3288 cerr << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<endl;3289 cerr << " PB calcul new on cuver points trop loin l=" << dp3290 << " v=" << nbv-1 << " " << nbv-2 << " Lcurve = " << Lcurve << AB <<v0.m<< v1.m <<endl;3291 }3292 3293 #endif3294 3056 e->ref = eeequi.ref; 3295 3057 e->adj[0]=PreviousNewEdge; … … 3593 3355 PreviousNewEdge = edges + nbe; 3594 3356 nbe++; 3595 #ifdef DEBUG13596 cout << " new points " << nbv-1 << " " << vb->r ;3597 cout << " new edge " << nbe-1 << " " ;3598 cout << va << vb << " kk0 = " << kk03599 << " " << kk1 << " ss=" << ss ;3600 cout << " " << sbb << endl;3601 cout << " " << aa << va->r << bb << vb->r3602 <<" length=" << Norme(va->r-vb->r) << endl;3603 cout << " s " << s << " lstep= " << lstep3604 << " ledge= " << ledge3605 << " lcurve= " << lcurve << endl;3606 #endif3607 3357 va = vb; 3608 3358 } … … 4193 3943 // turne around in positif sens 4194 3944 Int4 NbSwap =0; 4195 #ifdef DEBUG4196 Vertex * s = ns[i];4197 #endif4198 3945 Triangle *t = this; 4199 3946 int k=0,j =OppositeEdge[i]; … … 4202 3949 Triangle *tp= at[jp]; 4203 3950 jp = aa[jp]&3; 4204 #ifdef DEBUG4205 assert(tp->at[jp] == this);4206 #endif4207 3951 do { 4208 #ifdef DEBUG4209 assert(k++<20000);4210 assert( s == & (*t)[OppositeVertex[j]] );4211 #endif4212 3952 // cout << *t << " " << j << "\n\t try swap " ; 4213 3953 while (t->swap(j,koption)) … … 4218 3958 j= NextEdge[tp->aa[jp]&3]; 4219 3959 // cout << "\n\t s " << *t << " " << j << endl; 4220 #ifdef DEBUG4221 assert( s == & (*t)[OppositeVertex[j]] );4222 #endif4223 3960 } 4224 3961 // end on this Triangle … … 4724 4461 dete[jn] = t->det-dete[j] -dete[jp]; 4725 4462 4726 #ifdef DEBUG4727 const Vertex * s0 = (*t)(0);4728 const Vertex * s1 = (*t)(1);4729 const Vertex * s2 = (*t)(2);4730 assert(dete[0] == det(B ,*s1,*s2));4731 assert(dete[1] == det(*s0,B ,*s2));4732 assert(dete[2] == det(*s0,*s1,B ));4733 assert(t->det== (dete[0] + dete[1] +dete[2]));4734 #endif4735 4463 // count the number k of dete <0 4736 4464 int k=0,ii[3]; -
issm/trunk/src/c/Bamgx/Mesh2.h
r2801 r2802 450 450 t->aa[aa[a] & 3] &=55; // 23 + 32 451 451 aa[a] &=55 ;} 452 453 454 455 #ifdef DEBUG456 void inline checka(Int1 a);457 void inline check();458 #endif459 452 460 453 }; // end of Triangle class … … 1309 1302 } 1310 1303 1311 #ifdef DEBUG 1312 void inline Triangle::checka(Int1 a) { 1313 // verif de la coherence des adjacences de l arete a 1314 a = a%4; 1315 assert(a < 3 && a >= 0 ); 1316 Triangle *t1=this,*t2=at[a]; 1317 Int2 a1=a,a2=aa[a]%4; 1318 1319 assert(a2 < 3 && a2 >= 0 ); 1320 if (t2 && ( ((*t1).ns[VerticesOfTriangularEdge[a1][0]] != (*t2).ns[VerticesOfTriangularEdge[a2][1]]) 1321 || ((*t1).ns[VerticesOfTriangularEdge[a1][1]] != (*t2).ns[VerticesOfTriangularEdge[a2][0]]))) 1322 { 1323 if (CurrentTh) std::cerr << " In Triangles beetween Triangle " << CurrentTh->Number(t1) << " and " 1324 << CurrentTh->Number(t2) << std::endl; 1325 std::cerr << "---- t1="<< t1 << " " << a1 <<", t2="<< t2 << " " << a2 << std::endl; 1326 std::cerr <<"t1="<< t1 << " " << a1 << " " << t1->ns[VerticesOfTriangularEdge[a1][0]] 1327 << " " << t1->ns[VerticesOfTriangularEdge[a1][1]] <<std::endl; 1328 if (CurrentTh) 1329 std::cerr <<"t1="<< t1 << " " << a1 << " " << CurrentTh->Number(t1->ns[VerticesOfTriangularEdge[a1][0]]) 1330 << " " << CurrentTh->Number(t1->ns[VerticesOfTriangularEdge[a1][1]]) <<std::endl; 1331 if (t2) std::cerr <<"t2="<< t2 << " " << a2 << " " 1332 << t2->ns[VerticesOfTriangularEdge[a2][0]] 1333 << " " << t2->ns[VerticesOfTriangularEdge[a2][1]] <<std::endl; 1334 if (t2 &&CurrentTh) 1335 std::cerr <<"t2="<< t2 << " " << a2 << " " 1336 << CurrentTh->Number(t2->ns[VerticesOfTriangularEdge[a2][0]]) 1337 << " " << CurrentTh->Number(t2->ns[VerticesOfTriangularEdge[a2][1]]) <<std::endl; 1338 assert(0); 1339 } 1340 if (t2) assert(t1->aa[a1]/4 == t2->aa[a2]/4); // lock compatibite 1341 } 1342 1343 1344 void inline Triangle::check() { 1345 Icoor2 det2=0; 1346 // std::cout << " check " << this << std::endl; 1347 int infv=ns[0] ? (( ns[1] ? ( ns[2] ? -1 : 2) : 1 )) : 0; 1348 if (det<0) { 1349 if (infv<0 ) 1350 { if (CurrentTh) std::cerr << " In Triangles " << CurrentTh->Number(this) << std::endl; 1351 std::cerr << " det = " << det << " and " << infv << std::endl; 1352 MeshError(5); 1353 }} 1354 else if (infv>=0 ) 1355 { if (CurrentTh) std::cerr << " In Triangles " << CurrentTh->Number(this) << std::endl; 1356 std::cerr << " det = " << det << " and " << infv << std::endl; 1357 MeshError(5); 1358 } 1359 1360 if (det >=0) 1361 if( det != (det2=bamg::det(*ns[0],*ns[1],*ns[2]))) 1362 { // penthickness(4);Draw(); 1363 if (CurrentTh) std::cerr << " In Triangles" << CurrentTh->Number(this) 1364 << std::endl; 1365 std::cerr << *ns[0] << *ns[1] << " " << *ns[2] << " " << std::endl; 1366 std::cerr << " Bug in triangle " << this 1367 << ":" << det << " != " << det2 << std::endl; 1368 MeshError(5); 1369 } 1370 checka(0); 1371 checka(1); 1372 checka(2); 1373 // if (ns[0]) assert( ns[0] - Meshbegin >= 0 ); 1374 // if (ns[0]) assert( Meshend - ns[0] >= 0 ); 1375 // if (ns[1]) assert( ns[1] - Meshbegin >= 0 ); 1376 // if (ns[1]) assert( Meshend - ns[1] >= 0 ); 1377 // if (ns[2]) assert( ns[2] - Meshbegin >= 0 ); 1378 // if (ns[2]) assert( Meshend - ns[2] >= 0 ); 1379 assert(ns[0] != ns[2]); 1380 assert(ns[1] != ns[2]); 1381 assert(ns[0] != ns[1]); 1382 } 1383 1384 1304 } 1385 1305 #endif 1386 1387 }1388 #endif -
issm/trunk/src/c/Bamgx/MeshGeom.cpp
r2740 r2802 577 577 if (OppositeSens) 578 578 s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); 579 // find all the discretisation of the egde580 #ifdef DEBUG581 // we supposee edge on=[Ag,Bg] intersect edge [V0,V1];582 // => V0Ag.V0V1 > 0 || V0Bg.V0V1 >0583 // => V1Ag.V0V1 < 0 || V0Bg.V0V1 <0584 R2 V0V1 = V1-V0;585 Real8 cos0A = ((Ag-V0),V0V1);586 Real8 cos0B = ((Bg-V0),V0V1);587 Real8 cos1A = ((Ag-V1),V0V1);588 Real8 cos1B = ((Bg-V1),V0V1);589 if ( (cos0A < 0 && cos0B <0) || (cos1A> 0 && cos1B >0))590 {591 cerr << " Bug on pointer edge [" << V0 << " , " << V1 << " ] "592 << " on geometrical edge " << Number(on) << " = [" << Ag << " , " << Bg << " ] " << endl;593 cerr << cos0A << "> 0 || " << cos0B << " > 0 and ";594 cerr << cos1A << "< 0 || " << cos1B << " < 0 " << endl;595 596 exit (1);597 }598 599 #endif600 579 if(NbTry) cout << "bug: edge = " << v0.r << " -> " << v1.r << endl 601 580 << "sg 0 = " << vg0 -
issm/trunk/src/c/Bamgx/MeshQuad.cpp
r2800 r2802 178 178 // 1 first find a back ground edge contening the vertex A 179 179 // 2 walk n back gound boundary to find the final vertex B 180 181 182 #ifdef DEBUG183 // we suppose if the vertex A is on a background edge and184 // the abcisse is 0 or 1 when the related point is not185 // a end of curve <=> !IsRequiredVertex186 if (vA.vint == IsVertexOnEdge)187 if (tA<=0)188 assert(! (*vA.onbe->be)[0].on->IsRequiredVertex());189 else if (tA>=1)190 assert(!(*vA.onbe->be)[1].on->IsRequiredVertex());191 #endif192 180 193 181 if( vA.vint == IsVertexOnEdge) … … 419 407 // memcpy(newedges,edges,sizeof(Edge)*nbe); 420 408 SetOfEdges4 * edge4= new SetOfEdges4(nbe,nbv); 421 #ifdef DEBUG422 for (i=0;i<nbt;i++)423 triangles[i].check();424 #endif425 409 Int4 k=nbv; 426 410 Int4 kk=0; … … 512 496 k++; 513 497 } 514 #ifdef DEBUG515 assert(kvb == newNbVertexOnBThEdge);516 // verif edge517 { Vertex *v0 = vertices, *v1 = vertices+ k;518 for (Int4 i=0;i<ie;i++)519 {520 assert( &newedges[i][0] >= v0 && &newedges[i][0] < v1);521 assert( &newedges[i][1] >= v0 && &newedges[i][1] < v1);522 }523 }524 #endif525 498 if (edgesGtoB) delete [] edgesGtoB; 526 499 edgesGtoB=0; … … 897 870 FillHoleInMesh(); 898 871 899 #ifdef DEBUG900 for (i=0;i<nbt;i++)901 triangles[i].check();902 #endif903 904 872 if (verbosity>2) 905 873 cout << " (out) Nb of Quadrilaterals = " << NbOfQuad -
issm/trunk/src/c/Bamgx/Metric.cpp
r2800 r2802 97 97 else 98 98 s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1; 99 #ifdef DEBUG 100 assert(Abs(vp1y)+Abs(vp2y)>0); 101 #endif 102 V=D2xD2(vp1x,vp2x,vp1y,vp2y); 99 V=D2xD2(vp1x,vp2x,vp1y,vp2y); 103 100 } 104 101 return;
Note:
See TracChangeset
for help on using the changeset viewer.
