![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2854
- Timestamp:
- 01/15/10 11:23:20 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 9 edited
-
QuadTree.h (modified) (2 diffs)
-
objects/GeometricalEdge.cpp (modified) (3 diffs)
-
objects/Geometry.cpp (modified) (13 diffs)
-
objects/ListofIntersectionTriangles.cpp (modified) (3 diffs)
-
objects/MetricAnIso.cpp (modified) (3 diffs)
-
objects/QuadTree.cpp (modified) (4 diffs)
-
objects/Triangle.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (73 diffs)
-
objects/Vertex.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/QuadTree.h
r2841 r2854 1 #include "../shared/shared.h" 2 #include "../include/macros.h" 3 #include "../toolkits/toolkits.h" 4 1 5 namespace bamg { 2 6 … … 58 62 void Add( Vertex & w); 59 63 60 QuadTreeBox* NewQuadTreeBox() 61 { 62 ///cout << "NewQuadTreeBox " << sb << " " << sb->bc << " " 63 //<< sb->be << " " <<lenStorageQuadTreeBox <<endl; 64 QuadTreeBox* NewQuadTreeBox(){ 64 65 if(! (sb->bc<sb->be)) 65 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb);66 sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); 66 67 67 assert(sb && (sb->bc->n == 0)); 68 if (!sb || (sb->bc->n != 0)){ 69 throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)")); 70 } 68 71 NbQuadTreeBox++; 69 72 return sb->bc++; -
issm/trunk/src/c/Bamgx/objects/GeometricalEdge.cpp
r2850 r2854 8 8 #include "../QuadTree.h" 9 9 #include "../SetOfE4.h" 10 11 #include "../../shared/shared.h" 12 #include "../../include/macros.h" 13 #include "../../toolkits/toolkits.h" 10 14 11 15 #undef __FUNCT__ … … 27 31 // Real8 t1t1 = t1*t1; 28 32 Real8 tt = theta*theta; 29 assert( theta >=0); 30 assert( theta <=1); 33 if ( theta<0){ 34 throw ErrorException(__FUNCT__,exprintf("theta<0")); 35 } 36 if ( theta>1){ 37 throw ErrorException(__FUNCT__,exprintf("theta>1")); 38 } 31 39 if (TgA()) 32 40 if (TgB()) // interpolation d'hermite … … 92 100 { R2 A=v[0]->r,B=v[1]->r; 93 101 Real8 ca,cb,cta,ctb; 94 assert( theta >=-1e-12); 95 assert( theta <=1+1e-12); 102 if ( theta<-1e-12){ 103 throw ErrorException(__FUNCT__,exprintf("theta<-1e-12")); 104 } 105 if ( theta>1+1e-12){ 106 throw ErrorException(__FUNCT__,exprintf("theta>1+1e-12")); 107 } 96 108 if (TgA()) 97 109 if (TgB()) // interpolation d'hermite -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2850 r2854 42 42 for (i=0;i<NbSubDomains;i++) 43 43 subdomains[i].Set(Gh.subdomains[i],Gh,*this); 44 45 // for (i=0;i<nbt;i++) 46 // triangles[i].Set(Gh.triangles[i],Gh,*this); 47 assert(!nbt); 44 if (nbt); { 45 throw ErrorException(__FUNCT__,exprintf("nbt")); 46 } 48 47 } 49 48 /*}}}1*/ … … 52 51 long int verbosity=0; 53 52 54 assert(NbRef<=0); 53 if (NbRef>0){ 54 throw ErrorException(__FUNCT__,exprintf("NbRef>0")); 55 } 55 56 if(verbosity>9) 56 57 cout << "DELETE ~Geometry "<< this << endl; … … 512 513 // cout << P ; 513 514 int k=0; 514 while(pon != on) 515 {515 while(pon != on){ 516 k++; 516 517 pon = on; 517 assert(k++<100); 518 if (k>=100){ 519 throw ErrorException(__FUNCT__,exprintf("k>=100")); 520 } 518 521 R2 A= (*on)[0]; 519 522 R2 B= (*on)[1]; … … 538 541 retry: 539 542 s=save_s; 540 GeometricalEdge * on = e.on; 541 assert(on); 542 assert( e[0].on && e[1].on); 543 GeometricalEdge* on = e.on; 544 if (!on){ 545 throw ErrorException(__FUNCT__,exprintf("!on")); 546 } 547 if (!e[0].on || !e[1].on){ 548 throw ErrorException(__FUNCT__,exprintf("!e[0].on || !e[1].on")); 549 } 543 550 const Vertex &v0=e[0],&v1=e[1]; 544 551 V.m = Metric(1.0-s, v0,s, v1); … … 586 593 Number(eg0->Adj[1]) <<"," ; 587 594 ge[--bge] =eg0 = eg0->Adj[sens0]; 588 assert(bge>=0 && bge <= mxe); 595 if (bge<0 || bge>mxe){ 596 throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe")); 597 } 589 598 sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); 590 599 if(NbTry) … … 615 624 ge[++tge] =eg1 = eg1->Adj[sens1]; 616 625 sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; 617 assert(tge>=0 && tge <= mxe); 626 if (tge<0 || tge>mxe){ 627 throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)")); 628 } 618 629 if(NbTry) 619 630 cout << " Edge " << Number(eg1) << " " << sens1 << " S " … … 648 659 int i=bge; 649 660 Real8 ll=0; 650 for(i=bge;i<tge;i++) 651 { 652 assert( i>=0 && i <= mxe); 661 for(i=bge;i<tge;i++){ 662 if ( i<0 || i>mxe){ 663 throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); 664 } 653 665 BB = (*ge[i])[sensge[i]]; 654 666 lge[i]=ll += Norme2(AA-BB); … … 660 672 // <<" on = " << Number(ge[tge]) << " sens= " << sensge[tge] << endl; 661 673 // search the geometrical edge 662 assert(s <= 1.0); 674 if (s>1.0){ 675 throw ErrorException(__FUNCT__,exprintf("s>1.0")); 676 } 663 677 Real8 ls= s*ll; 664 678 on =0; … … 668 682 i=bge; 669 683 while ( (l1=lge[i]) < ls ) { 670 assert(i >= 0 && i <= mxe); 684 if (i<0 || i>mxe){ 685 throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); 686 } 671 687 i++,s0=1-(s1=sensge[i]),l0=l1;} 672 688 on=ge[i]; … … 679 695 sg = s0 * (1.0-s) + s * s1; 680 696 } 681 assert(on); 682 // assert(sg && sg-1); 697 if (!on){ 698 throw ErrorException(__FUNCT__,exprintf("!on")); 699 } 683 700 V.r= on->F(sg); 684 701 // if (eg0 != eg1) … … 716 733 GeometricalVertex * vg = (GeometricalVertex *) (void *) v; 717 734 int j = vg-v0g; 718 assert( v == & (Vertex &) vertices[j]); 735 if ( v != &(Vertex &) vertices[j]){ 736 throw ErrorException(__FUNCT__,exprintf(" v != &(Vertex &) vertices[j]")); 737 } 719 738 vertices[i].link = vertices + j; 720 739 k++; … … 949 968 } 950 969 951 }} 952 assert(nbgem && nbe); 970 } 971 } 972 if (nbgem==0 || nbe==0){ 973 throw ErrorException(__FUNCT__,exprintf("nbgem==0 || nbe==0")); 974 } 953 975 954 976 if(step==0) { … … 961 983 //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; 962 984 curves[i].master=true; 963 if(be->Equi() || be->ReverseEqui() ) 964 { 965 assert(eqbe); 985 if(be->Equi() || be->ReverseEqui() ){ 986 if (!eqbe){ 987 throw ErrorException(__FUNCT__,exprintf("!eqbe")); 988 } 966 989 int nc = eqbe->CurveNumber; 967 assert(i!=nc); 990 if (i==nc){ 991 throw ErrorException(__FUNCT__,exprintf("i==nc")); 992 } 968 993 curves[i].next=curves[nc].next; 969 994 curves[i].master=false; -
issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp
r2850 r2854 67 67 while (t->det <0) { // intersection boundary edge and a,b, 68 68 k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 69 assert(k>=0); 69 if (k<0){ 70 throw ErrorException(__FUNCT__,exprintf("k<0")); 71 } 70 72 ocut = OppositeEdge[k]; 71 73 i=VerticesOfTriangularEdge[ocut][0]; … … 114 116 j=VerticesOfTriangularEdge[iedge][1]; 115 117 Real8 dij = detj-deti; 116 assert(i+j+k == 0 + 1 +2); 118 if (i+j+k != 0 + 1 +2){ 119 throw ErrorException(__FUNCT__,exprintf("i+j+k != 0 + 1 +2")); 120 } 117 121 ba[j] = detj/dij; 118 122 ba[i] = -deti/dij; … … 268 272 Real8 ListofIntersectionTriangles::Length(){ 269 273 // cout << " n= " << Size << ":" ; 270 assert(Size>0); 274 if (Size<=0){ 275 throw ErrorException(__FUNCT__,exprintf("Size<=0")); 276 } 271 277 // computation of the length 272 278 R2 C; -
issm/trunk/src/c/Bamgx/objects/MetricAnIso.cpp
r2850 r2854 8 8 #include "../QuadTree.h" 9 9 #include "../SetOfE4.h" 10 11 #include "../../shared/shared.h" 12 #include "../../include/macros.h" 13 #include "../../toolkits/toolkits.h" 10 14 11 15 #undef __FUNCT__ … … 206 210 } 207 211 // warning for optimisation S is in [0:0.5] not in [0:1] 208 assert(i<512); 212 if (i>=512){ 213 throw ErrorException(__FUNCT__,exprintf("i>=512")); 214 } 209 215 LastMetricInterpole.lab=l; 210 216 LastMetricInterpole.opt=i; … … 243 249 r = 2*(S[i]*(L[j]-l)+ S[j]*(l-L[i]))/(L[j]-L[i]); 244 250 } 245 assert(r<=1 && r>=0); 251 if (r>1 || r<0){ 252 throw ErrorException(__FUNCT__,exprintf("r>1 || r<0")); 253 } 246 254 return r ; 247 255 } -
issm/trunk/src/c/Bamgx/objects/QuadTree.cpp
r2850 r2854 5 5 #include "../Mesh2.h" 6 6 #include "../QuadTree.h" 7 8 #include "../../shared/shared.h" 9 #include "../../include/macros.h" 10 #include "../../toolkits/toolkits.h" 7 11 8 12 #undef __FUNCT__ … … 33 37 sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); 34 38 root=NewQuadTreeBox(); 35 assert( MaxISize > MaxICoor); 39 if ( MaxISize <= MaxICoor){ 40 throw ErrorException(__FUNCT__,exprintf("MaxISize <= MaxICoor")); 41 } 36 42 for (Int4 i=0;i<nbv;i++) 37 43 Add(t->vertices[i]); … … 373 379 if (b->n > 0 && b->v[0] == &w) return; 374 380 } 375 assert(l); 381 if (l==0){ 382 throw ErrorException(__FUNCT__,exprintf("l==0")); 383 } 376 384 while ((b= *pb) && (b->n == 4)) // the QuadTreeBox is full 377 385 { … … 412 420 bc =b; 413 421 be = b +ll; 414 assert(b); 422 if (!b){ 423 throw ErrorException(__FUNCT__,exprintf("!b")); 424 } 415 425 } 416 426 /*}}}1*/ -
issm/trunk/src/c/Bamgx/objects/Triangle.cpp
r2850 r2854 37 37 t = ttc; 38 38 j = NextEdge[jc]; 39 assert(k<2000); 39 if (k>=2000){ 40 throw ErrorException(__FUNCT__,exprintf("k>=2000")); 41 } 40 42 } while ( (this!= t)); 41 43 return TriangleAdjacent(0,0); -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2853 r2854 2454 2454 Int4 BeginNewPoint[3]; 2455 2455 Int4 EndNewPoint[3]; 2456 #ifdef TRACETRIANGLE2457 Int4 trace=0;2458 #endif2459 2456 int step[3]; 2460 2457 Int4 *first_np_or_next_t = new Int4[nbtx]; … … 2480 2477 i=Headt; 2481 2478 next_t=-first_np_or_next_t[i]; 2482 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) 2483 { // for each triangle t 2479 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) { // for each triangle t 2484 2480 // we can change first_np_or_next_t[i] 2485 #ifdef TRACETRIANGLE 2486 trace = TRACETRIANGLE <0 ? 1 : i == TRACETRIANGLE; 2487 #endif 2488 assert(i>=0 && i < nbt ); 2481 if (i<0 || i >= nbt ){ 2482 throw ErrorException(__FUNCT__,exprintf("i<0 || i >= nbt")); 2483 } 2489 2484 first_np_or_next_t[i] = nbv; // to save the fist new point of triangle 2490 2485 for(j=0;j<3;j++) … … 2500 2495 Vertex & vB = * tj.EdgeVertex(1); 2501 2496 2502 #ifdef TRACETRIANGLE2503 if(trace) {2504 cout << " i " << Number(vA) <<" j "<< Number(vB)2505 << " " << t->Locked(j) ;2506 }2507 #endif2508 2497 if (!t->link) continue;// boundary 2509 2498 if (t->det <0) continue; … … 2523 2512 const Vertex & vC2 = *tadjj.OppositeVertex(); 2524 2513 2525 #ifdef TRACETRIANGLE2526 trace = trace || k == TRACETRIANGLE;2527 if(trace) {2528 cout << "Test Arete " << i << " AB = " << A << B2529 << "i " <<Number(vA)<< "j" <<Number(vB);2530 cout << " link" <<(int)t->link << " ta=" << Number( ta)2531 << " det " <<ta->det ;2532 cout << " hA = " <<vA.m.h << " hB = " <<vB.m.h ;2533 cout << " loc " << ta->Locked(j) << endl;2534 }2535 #endif2536 2537 2514 if(first_np_or_next_t[k]>0) { // this edge is done before 2538 2515 // find the color of the edge and begin , end of newpoint 2539 2516 register int kk = t->NuEdgeTriangleAdj(j); 2540 assert ((*t)(VerticesOfTriangularEdge[j][0]) == (*ta)(VerticesOfTriangularEdge[kk][1])); 2541 assert ((*t)(VerticesOfTriangularEdge[j][1]) == (*ta)(VerticesOfTriangularEdge[kk][0])); 2517 if ((*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])){ 2518 throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])")); 2519 } 2520 if ((*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])){ 2521 throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])")); 2522 } 2542 2523 register Int4 kolor =3*k + kk; 2543 2524 ColorEdge[j]=kolor; … … 2578 2559 EndNewPoint[j] = nbv-1; 2579 2560 } // end loop for each edge 2580 2581 #ifdef TRACETRIANGLE2582 if(trace) {2583 // verification des point cree2584 cout << "\n ------------ " << t->link << " " << t->det2585 << " b " << BeginNewPoint[0] << " " << BeginNewPoint[1]2586 << " " << BeginNewPoint[2] << " "2587 << " e " << EndNewPoint[0] << " " << EndNewPoint[1]2588 << " " << EndNewPoint[2] << " "2589 << " s " << step[0] << " " << step[1] << " " << step[2] << " "2590 << endl;2591 }2592 #endif2593 2561 2594 2562 if (!t->link) continue;// boundary … … 2664 2632 Int4 i1; 2665 2633 kstack =0; 2666 while( (i1=stack[kstack++]) >= 0) 2667 { // the two parameter is i0 and i1 2668 assert(i1 < nbv && i1 >= 0); 2669 assert(i0 < nbv && i0 >= 0); 2670 assert(i1 != i0); 2634 while( (i1=stack[kstack++]) >= 0){ // the two parameter is i0 and i1 2635 if (i1 >= nbv || i1 < 0){ 2636 throw ErrorException(__FUNCT__,exprintf("i1 >= nbv || i1 < 0")); 2637 } 2638 if (i0 >= nbv || i0 < 0){ 2639 throw ErrorException(__FUNCT__,exprintf("i0 >= nbv || i0 < 0")); 2640 } 2641 if (i1 == i0){ 2642 throw ErrorException(__FUNCT__,exprintf("i1 == i0")); 2643 } 2671 2644 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0]; 2672 2645 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01)); 2673 2674 2675 #ifdef TRACETRIANGLE 2676 if(trace) { 2677 cout << "\n test j" << j <<" " << i0 2678 << " " << i1 << " d01=" << d01 <<endl;} 2679 #endif 2680 assert (i0 >= nbvold); 2681 assert (i1 >= nbvold); 2682 assert(i0 != i1); 2646 if (i0< nbvold || i1<nbvold){ 2647 throw ErrorException(__FUNCT__,exprintf("i0< nbvold || i1<nbvold")); 2648 } 2683 2649 if (d01 == 0) 2684 2650 break; … … 2708 2674 = vertices[i0].r *w0 + vertices[i1].r* w1; 2709 2675 2710 #ifdef TRACETRIANGLE2711 if(trace) {2712 cout << "\n ref = "<< vertices[i0].ref << " " <<vertices[i1].ref <<endl;2713 }2714 #endif2715 2676 // update the new point points of the list 2716 2677 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) … … 2767 2728 2768 2729 quadtree->Add(*vi); // 2769 assert (tcvi->det >= 0) ;// internal 2730 if (tcvi->det<0){// internal 2731 throw ErrorException(__FUNCT__,exprintf("tcvi->det<0")); 2732 } 2770 2733 Add(*vi,tcvi,dete),NbSwap += vi->Optim(1); 2771 2734 } … … 2791 2754 if (first_np_or_next_t[kt]>0) 2792 2755 first_np_or_next_t[kt]=-Headt,Headt=kt; 2793 assert( ta.EdgeVertex(0) == s); 2756 if (ta.EdgeVertex(0) != s){ 2757 throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0) != s")); 2758 } 2794 2759 ta = Next(Adj(ta)); 2795 2760 } while ( (tbegin != (Triangle*) ta)); … … 3098 3063 // else if(mark[it] == -2 ) triangles[it].Draw(999); 3099 3064 it++;} // end white (it<nbt) 3100 assert(k== NbSubDomains); 3065 if (k!=NbSubDomains){ 3066 throw ErrorException(__FUNCT__,exprintf("k!=NbSubDomains")); 3067 } 3101 3068 if(OutSide) 3102 3069 { … … 3179 3146 // see routine MakeGeometricalEdgeToEdge 3180 3147 Edge &e = *GeometricalEdgetoEdge[Gh.Number(eg)]; 3181 assert(&e); 3148 if (&e==NULL){ 3149 throw ErrorException(__FUNCT__,exprintf("&e==NULL")); 3150 } 3182 3151 Vertex * v0 = e(0),*v1 = e(1); 3183 3152 Triangle *t = v0->t; … … 3186 3155 // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref 3187 3156 // << " ref edge =" << eg.ref << " sens " << sens ; 3188 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) 3189 sens = -sens ; 3157 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens ; 3190 3158 subdomains[i].sens = sens; 3191 3159 subdomains[i].edge = &e; 3192 // cout << " sens " << sens << " in geom " << eg[0].r << eg[1].r << " in mesh " << e[0].r << e[1].r << endl;3193 // cout << " v0 , v1 = " << Number(v0) << " " << Number(v1) << endl;3194 assert(t && sens);3160 if (!t || !sens){ 3161 throw ErrorException(__FUNCT__,exprintf("!t || !sens")); 3162 } 3195 3163 3196 3164 TriangleAdjacent ta(t,EdgesVertexTriangle[v0->vint][0]);// previous edges 3197 3165 3198 3166 while (1) { 3199 assert( v0 == ta.EdgeVertex(1) ); 3167 if ( v0 != ta.EdgeVertex(1) ){ 3168 throw ErrorException(__FUNCT__,exprintf("v0 != ta.EdgeVertex(1)")); 3169 } 3200 3170 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl; 3201 3171 if (ta.EdgeVertex(0) == v1) { // ok we find the edge … … 3225 3195 { 3226 3196 kkk++; 3227 assert(mark[Number(tt)]<0); 3197 if (mark[Number(tt)]>=0){ 3198 throw ErrorException(__FUNCT__,exprintf("mark[Number(tt)]>=0")); 3199 } 3228 3200 mark[Number(tt)]=i; 3229 3201 tt=tt->link; … … 3344 3316 { 3345 3317 t=t0=subdomains[i].head; 3346 assert(t0); // no empty sub domain 3318 if (!t0){ // not empty sub domain 3319 throw ErrorException(__FUNCT__,exprintf("!t0")); 3320 } 3347 3321 do { 3348 3322 Int4 kt = Number(t); 3349 assert(kt>=0 && kt < nbt ); 3350 assert(renu[kt]==-1); 3323 if (kt<0 || kt >= nbt ){ 3324 throw ErrorException(__FUNCT__,exprintf("kt<0 || kt >= nbt")); 3325 } 3326 if (renu[kt]!=-1){ 3327 throw ErrorException(__FUNCT__,exprintf("renu[kt]!=-1")); 3328 } 3351 3329 renu[kt]=k++; 3352 3330 } … … 3362 3340 3363 3341 // put the outside triangles at the end 3364 for ( it=0;it<nbt;it++) 3365 if (renu[it]==-1) 3366 renu[it]=k++; 3367 3368 assert(k == nbt); 3342 for ( it=0;it<nbt;it++){ 3343 if (renu[it]==-1) renu[it]=k++; 3344 } 3345 if (k != nbt){ 3346 throw ErrorException(__FUNCT__,exprintf("k != nbt")); 3347 } 3369 3348 // do the change on all the pointeur 3370 3349 for ( it=0;it<nbt;it++) … … 3396 3375 /*}}}1*/ 3397 3376 /*FUNCTION Triangles::ConsRefTriangle{{{1*/ 3398 Int4 Triangles::ConsRefTriangle(Int4 *reft) const {3377 Int4 Triangles::ConsRefTriangle(Int4* reft) const { 3399 3378 long int verbosity=0; 3400 assert(reft);3401 3379 register Triangle *t0,*t; 3402 3380 register Int4 k=0, num; 3403 for (Int4 it=0;it<nbt;it++) 3404 reft[it]=-1; // outside triangle 3405 for (Int4 i=0;i<NbSubDomains;i++) 3406 { 3381 for (Int4 it=0;it<nbt;it++) reft[it]=-1; // outside triangle 3382 for (Int4 i=0;i<NbSubDomains;i++){ 3407 3383 t=t0=subdomains[i].head; 3408 assert(t0); // no empty sub domain 3384 if (!t0){ // no empty sub domai{ 3385 throw ErrorException(__FUNCT__,exprintf("!t0")); 3386 } 3409 3387 // register Int4 color=i+1;// because the color 0 is outside triangle 3410 3388 do { k++; 3411 3389 num = Number(t); 3412 assert(num>=0 &&num < nbt); 3390 if (num<0 || num>=nbt){ 3391 throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt")); 3392 } 3413 3393 reft[num]=i; 3414 // cout << Number(t0) << " " <<Number(t)<< " " << i << endl;3415 3394 } 3416 3395 while (t0 != (t=t->link)); … … 3446 3425 1 on GeometricalEdge + abcisse 3447 3426 2 internal 3448 3449 3427 *************************************************************************/ 3450 assert(&BTh.Gh == &Gh); 3428 if (&BTh.Gh != &Gh){ 3429 throw ErrorException(__FUNCT__,exprintf("&BTh.Gh != &Gh")); 3430 } 3451 3431 3452 3432 long int verbosity=0; … … 3481 3461 throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); 3482 3462 } 3483 assert(vertices); 3484 for (i=0;i<Gh.nbv;i++) 3485 if (Gh[i].Required()) {//Gh vertices Required 3486 vertices[nbv] = Gh[i]; 3487 vertices[nbv].i = I2(0,0); 3488 Gh[i].to = vertices + nbv;// save Geom -> Th 3489 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); 3490 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl; 3491 nbv++;} 3492 else Gh[i].to=0; 3493 // 3494 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++) 3495 { 3463 if (vertices==NULL){ 3464 throw ErrorException(__FUNCT__,exprintf("vertices==NULL")); 3465 } 3466 for (i=0;i<Gh.nbv;i++){ 3467 if (Gh[i].Required()) {//Gh vertices Required 3468 vertices[nbv] = Gh[i]; 3469 vertices[nbv].i = I2(0,0); 3470 Gh[i].to = vertices + nbv;// save Geom -> Th 3471 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); 3472 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl; 3473 nbv++; 3474 } 3475 else Gh[i].to=0; 3476 } 3477 for (i=0;i<BTh.NbVerticesOnGeomVertex;i++){ 3496 3478 VertexOnGeom & vog = BTh.VerticesOnGeomVertex[i]; 3497 3479 if (vog.IsRequiredVertex()) { 3498 3480 GeometricalVertex * gv = vog; 3499 3481 Vertex *bv = vog; 3500 assert(gv->to);// use of Geom -> Th 3482 if (!gv->to){// use of Geom -> Th 3483 throw ErrorException(__FUNCT__,exprintf("!gv->to")); 3484 } 3501 3485 VertexOnBThVertex[NbVertexOnBThVertex++] = VertexOnVertex(gv->to,bv); 3502 3486 gv->to->m = bv->m; // for taking the metrix of the background mesh 3503 ;} 3487 ; 3488 } 3504 3489 } 3505 assert(NbVertexOnBThVertex == NbVerticesOnGeomVertex); 3490 if (NbVertexOnBThVertex != NbVerticesOnGeomVertex){ 3491 throw ErrorException(__FUNCT__,exprintf("NbVertexOnBThVertex != NbVerticesOnGeomVertex")); 3492 } 3506 3493 // new stuff FH with curve 3507 3494 // find the begin of the curve in BTh 3508 {3509 3495 Gh.UnMarkEdges(); 3510 3496 int bfind=0; … … 3543 3529 } 3544 3530 } 3545 assert( bfind==Gh.NbOfCurves); 3546 } 3531 if ( bfind!=Gh.NbOfCurves){ 3532 throw ErrorException(__FUNCT__,exprintf("bfind!=Gh.NbOfCurves")); 3533 } 3534 3547 3535 // method in 2 + 1 step 3548 3536 // 0.0) compute the length and the number of vertex to do allocation … … 3610 3598 VertexOnGeom *GA1; 3611 3599 Edge * PreviousNewEdge = 0; 3612 // cout << " --------------New Curve phase " << phase 3613 // << "---------- A0=" << *A0 << ei[k0] <<endl; 3614 assert (A0-vertices>=0 && A0-vertices <nbv); 3600 // New Curve phase 3601 if (A0-vertices<0 || A0-vertices>=nbv){ 3602 throw ErrorException(__FUNCT__,exprintf("A0-vertices<0 || A0-vertices>=nbv")); 3603 } 3615 3604 if(ongequi->Required() ) 3616 3605 { … … 3619 3608 } 3620 3609 else 3621 for(;;) 3622 { 3623 // assert(pe && BTh.Number(pe)>=0 && BTh.Number(pe)<=BTh.nbe); 3610 for(;;){ 3624 3611 Edge &ee=*pe; 3625 3612 Edge &eeequi=*peequi; … … 3627 3614 k1equi= 1 - k0equi; 3628 3615 3629 assert(pe && ee.on); 3616 if (!pe || !ee.on){ 3617 throw ErrorException(__FUNCT__,exprintf("!pe || !ee.on")); 3618 } 3630 3619 ee.on->SetMark(); 3631 3620 Vertex & v0=ee[0], & v1=ee[1]; … … 3636 3625 if (phase) {// computation of the new points 3637 3626 while ((i!=NbCreatePointOnCurve) && sNew <= L) { 3638 // cout << " L0= " << L0 << " L " << L << " sN=" 3639 // << sNew << " LAB=" << LAB << " NBPC =" <<NbCreatePointOnCurve<< " i " << i << endl; 3640 assert (sNew >= L0); 3641 assert(LAB); 3642 3643 3644 assert(vertices && nbv<nbvx); 3645 assert(edges && nbe < nbex); 3646 assert(VerticesOnGeomEdge && NbVerticesOnGeomEdge < NbVerticesOnGeomEdgex); 3627 if (sNew<L0){ 3628 throw ErrorException(__FUNCT__,exprintf("sNew<L0")); 3629 } 3630 if (!LAB){ 3631 throw ErrorException(__FUNCT__,exprintf("!LAB")); 3632 } 3633 if (!vertices || nbv>=nbvx){ 3634 throw ErrorException(__FUNCT__,exprintf("!vertices || nbv>=nbvx")); 3635 } 3636 if (!edges || nbe>=nbex){ 3637 throw ErrorException(__FUNCT__,exprintf("!edges || nbe>=nbex")); 3638 } 3639 if (!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex){ 3640 throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge || NbVerticesOnGeomEdge>=NbVerticesOnGeomEdgex")); 3641 } 3647 3642 // new vertex on edge 3648 3643 A1=vertices+nbv++; … … 3650 3645 Edge *e = edges + nbe++; 3651 3646 Real8 se= (sNew-L0)/LAB; 3652 assert(se>=0 && se < 1.000000001); 3647 if (se<0 || se>=1.000000001){ 3648 throw ErrorException(__FUNCT__,exprintf("se<0 || se>=1.000000001")); 3649 } 3653 3650 se = abscisseInterpole(v0.m,v1.m,AB,se,1); 3654 assert(se>=0 && se <= 1); 3651 if (se<0 || se>1){ 3652 throw ErrorException(__FUNCT__,exprintf("se<0 || se>1")); 3653 } 3655 3654 //((k1==1) != (k1==k1equi)) 3656 3655 se = k1 ? se : 1. - se; … … 3683 3682 3684 3683 } 3685 assert(ee.on->CurveNumber==ei.on->CurveNumber); 3684 if (ee.on->CurveNumber!=ei.on->CurveNumber){ 3685 throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber")); 3686 } 3686 3687 if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl; 3687 3688 if ( ee[k1].on->IsRequiredVertex()) { 3688 assert(eeequi[k1equi].on->IsRequiredVertex()); 3689 if (!eeequi[k1equi].on->IsRequiredVertex()){ 3690 throw ErrorException(__FUNCT__,exprintf("!eeequi[k1equi].on->IsRequiredVertex()")); 3691 } 3689 3692 register GeometricalVertex * GA1 = *eeequi[k1equi].on; 3690 3693 A1=GA1->to;// the vertex in new mesh 3691 assert (A1-vertices>=0 && A1-vertices <nbv); 3694 if (A1-vertices<0 || A1-vertices>=nbv){ 3695 throw ErrorException(__FUNCT__,exprintf("A1-vertices<0 || A1-vertices>=nbv")); 3696 } 3692 3697 break;} 3693 3698 if (!ee.adj[k1]) { … … 3717 3722 PreviousNewEdge = e; 3718 3723 3719 assert(i==NbCreatePointOnCurve); 3724 if (i!=NbCreatePointOnCurve){ 3725 throw ErrorException(__FUNCT__,exprintf("i!=NbCreatePointOnCurve")); 3726 } 3720 3727 3721 3728 } … … 3763 3770 } 3764 3771 } // for(step;;) 3765 assert(nbe); 3766 3772 if (nbe==0){ 3773 throw ErrorException(__FUNCT__,exprintf("nbe==0")); 3774 } 3767 3775 delete [] bcurve; 3768 3776 … … 3809 3817 nbv++; 3810 3818 } 3811 // assert( Gh.nbv < nbvx);3812 3819 3813 3820 // Method in 2 step: 0 and 1 … … 3815 3822 // 2) construct the edge 3816 3823 // generation of the curves 3817 assert(! edges); 3824 if (edges){ 3825 throw ErrorException(__FUNCT__,exprintf("edges")); 3826 } 3818 3827 // 2 step 3819 3828 // --step=0 to compute the number of edges + alloc at end … … 3849 3858 a=ei(0)->The(); 3850 3859 b=ei(1)->The(); 3851 assert(edges); 3860 if (!edges){ 3861 throw ErrorException(__FUNCT__,exprintf("!edges")); 3862 } 3852 3863 edges[nbe].v[0]=a->to; 3853 3864 edges[nbe].v[1]=b->to;; … … 3868 3879 NbNewPoints=0; 3869 3880 NbEdgeCurve=0; 3870 assert(nbvend < nbvx); 3881 if (nbvend>=nbvx){ 3882 throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx")); 3883 } 3871 3884 lcurve =0; 3872 3885 s = lstep; … … 3934 3947 else 3935 3948 kk0=kkk,ll0=llk;} 3936 assert(kk1 != kk0); 3949 if (kk1 == kk0){ 3950 throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); 3951 } 3937 3952 3938 3953 Real8 sbb = (ss-ll0 )/(ll1-ll0); … … 3972 3987 k = e->SensAdj[kprev];// next vertices 3973 3988 e = e->Adj[kprev]; 3974 assert(e); 3989 if (!e){ 3990 throw ErrorException(__FUNCT__,exprintf("!e")); 3991 } 3975 3992 }// for(;;) 3976 3993 vb = b->to; … … 4006 4023 } // for (i=0;i<nbe;i++) 4007 4024 if(!step) { 4008 // cout << "edges " << edges << " VerticesOnGeomEdge " <<VerticesOnGeomEdge << endl; 4009 assert(!edges); 4010 assert(!VerticesOnGeomEdge); 4025 if (edges){ 4026 throw ErrorException(__FUNCT__,exprintf("edges")); 4027 } 4028 if (VerticesOnGeomEdge){ 4029 throw ErrorException(__FUNCT__,exprintf("VerticesOnGeomEdge")); 4030 } 4011 4031 edges = new Edge[nbex=nbe]; 4012 4032 if(NbVerticesOnGeomEdge0) 4013 4033 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; 4014 assert(edges); 4015 assert(VerticesOnGeomEdge || NbVerticesOnGeomEdge0 ==0); 4034 if (!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0){ 4035 throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0")); 4036 } 4016 4037 // do the vertex on a geometrical vertex 4017 4038 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge; 4018 4039 } 4019 else 4020 assert(NbVerticesOnGeomEdge == NbVerticesOnGeomEdge0); 4021 // cout << " Nb of Curves = " << NbOfCurves << "nbe = " << nbe 4022 // << "== " << nbex << " nbv = " << nbv << endl; 4023 assert(nbex=nbe); 4040 else if (NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0){ 4041 throw ErrorException(__FUNCT__,exprintf("NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0")); 4042 } 4024 4043 } // for (step=0;step<2;step++) 4025 4044 … … 4035 4054 /*FUNCTION Triangles::MakeGeometricalEdgeToEdge{{{1*/ 4036 4055 Edge** Triangles::MakeGeometricalEdgeToEdge() { 4037 assert(Gh.nbe); 4056 if (!Gh.nbe){ 4057 throw ErrorException(__FUNCT__,exprintf("!Gh.nbe")); 4058 } 4038 4059 Edge **e= new (Edge* [Gh.nbe]); 4039 4060 … … 4072 4093 cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber 4073 4094 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl; 4074 // assert( e[i]);4075 4095 } 4076 4096 if(kk) throw ErrorException(__FUNCT__,exprintf("See above")); … … 4096 4116 pmax = pmax+DD; 4097 4117 coefIcoor= (MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y)); 4098 assert(coefIcoor >0); 4118 if (coefIcoor<=0){ 4119 throw ErrorException(__FUNCT__,exprintf("coefIcoor<=0")); 4120 } 4099 4121 4100 4122 // generation of integer coord … … 4142 4164 4143 4165 //initialize ordre 4144 assert(ordre); 4166 if (!ordre){ 4167 throw ErrorException(__FUNCT__,exprintf("!ordre")); 4168 } 4145 4169 for (i=0;i<nbv;i++) ordre[i]=0; 4146 4170 … … 4174 4198 } 4175 4199 else if(st[k]>=0) { 4176 assert( ! triangles[i].TriangleAdj(j) && !triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))); 4200 if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){ 4201 throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))")); 4202 } 4177 4203 4178 4204 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); … … 4329 4355 Vertex *v1= ta.EdgeVertex(1); 4330 4356 Int4 k =edge4->addtrie(v0?Number(v0):nbv,v1? Number(v1):nbv); 4331 assert(st[k] >=0); 4357 if (st[k]<0){ 4358 throw ErrorException(__FUNCT__,exprintf("st[k]<0")); 4359 } 4332 4360 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); 4333 4361 ta.SetLock(); … … 4348 4376 } 4349 4377 } 4350 // cout << savenbt+NbTfillHoll << " " << savenbtx << endl; 4351 assert(savenbt+NbTfillHoll <= savenbtx ); 4378 if (savenbt+NbTfillHoll>savenbtx ){ 4379 throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx")); 4380 } 4352 4381 // copy of the outside triangles in saveTriangles 4353 4382 for (i=0;i<nbt;i++){ … … 4442 4471 { 4443 4472 NbSwap++; 4444 assert(k++<20000); 4473 k++; 4474 if (k>=20000){ 4475 throw ErrorException(__FUNCT__,exprintf("k>=20000")); 4476 } 4445 4477 t= tp->at[jp]; // set unchange t qnd j for previous triangles 4446 4478 j= NextEdge[tp->aa[jp]&3]; … … 4641 4673 /*FUNCTION Triangles::Crack{{{1*/ 4642 4674 int Triangles::Crack() { 4643 assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 4644 for (int i=0;i<NbCrackedEdges;i++) 4645 CrackedEdges[i].Crack(); 4675 if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{ 4676 throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges!=0 && NbCrackedVertices<=0")); 4677 } 4678 for (int i=0;i<NbCrackedEdges;i++) CrackedEdges[i].Crack(); 4646 4679 return NbCrackedEdges; 4647 4680 } … … 4649 4682 /*FUNCTION Triangles::UnCrack{{{1*/ 4650 4683 int Triangles::UnCrack() { 4651 assert(NbCrackedEdges ==0 || NbCrackedVertices >0); 4684 if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{ 4685 throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges ==0 || NbCrackedVertices >0")); 4686 } 4652 4687 for (int i=0;i<NbCrackedEdges;i++) 4653 4688 CrackedEdges[i].UnCrack(); … … 4703 4738 Triangle * tbegin = v.t; 4704 4739 int i = v.vint; 4705 assert(tbegin && (i >= 0 ) && (i <3)); 4740 if (!tbegin || (i<0) || (i>=3)){ 4741 throw ErrorException(__FUNCT__,exprintf("!tbegin || (i<0) || (i>=3)")); 4742 } 4706 4743 // turn around the vertex v 4707 4744 TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge … … 4714 4751 { 4715 4752 TriangleAdjacent tta=(ta.Adj()); 4716 assert(tta.Cracked()); 4753 if (!tta.Cracked()){ 4754 throw ErrorException(__FUNCT__,exprintf("!tta.Cracked()")); 4755 } 4717 4756 if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart 4718 4757 if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use … … 4721 4760 if ( tt->link ) { // if good triangles store the value 4722 4761 int it = Number(tt); 4723 assert(it < nt); 4762 if (it>=nt){ 4763 throw ErrorException(__FUNCT__,exprintf("(it>=nt)")); 4764 } 4724 4765 (*tt)(kv)= vv; // Change the vertex of triangle 4725 4766 if(vv<vend) {*vv= v;vv->ReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref … … 4732 4773 ta = Next(ta).Adj(); 4733 4774 } while ( (tbegin != ta)); 4734 assert(k); 4775 if (!k){ 4776 throw ErrorException(__FUNCT__,exprintf("!k")); 4777 } 4735 4778 if (kc) nbcrakev++; 4736 4779 } … … 4799 4842 Triangle * t=0; 4800 4843 int j,jp,jn,jj; 4801 if (tstart) 4802 t=tstart;4803 else4804 {4805 assert(quadtree);4844 if (tstart) t=tstart; 4845 else { 4846 if (!quadtree){ 4847 throw ErrorException(__FUNCT__,exprintf("!quadtree")); 4848 } 4806 4849 Vertex *a = quadtree->NearestVertex(B.x,B.y) ; 4807 4850 … … 4812 4855 throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening")); 4813 4856 } 4814 assert(a>= vertices && a < vertices+nbv); 4815 // int k=0; 4857 if (a<vertices || a>=vertices+nbv){ 4858 throw ErrorException(__FUNCT__,exprintf("a<vertices || a>=vertices+nbv")); 4859 } 4816 4860 t = a->t; 4817 assert(t>= triangles && t < triangles+nbt); 4818 4819 } 4861 if (t<triangles || t>=triangles+nbt){ 4862 throw ErrorException(__FUNCT__,exprintf("t<triangles || t>=triangles+nbt")); 4863 } 4864 } 4820 4865 Icoor2 detop ; 4821 4866 int kkkk =0; // number of test triangle 4822 4867 4823 while ( t->det < 0) 4824 { // the initial triangles is outside 4868 while ( t->det < 0){ // the initial triangles is outside 4825 4869 int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 4826 assert(k0>=0); // k0 the NULL vertex 4870 if (k0<0){ // k0 the NULL vertex 4871 throw ErrorException(__FUNCT__,exprintf("k0<0")); 4872 } 4827 4873 int k1=NextVertex[k0],k2=PreviousVertex[k0]; 4828 4874 dete[k0]=det(B,(*t)[k1],(*t)[k2]); … … 4831 4877 return t; 4832 4878 t = t->TriangleAdj(OppositeEdge[k0]); 4833 assert(kkkk++ < 2); 4879 kkkk++; 4880 if (kkkk>=2){ 4881 throw ErrorException(__FUNCT__,exprintf("kkkk>=2")); 4882 } 4834 4883 } 4835 4884 … … 4837 4886 detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B); 4838 4887 4839 while(t->det > 0 ) 4840 { 4841 assert( kkkk++ < 2000 ); 4888 while(t->det > 0 ) { 4889 kkkk++; 4890 if (kkkk>=2000){ 4891 throw ErrorException(__FUNCT__,exprintf("kkkk>=2000")); 4892 } 4842 4893 j= OppositeVertex[jj]; 4843 4894 dete[j] = detop; //det(*b,*s1,*s2); … … 4856 4907 // 2 => two way go in way 1 or 2 randomly 4857 4908 4858 if (k==0) 4859 break;4860 if ( k == 2 && BinaryRand())4861 Exchange(ii[0],ii[1]);4862 assert ( k < 3);4909 if (k==0) break; 4910 if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]); 4911 if ( k>=3){ 4912 throw ErrorException(__FUNCT__,exprintf("k>=3")); 4913 } 4863 4914 TriangleAdjacent t1 = t->Adj(jj=ii[0]); 4864 4915 if ((t1.det() < 0 ) && (k == 2)) … … 4887 4938 Real8 hmin = Gh.MinimalHmin(); 4888 4939 Real8 maxaniso = 1e6; 4889 assert(hmax>0); 4940 if (hmax<=0){ 4941 throw ErrorException(__FUNCT__,exprintf("hmax<=0")); 4942 } 4890 4943 SetVertexFieldOn(); 4891 4944 if (errC > 1) errC = 1; … … 4902 4955 Real8 s = GV; 4903 4956 R2 tg; 4904 // cerr << i << " " << j << " " << Number(V) << " on = "4905 // << Gh.Number(eg) << " at s = " << s << " " << endl;4906 4957 Real8 R1= eg->R1tg(s,tg); 4907 // cerr << " R = " << 1/Max(R1,1e-20) << tg << " on x "4908 // << V.r << errC/ Max(R1,1e-20) << " hold=" <<V.m(tg) << " " << endl;4909 4958 Real8 ht = hmax; 4910 if (R1>1.0e-20) 4911 { // err relative to the length of the edge 4959 if (R1>1.0e-20) { // err relative to the length of the edge 4912 4960 ht = Min(Max(errC/R1,hmin),hmax); 4913 4961 } 4914 4962 Real8 hn = iso? ht : Min(hmax,ht*maxaniso); 4915 //cerr << ht << " " << hn << "m=" << edges[i][j].m << endl; 4916 assert(ht>0 && hn>0); 4963 if (ht<=0 || hn<=0){ 4964 throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); 4965 } 4917 4966 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); 4918 4967 //cerr << " : " ; … … 5217 5266 Real8 cAB = det3x3(taa[0],taa[1],bb); 5218 5267 5219 assert(det33); 5220 // det33=1; 5221 // verif 5222 // cout << " " << (taa[0][0]*cBC + taa[1][0]*cCA + taa[2][0] * cAB)/det33 << " == " << bb[0] ; 5223 // cout << " " << (taa[0][1]*cBC + taa[1][1]*cCA + taa[2][1] * cAB)/det33 << " == " << bb[1]; 5224 // cout << " " << (taa[0][2]*cBC + taa[1][2]*cCA + taa[2][2] * cAB)/det33 << " == " << bb[2] 5225 // << " -- " ; 5226 //cout << lla*sA + llb*sB+llc*sC+ (lla*llb* cAB + llb*llc* cBC + llc*lla*cCA)/det33 5227 // << " == " << llf << endl; 5268 if (!det33){ 5269 throw ErrorException(__FUNCT__,exprintf("!det33")); 5270 } 5228 5271 // computation of the gradient in the element 5229 5272 … … 5514 5557 while (Head0>=0&& kk++<100) { 5515 5558 kch=0; 5516 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) 5517 {// pour tous les triangles autour du sommet s5559 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) { 5560 // pour tous les triangles autour du sommet s 5518 5561 // cout << kk << " i = " << i << " " << ip << endl; 5519 register Triangle * t= vertices[i].t; 5520 assert(t); 5562 register Triangle* t= vertices[i].t; 5563 if (!t){ 5564 throw ErrorException(__FUNCT__,exprintf("!t")); 5565 } 5521 5566 Vertex & vi = vertices[i]; 5522 5567 TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]); 5523 5568 Vertex *pvj0 = ta.EdgeVertex(0); 5524 5569 while (1) { 5525 // cout << i << " " << Number(ta.EdgeVertex(0)) << " "5526 // << Number(ta.EdgeVertex(1)) << " ---> " ;5527 5570 ta=Previous(Adj(ta)); 5528 // cout << Number(ta.EdgeVertex(0)) << " " << Number(ta.EdgeVertex(1)) << endl; 5529 assert(vertices+i == ta.EdgeVertex(1)); 5571 if (vertices+i != ta.EdgeVertex(1)){ 5572 throw ErrorException(__FUNCT__,exprintf("vertices+i != ta.EdgeVertex(1)")); 5573 } 5530 5574 Vertex & vj = *(ta.EdgeVertex(0)); 5531 5575 if ( &vj ) { 5532 5576 j= &vj-vertices; 5533 assert(j>=0 && j < nbv); 5577 if (j<0 || j >= nbv){ 5578 throw ErrorException(__FUNCT__,exprintf("j<0 || j >= nbv")); 5579 } 5534 5580 R2 Aij = (R2) vj - (R2) vi; 5535 5581 Real8 ll = Norme2(Aij); … … 5611 5657 if (inbvx) { 5612 5658 vertices=new Vertex[nbvx]; 5613 assert(vertices); 5659 if (!vertices){ 5660 throw ErrorException(__FUNCT__,exprintf("!vertices")); 5661 } 5614 5662 ordre=new (Vertex* [nbvx]); 5615 assert(ordre); 5663 if (!ordre){ 5664 throw ErrorException(__FUNCT__,exprintf("!ordre")); 5665 } 5616 5666 triangles=new Triangle[nbtx]; 5617 assert(triangles);} 5667 if (!triangles){ 5668 throw ErrorException(__FUNCT__,exprintf("!triangles")); 5669 } 5670 } 5618 5671 else { 5619 5672 vertices=0; … … 5622 5675 nbtx=0; 5623 5676 } 5624 if ( name || inbvx) 5625 { 5677 if ( name || inbvx) { 5626 5678 time_t timer =time(0); 5627 5679 char buf[70]; … … 5713 5765 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){ 5714 5766 CurrentTh=&Th; 5715 assert(&Th);5716 5767 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 5717 5768 Icoor2 dete[3]; … … 5738 5789 k = ta; 5739 5790 Exchange(aa,bb); 5740 assert(tc->link); 5791 if (!tc->link){ 5792 throw ErrorException(__FUNCT__,exprintf("!tc->link")); 5793 } 5741 5794 } 5742 5795 a[VerticesOfTriangularEdge[k][0]] = aa; … … 5751 5804 int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 5752 5805 int dir=0; 5753 assert(k>=0); 5806 if (k<0){ 5807 throw ErrorException(__FUNCT__,exprintf("k<0")); 5808 } 5754 5809 int kkk=0; 5755 5810 Icoor2 IJ_IA,IJ_AJ; 5756 5811 TriangleAdjacent edge(t,OppositeEdge[k]); 5757 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) 5758 { 5759 5760 assert(kkk++<1000); 5812 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { 5813 kkk++; 5814 if (kkk>=1000){ 5815 throw ErrorException(__FUNCT__,exprintf("kkk>=1000")); 5816 } 5761 5817 Vertex &vI = *edge.EdgeVertex(0); 5762 5818 Vertex &vJ = *edge.EdgeVertex(1); … … 5774 5830 continue;}}// go in direction j 5775 5831 double IJ2 = IJ_IA + IJ_AJ; 5776 assert(IJ2); 5832 if (IJ2==0){ 5833 throw ErrorException(__FUNCT__,exprintf("IJ2==0")); 5834 } 5777 5835 a= IJ_AJ/IJ2; 5778 5836 b= IJ_IA/IJ2; … … 5790 5848 // restart: 5791 5849 // int dir=0; 5792 assert(t->link == 0); 5850 if (t->link != 0){ 5851 throw ErrorException(__FUNCT__,exprintf("t->link != 0")); 5852 } 5793 5853 // to have a starting edges 5794 5854 // try the 3 edge bourna-- in case of internal hole … … 5864 5924 } 5865 5925 } 5866 assert(cas !=-2); 5926 if (cas ==-2){ 5927 throw ErrorException(__FUNCT__,exprintf("cas==-2")); 5928 } 5867 5929 // l1 = ||C s1|| , l0 = ||C s0|| 5868 5930 // where s0,s1 are the vertex of the edge er … … 5877 5939 5878 5940 Triangle * tt=t=edge=Adj(Previous(edge)); 5879 do { // loop around vertex s 5880 5881 assert(edge.EdgeVertex(0)==s && kkk++<10000); 5941 do { // loop over vertex s 5942 kkk++; 5943 if (edge.EdgeVertex(0)!=s && kkk>=10000){ 5944 throw ErrorException(__FUNCT__,exprintf("edge.EdgeVertex(0)!=s && kkk>=10000")); 5945 } 5882 5946 5883 5947 int link = tt->link == 0; … … 5909 5973 } while (t!=tt); 5910 5974 5911 assert((Triangle *) er); 5975 if (!(Triangle *) er){ 5976 throw ErrorException(__FUNCT__,exprintf("!(Triangle *) er")); 5977 } 5912 5978 I2 A((I2)*er.EdgeVertex(0)); 5913 5979 I2 B((I2)*er.EdgeVertex(1)); … … 5980 6046 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 5981 6047 Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 5982 assert ( a1 >= 0 && a1 < 3 ); 6048 if ( a1<0 || a1>=3 ){ 6049 throw ErrorException(__FUNCT__,exprintf("a1<0 || a1>=3")); 6050 } 5983 6051 5984 6052 Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; … … 5990 6058 Icoor2 det1=t1->det , det2=t2->det ; 5991 6059 Icoor2 detT = det1+det2; 5992 assert((det1>0 ) && (det2 > 0)); 5993 assert ( (detsa < 0) && (detsb >0) ); // [a,b] cut infinite line va,bb 6060 if ((det1<=0 ) || (det2<=0)){ 6061 throw ErrorException(__FUNCT__,exprintf("(det1<=0 ) || (det2<=0)")); 6062 } 6063 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 6064 throw ErrorException(__FUNCT__,exprintf("(detsa>=0) || (detsb<=0)")); 6065 } 5994 6066 Icoor2 ndet1 = bamg::det(s1,sa,s2); 5995 6067 Icoor2 ndet2 = detT - ndet1; … … 6052 6124 int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) { 6053 6125 int NbSwap =0; 6054 assert(a.t && b.t); // the 2 vertex is in a mesh 6126 if (!a.t || !b.t){ // the 2 vertex is in a mesh 6127 throw ErrorException(__FUNCT__,exprintf("!a.t || !b.t")); 6128 } 6055 6129 int k=0; 6056 6130 taret=TriangleAdjacent(0,0); // erreur … … 6067 6141 v2 = tta.EdgeVertex(0); 6068 6142 vbegin =v2; 6069 assert(v2); 6143 if (!v2){ 6144 throw ErrorException(__FUNCT__,exprintf("!v2")); 6145 } 6070 6146 det2 = det(*v2,a,b); 6071 6147 // cout << " No Change try the next" << endl; … … 6083 6159 Vertex * va = &a, *vb = &b; 6084 6160 tc = Previous(tc); 6085 assert ( v1 && v2); 6161 if (!v1 || !v2){ 6162 throw ErrorException(__FUNCT__,exprintf("!v1 || !v2")); 6163 } 6086 6164 Icoor2 detss = 0,l=0,ks; 6087 6165 // cout << "Real ForcingEdge " << *va << *vb << detss << endl; … … 6105 6183 } 6106 6184 tta = tc; 6107 assert(k++<2000); 6185 k++; 6186 if (k>=2000){ 6187 throw ErrorException(__FUNCT__,exprintf("k>=2000")); 6188 } 6108 6189 if ( vbegin == v2 ) return -1;// error 6109 6190 } -
issm/trunk/src/c/Bamgx/objects/Vertex.cpp
r2850 r2854 28 28 register int k=0,kk=0,j = EdgesVertexTriangle[vint][0],jc; 29 29 R2 P(s->r),PNew(0,0); 30 // cout << BTh.quadtree << " " << BTh.quadtree->root << endl;31 // assert(BTh.quadtree && BTh.quadtree->root);32 30 do { 33 31 k++; … … 46 44 tria = ttc; 47 45 j = NextEdge[jc]; 48 assert(k<2000); 46 if (k>=2000){ 47 throw ErrorException(__FUNCT__,exprintf("k>=2000")); 48 } 49 49 } while ( tbegin != tria); 50 50 if (kk<4) return 0; … … 117 117 tria = ttc; 118 118 j = NextEdge[jc]; 119 assert(k<2000); 119 if (k>=2000){ 120 throw ErrorException(__FUNCT__,exprintf("k>=2000")); 121 } 120 122 } while ( tbegin != tria); 121 123 if (ok && loop) vP=vPsave; // no move
Note:
See TracChangeset
for help on using the changeset viewer.
