![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2858
- Timestamp:
- 01/15/10 14:05:28 (15 years ago)
- Location:
- issm/trunk/src/c/Bamgx
- Files:
-
- 5 edited
-
Bamgx.cpp (modified) (1 diff)
-
Mesh2.h (modified) (1 diff)
-
Metric.h (modified) (1 diff)
-
objects/Geometry.cpp (modified) (1 diff)
-
objects/Triangles.cpp (modified) (94 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Bamgx.cpp
r2855 r2858 69 69 // some verification 70 70 if ( maxsubdiv > boundmaxsubdiv || maxsubdiv <= 1.0){ 71 cerr << " -maxsubdiv " << maxsubdiv << " is not in ]1,"<< boundmaxsubdiv << "]" << endl; 72 exit(3); 71 throw ErrorException(__FUNCT__,exprintf("maxsubdiv %g is should be in ]1,%g]",maxsubdiv,boundmaxsubdiv)); 73 72 } 74 73 if (iso) anisomax=1; -
issm/trunk/src/c/Bamgx/Mesh2.h
r2856 r2858 815 815 void NewPoints( Triangles &,int KeepBackVertex =1 ); 816 816 Int4 InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) ; 817 void NewPointsOld( Triangles & );818 817 void NewPoints(int KeepBackVertex=1){ NewPoints(*this,KeepBackVertex);} 819 818 void ReNumberingTheTriangleBySubDomain(bool justcompress=false); -
issm/trunk/src/c/Bamgx/Metric.h
r2837 r2858 179 179 180 180 181 inline void MetricAnIso::Box(Real4 &hx,Real4 &hy) const 182 { 181 inline void MetricAnIso::Box(Real4 &hx,Real4 &hy) const { 183 182 Real8 d= a11*a22-a21*a21; 184 183 hx = sqrt(a22/d); 185 184 hy = sqrt(a11/d); 186 // cerr << " hx = " << hx << " hy = " << hy << endl;187 185 } 188 186 -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r2856 r2858 708 708 GeometricalEdge & e1=edges[i]; 709 709 GeometricalEdge & e2=*e1.link; 710 cerr << i << " " << e1[0].The() << " " << e2[0].The() << " " << e1[1].The() << " " << e2[1].The() << endl;711 710 if ( e1[0].The() == e2[0].The() && e1[1].The() == e2[1].The() ) 712 711 { -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2857 r2858 1546 1546 } // tt 1547 1547 else 1548 cerr <<endl << " Bug " <<i<< " " << j << " t=" << t << endl; 1549 1548 throw ErrorException(__FUNCT__,exprintf("Bug...")); 1550 1549 } // ke<0 1551 1550 else … … 2349 2348 } 2350 2349 /*}}}1*/ 2351 /*FUNCTION Triangles::NewPointsOld{{{1*/2352 void Triangles::NewPointsOld(Triangles & Bh) {2353 long int verbosity=2;2354 Real8 seuil= 0.7 ;// for two neart point2355 if (verbosity>1)2356 cout << " begin : Triangles::NewPointsOld " << endl;2357 Int4 i,k;2358 int j;2359 Int4 BeginNewPoint[3];2360 Int4 EndNewPoint[3];2361 int step[3];2362 Int4 *first_np_or_next_t = new Int4[nbtx];2363 Int4 ColorEdge[3];2364 Int4 color=-1;2365 Triangle *t;2366 // generation of the list of next Triangle2367 // at 1 time we test all the triangles2368 Int4 Headt =0,next_t;2369 for(i=0;i<nbt;i++)2370 first_np_or_next_t[i]=-(i+1);2371 // end list i >= nbt2372 // the list of test triangle is2373 // the next Triangle on i is -first_np_or_next_t[i]2374 Int4 nbtold,nbvold;2375 2376 // Big loop2377 do {2378 nbtold = nbt;2379 nbvold = nbv;2380 // default size of IntersectionTriangle2381 2382 i=Headt;2383 next_t=-first_np_or_next_t[i];2384 for(t=&triangles[i];i<nbt;t=&triangles[i=next_t],next_t=-first_np_or_next_t[i]) { // for each triangle t2385 // we can change first_np_or_next_t[i]2386 if (i<0 || i >= nbt ){2387 throw ErrorException(__FUNCT__,exprintf("i<0 || i >= nbt"));2388 }2389 first_np_or_next_t[i] = nbv; // to save the fist new point of triangle2390 for(j=0;j<3;j++)2391 { // for each edge2392 TriangleAdjacent tj(t,j);2393 // color++;// the color is 3i+j2394 color = 3*i + j ;;2395 ColorEdge[j]=color;2396 BeginNewPoint[j]=nbv;2397 EndNewPoint[j]=nbv-1;2398 step[j]=1;// right sens2399 Vertex & vA = * tj.EdgeVertex(0);2400 Vertex & vB = * tj.EdgeVertex(1);2401 2402 if (!t->link) continue;// boundary2403 if (t->det <0) continue;2404 if (t->Locked(j)) continue;2405 2406 TriangleAdjacent tadjj = t->Adj(j);2407 Triangle * ta= tadjj;2408 2409 if (ta->det <0) continue;2410 2411 R2 A = vA;2412 R2 B = vB;2413 2414 k=Number(ta);2415 // the 2 opposite vertices2416 const Vertex & vC1 = *tj.OppositeVertex();2417 const Vertex & vC2 = *tadjj.OppositeVertex();2418 2419 if(first_np_or_next_t[k]>0) { // this edge is done before2420 // find the color of the edge and begin , end of newpoint2421 register int kk = t->NuEdgeTriangleAdj(j);2422 if ((*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])){2423 throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][0]) != (*ta)(VerticesOfTriangularEdge[kk][1])"));2424 }2425 if ((*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])){2426 throw ErrorException(__FUNCT__,exprintf("(*t)(VerticesOfTriangularEdge[j][1]) != (*ta)(VerticesOfTriangularEdge[kk][0])"));2427 }2428 register Int4 kolor =3*k + kk;2429 ColorEdge[j]=kolor;2430 register Int4 kkk= 1;2431 step[j]=-1;// other sens2432 BeginNewPoint[j]=0;2433 EndNewPoint[j]=-1; // empty list2434 for (Int4 iv=first_np_or_next_t[k];iv<nbv;iv++)2435 if (vertices[iv].color > kolor)2436 break; // the color is passed2437 else if (vertices[iv].color == kolor) {2438 EndNewPoint[j]=iv;2439 if (kkk) // one time test2440 kkk=0,BeginNewPoint[j]=iv;}2441 continue; // next edge of the triangle2442 } // end if( k < i)2443 2444 2445 const Int4 NbvOld = nbv;2446 lIntTria.SplitEdge(Bh,A,B);2447 // Int4 NbvNp =2448 lIntTria.NewPoints(vertices,nbv,nbvx);2449 Int4 nbvNew=nbv;2450 nbv = NbvOld;2451 for (Int4 iv=NbvOld;iv<nbvNew;iv++) {2452 vertices[nbv].color = color;2453 vertices[nbv].ReferenceNumber=nbv;// circular link2454 R2 C = vertices[iv].r;2455 vertices[nbv].r = C;2456 vertices[nbv].m = vertices[iv].m ;2457 // test if the new point is not to close to the 2 opposite vertex2458 R2 CC1 = C-vC1 , CC2 = C-vC2;2459 if ( ( (vC1.m(CC1) + vertices[nbv].m(CC1)) > seuil)2460 && ( (vC2.m(CC2) + vertices[nbv].m(CC2)) > seuil) )2461 nbv++;2462 }2463 2464 EndNewPoint[j] = nbv-1;2465 } // end loop for each edge2466 2467 if (!t->link) continue;// boundary2468 if (t->det<=0) continue;// outside2469 // continue;2470 for(int j0=0;j0<3;j0++)2471 for (Int4 i0= BeginNewPoint[j0]; i0 <= EndNewPoint[j0];i0++)2472 {2473 // find the neart point one the opposite edge2474 // to compute i12475 Vertex & vi0 = vertices[i0];2476 int kstack = 0;2477 Int4 stack[10];2478 // Int4 savRef[10];2479 int j1 = j0;2480 while (j0 != (j1 = NextEdge[j1])) {//loop on the 2 other edge2481 // computation of the intersection of edge j1 and DOrto2482 // take the good sens2483 2484 if (BeginNewPoint[j1]> EndNewPoint[j1])2485 continue; //2486 else if (EndNewPoint[j1] - BeginNewPoint[j1] <1) {2487 for (Int4 ii1= BeginNewPoint[j1];ii1<=EndNewPoint[j1];ii1++)2488 stack[kstack++] = ii1;2489 continue;}2490 2491 2492 int k0,k1;2493 if (step[j1]<0) k0=1,k1=0; // reverse2494 else k0=0,k1=1;2495 R2 V10 = (R2)(*t)[VerticesOfTriangularEdge[j1][k0]];2496 R2 V11 = (R2)(*t)[VerticesOfTriangularEdge[j1][k1]];2497 R2 D = V11-V10;2498 Real8 c0 = vi0.m(D,(R2) vi0);2499 2500 Real8 c10 = vi0.m(D,V10);2501 Real8 c11 = vi0.m(D,V11);2502 2503 Real8 s;2504 //cout << " --i0 = " << i0 << D << V10 << V11 << endl ;2505 //cout << " c10 " << c10 << " c0 " << c0 << " c11 " << c11 << endl;2506 if (( c10 < c0 ) && (c0 < c11))2507 s = (c11-c0)/(c11-c10);2508 else if (( c11 < c0 ) && (c0 < c10))2509 s = (c11-c0) /(c11-c10);2510 else break;2511 R2 VP = V10*s + V11*(1-s);2512 int sss = (c11-c10) >0 ? 1 : -1;2513 // find the 2 point by dichotomie2514 //cout << " t =" << Number(t) << " c0 " << c0 ;2515 Int4 ii0 = BeginNewPoint[j1];2516 Int4 ii1 = EndNewPoint[j1];2517 Real8 ciii=-1,cii0=-1,cii1=-1 ;2518 if ( sss * ((cii0=vi0.m(D,(R2) vertices[ii0]))- c0) >0 )2519 stack[kstack++] = ii0;//cout << " add+0 " << ii0;2520 else if ( sss * ((cii1= vi0.m(D ,(R2) vertices[ii1]))- c0) < 0 )2521 stack[kstack++] = ii1;//cout << " add+1 " << ii1;2522 else {2523 while ((ii1-ii0)> 1) {2524 Int4 iii = (ii0+ii1)/2;2525 ciii = vi0.m( D ,(R2) vertices[iii]);2526 //cout << " (iii = " << iii << " " << ciii << ") ";2527 if ( sss * (ciii - c0) <0 ) ii0 = iii;2528 else ii1 = iii;}2529 stack[kstack++] = ii0;// cout << " add0 " << ii0;2530 if (ii1 != ii0) stack[kstack++] = ii1;//cout << " add1 " << ii1;2531 }2532 if (kstack >5) // bug ?2533 cout << "NewPoints: bug????? " << kstack << " stack " << stack[kstack]<< endl;2534 }2535 2536 stack[kstack++] = -1; // to stop2537 Int4 i1;2538 kstack =0;2539 while( (i1=stack[kstack++]) >= 0){ // the two parameter is i0 and i12540 if (i1 >= nbv || i1 < 0){2541 throw ErrorException(__FUNCT__,exprintf("i1 >= nbv || i1 < 0"));2542 }2543 if (i0 >= nbv || i0 < 0){2544 throw ErrorException(__FUNCT__,exprintf("i0 >= nbv || i0 < 0"));2545 }2546 if (i1 == i0){2547 throw ErrorException(__FUNCT__,exprintf("i1 == i0"));2548 }2549 R2 v01 = (R2) vertices[i1]- (R2) vertices[i0];2550 Real8 d01 = (vertices[i0].m(v01) + vertices[i1].m(v01));2551 if (i0< nbvold || i1<nbvold){2552 throw ErrorException(__FUNCT__,exprintf("i0< nbvold || i1<nbvold"));2553 }2554 if (d01 == 0)2555 break;2556 if ( d01 < seuil)2557 if (i1<nbvold) {2558 // remove all the points i0;2559 register Int4 ip,ipp;2560 for (ip=i0;i0 != (ipp = vertices[ip].ReferenceNumber);ip=ipp)2561 vertices[ip].ReferenceNumber = -1;// mark remove2562 vertices[ip].ReferenceNumber = -1;// mark remove2563 }2564 else {2565 // remove on of two points2566 register Int4 ip0, ip1, ipp0,ipp1;2567 register int kk0=1,kk1=1;2568 // count the number of common points to compute weight w0,w12569 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0) kk0++;2570 for (ip1=i1;i1 != (ipp1 = vertices[ip1].ReferenceNumber);ip1=ipp1) kk1++;2571 2572 register Real8 w0 = ((Real8) kk0)/(kk0+kk1);2573 register Real8 w1 = ((Real8) kk1)/(kk0+kk1);2574 2575 // make a circular link2576 Exchange(vertices[i0].ReferenceNumber,vertices[i1].ReferenceNumber);2577 // the new coordinate2578 R2 C //= vertices[i0] ;2579 = vertices[i0].r *w0 + vertices[i1].r* w1;2580 2581 // update the new point points of the list2582 for (ip0=i0;i0 != (ipp0 = vertices[ip0].ReferenceNumber);ip0=ipp0)2583 vertices[ip0].r = C;2584 vertices[ip0].r = C;2585 }2586 }2587 } // for (i0= ....2588 }// for triangle2589 2590 // remove of all the double points2591 2592 Int4 ip,ipp,kkk=nbvold;2593 for (i=nbvold;i<nbv;i++)2594 if (vertices[i].ReferenceNumber>=0) {// good points2595 // cout <<" i = " << i ;2596 for (ip=i;i != (ipp = vertices[ip].ReferenceNumber);ip=ipp)2597 vertices[ip].ReferenceNumber = -1;// mark remove2598 vertices[ip].ReferenceNumber = -1;// mark remove2599 // cout << i << " ---> " << kkk << endl;2600 vertices[kkk] = vertices[i];2601 vertices[kkk].i = toI2(vertices[kkk].r);2602 vertices[kkk++].ReferenceNumber = 0;2603 2604 }2605 2606 // insertion part ---2607 2608 const Int4 nbvnew = kkk-nbvold;2609 2610 cout << " Remove " << nbv - kkk << " to close vertex " ;2611 cout << " and Insert the " <<nbvnew<< " new points " << endl;2612 nbv = kkk;2613 Int4 NbSwap=0;2614 Icoor2 dete[3];2615 2616 // construction d'un ordre aleatoire2617 if (! nbvnew)2618 break;2619 if (nbvnew) {2620 const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ;2621 Int4 k3 = rand()%nbvnew ;2622 for (Int4 is3=0; is3<nbvnew; is3++)2623 ordre[nbvold+is3]= &vertices[nbvold +(k3 = (k3 + PrimeNumber)% nbvnew)];2624 2625 for (i=nbvold;i<nbv;i++)2626 { Vertex * vi = ordre[i];2627 Triangle *tcvi = FindTriangleContening(vi->i,dete);2628 // Vertex * nv = quadtree->NearestVertex(vi->i.x,vi->i.y);2629 // cout << " Neart Vertex of " << Number(vi)<< vi->i << " is "2630 // << Number(nv) << nv->i << endl;2631 // Int4 kt = Number(tcvi);2632 //2633 2634 quadtree->Add(*vi); //2635 if (tcvi->det<0){// internal2636 throw ErrorException(__FUNCT__,exprintf("tcvi->det<0"));2637 }2638 Add(*vi,tcvi,dete),NbSwap += vi->Optim(1);2639 }2640 }2641 cout << " Nb swap = " << NbSwap << " after " ;2642 2643 for (i=nbvold;i<nbv;i++)2644 NbSwap += vertices[i].Optim(1);2645 cout << NbSwap << endl;2646 2647 for (i=nbtold;i<nbt;i++)2648 first_np_or_next_t[i]=1;2649 2650 Headt = nbt; // empty list2651 for (i=nbvold;i<nbv;i++)2652 { // for all the triangle contening the vertex i2653 Vertex * s = vertices + i;2654 TriangleAdjacent ta(s->t, EdgesVertexTriangle[s->vint][1]);2655 Triangle * tbegin= (Triangle*) ta;2656 Int4 kt;2657 do {2658 kt = Number((Triangle*) ta);2659 if (first_np_or_next_t[kt]>0)2660 first_np_or_next_t[kt]=-Headt,Headt=kt;2661 if (ta.EdgeVertex(0) != s){2662 throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0) != s"));2663 }2664 ta = Next(Adj(ta));2665 } while ( (tbegin != (Triangle*) ta));2666 }2667 2668 2669 } while (nbv!=nbvold);2670 delete [] first_np_or_next_t;2671 cout << " end : Triangles::NewPoints old nbv=" << nbv << endl;2672 2673 }2674 /*}}}1*/2675 2350 /*FUNCTION Triangles::Insert{{{1*/ 2676 2351 void Triangles::Insert() { 2677 2352 long int verbosity=2; 2678 if (verbosity>2) cout << " -- Insert initial " << nbv << " vertices " << endl;2353 if (verbosity>2) printf(" Insert initial %i vertices\n",nbv); 2679 2354 Triangles * OldCurrentTh =CurrentTh; 2680 2355 2681 2356 CurrentTh=this; 2682 double time0=CPUtime(),time1,time2,time3;2683 2357 SetIntCoor(); 2684 2358 Int4 i; … … 2737 2411 quadtree->Add(*v1); 2738 2412 2739 // on ajoute les sommets un Ò un2413 //the vertices are added one by one 2740 2414 Int4 NbSwap=0; 2741 2415 2742 time1=CPUtime(); 2743 2744 if (verbosity>3) cout << " -- Begin of insertion process " << endl; 2416 2417 if (verbosity>3) printf(" Begining of insertion process...\n"); 2745 2418 2746 2419 for (Int4 icount=2; icount<nbv; icount++) { 2747 2420 Vertex *vi = ordre[icount]; 2748 // cout << " Insert " << Number(vi) << endl;2749 2421 Icoor2 dete[3]; 2750 2422 Triangle *tcvi = FindTriangleContening(vi->i,dete); … … 2752 2424 Add(*vi,tcvi,dete); 2753 2425 NbSwap += vi->Optim(1,0); 2754 2755 }// fin de boucle en icount 2756 time2=CPUtime(); 2757 if (verbosity>3) 2758 cout << " NbSwap of insertion " << NbSwap 2759 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2760 << " NbUnSwap " << NbUnSwap << " Nb UnSwap/Nbv " 2761 << (float)NbUnSwap /(float) nbv 2762 <<endl; 2426 } 2427 if (verbosity>3) { 2428 printf(" NbSwap of insertion: %i\n",NbSwap); 2429 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 2430 printf(" NbUnSwap: %i\n",NbUnSwap); 2431 printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); 2432 } 2763 2433 NbUnSwap = 0; 2764 2434 // construction d'un ordre aleatoire … … 2770 2440 ordre[is4]= &vertices[k3 = (k3 + PrimeNumber)% nbv]; 2771 2441 2772 double timeloop = time2 ;2773 2442 for(int Nbloop=0;Nbloop<NBLOOPOPTIM;Nbloop++) 2774 2443 { 2775 double time000 = timeloop;2776 2444 Int4 NbSwap = 0; 2777 2445 for (int is1=0; is1<nbv; is1++) 2778 2446 NbSwap += ordre[is1]->Optim(0,0); 2779 timeloop = CPUtime(); 2780 if (verbosity>3) 2781 cout << " Optim Loop "<<Nbloop<<" NbSwap: " << NbSwap 2782 << " NbSwap/Nbv " << (float) NbSwap / (float) nbv 2783 << " CPU=" << timeloop - time000 << " s, " 2784 << " NbUnSwap/Nbv " << (float)NbUnSwap /(float) nbv 2785 << endl; 2447 if (verbosity>3) { 2448 printf(" Optim Loop: %i\n",Nbloop); 2449 printf(" NbSwap/nbv: %i\n",NbSwap/nbv); 2450 printf(" NbUnSwap: %i\n",NbUnSwap); 2451 printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); 2452 } 2786 2453 NbUnSwap = 0; 2787 2454 if(!NbSwap) break; … … 2790 2457 // because we break the TriangleContainingTheVertex 2791 2458 #endif 2792 time3=CPUtime();2793 if (verbosity>4)2794 cout << " init " << time1 - time0 << " initialisation, "2795 << time2 - time1 << "s, insert point "2796 << time3 -time2 << "s, optim " << endl2797 << " Init Total Cpu Time = " << time3 - time0 << "s " << endl;2798 2799 2459 CurrentTh=OldCurrentTh; 2800 2460 } … … 2803 2463 void Triangles::ForceBoundary() { 2804 2464 long int verbosity=2; 2805 if (verbosity > 2) 2806 cout << " -- ForceBoundary nb of edge " << nbe << endl; 2465 if (verbosity > 2) printf(" ForceBoundary nb of edge: %i\n",nbe); 2807 2466 int k=0; 2808 2467 Int4 nbfe=0,nbswp=0,Nbswap=0; 2809 for (Int4 t = 0; t < nbt; t++) 2810 if (!triangles[t].det)2811 k++,cerr << " det T" << t << " = " << 0 << endl;2468 for (Int4 t = 0; t < nbt; t++){ 2469 if (!triangles[t].det) k++; 2470 } 2812 2471 if (k!=0) { 2813 2472 throw ErrorException(__FUNCT__,exprintf("there is %i triangles of mes = 0",k)); … … 2822 2481 else Nbswap += nbswp; 2823 2482 if (nbswp) nbfe++; 2824 if ( nbswp < 0 && k < 5) 2825 { 2826 cerr << " Missing Edge " << i << " v0 = " << Number(edges[i][0]) << edges[i][0].r 2827 <<" v1= " << Number(edges[i][1]) << edges[i][1].r << " " << edges[i].on->Cracked() << " " << (Triangle *) ta ; 2828 if(ta.t) 2829 { 2830 Vertex *aa = ta.EdgeVertex(0), *bb = ta.EdgeVertex(1); 2831 cerr << " crossing with [" << Number(aa) << ", " << Number(bb) << "]\n"; 2832 } 2833 else cerr << endl; 2834 2483 if ( nbswp < 0 && k < 5){ 2484 throw ErrorException(__FUNCT__,exprintf("Missing Edge %i, v0=%i,v1=%i",i ,Number(edges[i][0]),Number(edges[i][1]))); 2835 2485 } 2836 2486 if ( nbswp >=0 && edges[i].on->Cracked()) … … 2844 2494 for (Int4 j=0;j<nbv;j++) 2845 2495 Nbswap += vertices[j].Optim(1,0); 2846 if (verbosity > 3) 2847 cout << " Nb of inforced edge = " << nbfe << " Nb of Swap " << Nbswap << endl; 2848 2496 if (verbosity > 3) printf(" number of inforced edge = %i, number of swap= %i\n",nbfe,Nbswap); 2849 2497 } 2850 2498 /*}}}1*/ … … 2853 2501 long int verbosity=0; 2854 2502 2855 if (verbosity >2) 2856 { 2857 if (OutSide) 2858 cout << " -- Find all external sub-domain "; 2859 else 2860 cout << " -- Find all internal sub-domain "; 2861 if(verbosity>99) 2862 { 2863 2864 for(int i=0;i<nbt;++i) 2865 cout << i<< " " << triangles[i] << endl; 2866 } 2867 2503 if (verbosity >2){ 2504 if (OutSide) printf(" Find all external sub-domain\n"); 2505 else printf(" Find all internal sub-domain\n"); 2868 2506 } 2869 // if (verbosity > 4) cout << " OutSide=" << OutSide << endl;2870 2507 short * HeapArete = new short[nbt]; 2871 2508 Triangle ** HeapTriangle = new Triangle* [nbt]; … … 2919 2556 // infini triangle 2920 2557 NbSubDomTot --; 2921 // cout << " triangle infini " << it << triangles[it] << endl;2922 2558 t=&triangles[it]; 2923 2559 NbOutT--; // on fait un coup de trop. 2924 while (t){ // cout << Number(t) << " " << endl;2560 while (t){ 2925 2561 NbOutT++; 2926 2562 t1=t; 2927 2563 t=t->link; 2928 t1->link=0;}//while (t) 2564 t1->link=0; 2565 } 2929 2566 } 2930 2567 } … … 2960 2597 mark[it]=k; 2961 2598 subdomains[k].head = t1; 2962 // cout << " New -- " << Number(t1) << " " << it << endl; 2963 do {// cout << " k " << k << " " << Number(t) << endl; 2599 do { 2964 2600 mark[Number(t)]=k; 2965 2601 t=t->link; … … 2985 2621 Int4 kr = mark[it]; 2986 2622 if(kr !=kl) { 2987 //cout << kl << " " << kr << " rl " << subdomains[kl].ref2988 // << " rr " << subdomains[kr].ref ;2989 2623 if (kl >=0 && subdomains[kl].ref <0 && kr >=0 && subdomains[kr].ref>=0) 2990 2624 nbk--,subdomains[kr].ref=subdomains[kl].ref-1; … … 2995 2629 if(kl<0 && kr >=0 && subdomains[kr].ref>=0) 2996 2630 nbk--,subdomains[kr].ref=-1; 2997 // cout << " after \t "2998 // << kl << subdomains[kl].ref << " rr " << kr2999 // << subdomains[kr].ref << endl;3000 2631 } 3001 2632 } 3002 // cout << subdomains[0].ref << subdomains[1].ref << endl;3003 2633 Int4 j=0; 3004 2634 for ( i=0;i<NbSubDomains;i++) 3005 2635 if((-subdomains[i].ref) %2) { // good 3006 //cout << " sudom ok = " << i << " " << subdomains[i].ref3007 // << " " << (-subdomains[i].ref) %2 << endl;3008 2636 if(i != j) 3009 2637 Exchange(subdomains[i],subdomains[j]); 3010 2638 j++;} 3011 else 3012 { //cout << " remove sub domain " << i << endl; 2639 else{ 3013 2640 t= subdomains[i].head; 3014 while (t){// cout << Number(t) << " " << endl;2641 while (t){ 3015 2642 NbOutT++; 3016 2643 t1=t; 3017 2644 t=t->link; 3018 t1->link=0;}//while (t) 2645 t1->link=0; 2646 }//while (t) 3019 2647 } 3020 3021 if(verbosity>4) 3022 cout << " Number of remove sub domain (OutSideMesh) =" << NbSubDomains-j << endl; 2648 if(verbosity>4) printf(" Number of removes subdomains (OutSideMesh) = %i\n",NbSubDomains-j); 3023 2649 NbSubDomains=j; 3024 2650 } … … 3034 2660 subdomains = new SubDomain[ Gh.NbSubDomains]; 3035 2661 NbSubDomains =Gh.NbSubDomains; 3036 if(verbosity>4)3037 cout << " find the " << NbSubDomains << " sub domain " << endl;3038 2662 Int4 err=0; 3039 2663 ReMakeTriangleContainingTheVertex(); … … 3058 2682 int sens = Gh.subdomains[i].sens; 3059 2683 // test if ge and e is in the same sens 3060 // cout << " geom edge = " << Gh.Number(eg) <<" @" << &eg << " ref = " << subdomains[i].ref3061 // << " ref edge =" << eg.ref << " sens " << sens ;3062 2684 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) sens = -sens ; 3063 2685 subdomains[i].sens = sens; … … 3073 2695 throw ErrorException(__FUNCT__,exprintf("v0 != ta.EdgeVertex(1)")); 3074 2696 } 3075 // cout << " recherche " << Number( ta.EdgeVertex(0)) << endl;3076 2697 if (ta.EdgeVertex(0) == v1) { // ok we find the edge 3077 2698 if (sens>0) … … 3079 2700 else 3080 2701 subdomains[i].head=t=ta; 3081 //cout << " triangle =" << Number(t) << " = " << (*t)[0].r << (*t)[1].r << (*t)[2].r << endl;3082 2702 if(t<triangles || t >= triangles+nbt || t->det < 0 || t->link == 0) { 3083 2703 throw ErrorException(__FUNCT__,exprintf("bad definition of SubSomain %i",i)); … … 3085 2705 Int4 it = Number(t); 3086 2706 if (mark[it] >=0) { 3087 if(verbosity>10){3088 cerr << " Warning: the sub domain " << i << " ref = " << subdomains[i].ref3089 << " is previouly defined with " <<mark[it] << " ref = " << subdomains[mark[it]].ref3090 << " skip this def " << endl;3091 }3092 2707 break; 3093 2708 } … … 3106 2721 tt=tt->link; 3107 2722 } while (tt!=t); 3108 if(verbosity>7)3109 cout << " Nb de triangles dans le sous domaine " << i << " de ref " << subdomains[i].ref << " = " << kkk << endl;3110 2723 break; 3111 2724 } … … 3118 2731 3119 2732 if (inew < NbSubDomains) { 3120 if (verbosity>5) 3121 cout << " Warning: We remove " << NbSubDomains-inew << " SubDomains " << endl; 2733 if (verbosity>5) printf("WARNING: %i SubDomains are being removed\n",NbSubDomains-inew); 3122 2734 NbSubDomains=inew;} 3123 2735 … … 3133 2745 for (it=0;it<nbt;it++) 3134 2746 if(!triangles[it].link) NbOutT++; 3135 if (verbosity> 4)3136 cout << " " ;3137 if (verbosity> 2)3138 cout << " Nb of Sub borned Domain = " << NbSubDomTot << " NbOutTriangles = " << NbOutT <<endl;3139 3140 3141 2747 } 3142 2748 /*}}}1*/ … … 3236 2842 while (t0 != (t=t->link)); 3237 2843 } 3238 if (verbosity>9)3239 cout << " number of inside triangles " << k << " nbt = " << nbt << endl;3240 2844 // take is same numbering if possible 3241 2845 if(justcompress) … … 3300 2904 while (t0 != (t=t->link)); 3301 2905 } 3302 // NbOutT = nbt - k;3303 if (verbosity>5)3304 cout << " Nb of Sub Domain =" << NbSubDomains << " Nb of In Triangles " << k3305 << " Nbt = " << nbt << " Out Triangles = " << nbt - k << endl;3306 3307 2906 return k; 3308 3309 2907 } 3310 2908 /*}}}1*/ … … 3375 2973 Gh[i].to = vertices + nbv;// save Geom -> Th 3376 2974 VerticesOnGeomVertex[nbv]= VertexOnGeom(vertices[nbv],Gh[i]); 3377 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;3378 2975 nbv++; 3379 2976 } … … 3400 2997 Gh.UnMarkEdges(); 3401 2998 int bfind=0; 3402 /*3403 cout << " nb curves = " << Gh.NbOfCurves << endl;3404 for(int i=0;i<Gh.NbOfCurves ;i++)3405 {3406 cout << " Curve " << i << " begin e=" << Gh.Number(Gh.curves[i].be) << " k=" << Gh.curves[i].kb3407 << " end e= " << Gh.Number(Gh.curves[i].ee) << " k=" << Gh.curves[i].ke << endl;3408 }*/3409 2999 for (int i=0;i<Gh.NbOfCurves;i++) 3410 3000 { … … 3420 3010 // a begin of curve 3421 3011 int nc = ei.on->CurveNumber; 3422 3423 //cout << "curve " << nc << " v " << Gh.Number((GeometricalVertex *) *ei[je].on) << " "3424 // << " e " << " " << Gh.Number(ei.on) << " vc " << Gh.Number((*Gh.curves[nc].be)[Gh.curves[nc].kb]) << endl;3425 3012 if( 3426 3013 ei.on==Gh.curves[nc].be && … … 3428 3015 ) 3429 3016 { 3430 // cout << " find " << endl;3431 3017 bcurve[nc]=iedge*2+je; 3432 3018 bfind++; … … 3449 3035 Int4 iedge; 3450 3036 Gh.UnMarkEdges(); 3451 /* add Curve loop FH3452 // find a starting back groud edges to walk3453 for (iedge=0;iedge<BTh.nbe;iedge++) {3454 Edge & ei = BTh.edges[iedge];3455 for(int jedge=0;jedge<2;jedge++) // for the 2 extremites3456 if (!ei.on->Mark() && ei[jedge].on->IsRequiredVertex() )3457 {3458 */3459 // new code FH 20043460 3037 Real8 L=0; 3461 3038 for (int icurve=0;icurve<Gh.NbOfCurves;icurve++) … … 3472 3049 Int4 NbCreatePointOnCurve=0;// Nb of new points on curve (phase==1) 3473 3050 3474 3475 // cout.precision(16);3476 3051 for(int phase=0;phase<=step;phase++) 3477 3052 { … … 3564 3139 A1->ReferenceNumber = eeequi.ref; 3565 3140 A1->DirOfSearch =NoDirOfSearch; 3566 //cout << icurveequi << " " << i << " " << *A1 << endl;3567 3141 e->on = ongequi; 3568 3142 e->v[0]= A0; 3569 3143 e->v[1]= A1; 3570 if(verbosity>99)3571 cout << i << "+ New P "<< nbv-1 << " " <<sNew<< " L0=" << L03572 << " AB=" << LAB << " s=" << (sNew-L0)/LAB << " se= "3573 << se <<" B edge " << BTh.Number(ee) << " signe = " << k1 <<" " << A1->r <<endl;3574 3575 3144 e->ref = eeequi.ref; 3576 3145 e->adj[0]=PreviousNewEdge; … … 3581 3150 A0=A1; 3582 3151 sNew += Lstep; 3583 // cout << " sNew = " << sNew << " L = " << L3584 // << " ------" <<NbCreatePointOnCurve << " " << i << endl;3585 3152 if (++i== NbCreatePointOnCurve) break; 3586 3153 } … … 3590 3157 throw ErrorException(__FUNCT__,exprintf("ee.on->CurveNumber!=ei.on->CurveNumber")); 3591 3158 } 3592 if(verbosity>98) cout << BTh.Number(ee) << " " << " on=" << *ee[k1].on << " "<< ee[k1].on->IsRequiredVertex() << endl;3593 3159 if ( ee[k1].on->IsRequiredVertex()) { 3594 3160 if (!eeequi[k1equi].on->IsRequiredVertex()){ … … 3614 3180 { 3615 3181 Edge *e = edges + nbe++; 3616 if (verbosity>10)3617 cout << " Fin curve A1" << *A1 << " " << icurve << " <=> " << icurveequi <<"-----" <<3618 NbCreatePointOnCurve << " == " <<i<<endl;3619 3182 e->on = ongequi; 3620 3183 e->v[0]= A0; … … 3645 3208 NbOfNewPoints += NbCreatePointOnCurve; 3646 3209 } 3647 if(verbosity>5)3648 cout << icurve << " NbSegOnCurve = " << NbSegOnCurve << " Lstep="3649 << Lstep <<" " << NbOfNewPoints<< " NBPC= " << NbCreatePointOnCurve <<endl;3650 3210 // do'nt 3651 3211 // if(NbCreatePointOnCurve<1) break; … … 3719 3279 Gh[i].to = vertices + nbv;// save Geom -> Th 3720 3280 VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); 3721 // cout << "--------- " <<Number(Gh[i].to) << " " << Gh[i].to << " " << i << endl;3722 3281 nbv++; 3723 3282 } … … 3738 3297 nbe = 0; 3739 3298 Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; 3740 // cout << " -------------- step =" << step << endl;3741 3299 Gh.UnMarkEdges(); 3742 3300 NbOfCurves = 0; … … 3747 3305 if (!ei.Mark() && ei[j].Required()) { 3748 3306 // warning ei.Mark() can be change in loop for(j=0;j<2;j++) 3749 // cout << " New curve = " << NbOfCurves << endl;3750 3307 Int4 nbvend =0; 3751 3308 … … 3794 3351 va = a->to; 3795 3352 e->SetMark(); 3796 // cout << " New curve " ;3797 3353 3798 3354 // if SameGeo We have go in the background geometry … … 3818 3374 A= a->r; 3819 3375 Metric MAs =MA,MBs; 3820 // cout << " lSubEdge old=" << ledge3821 // << " new " << A << MA << endl;3822 3376 ledge = 0; 3823 3377 Real8 x =0, xstep= 1. / NbSubEdge; … … 3828 3382 AB = A-B; 3829 3383 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); 3830 // cout << " " << lSubEdge[kk] << " x " << x3831 // << " " << A << B << MA << MB<< endl ;3832 3384 } 3833 // cout << endl;3834 3385 } 3835 3386 … … 3884 3435 lcurve = lcurveb; 3885 3436 e->SetMark(); 3886 // cout << e-Gh.edges << ", " << k << " "3887 // <<(*e)[k].r <<" " <<(*e)[1-k].r <<" "3888 // << lcurve<< ";; " ;3889 3437 a=b; 3890 3438 if (b->Required() ) break; … … 3897 3445 }// for(;;) 3898 3446 vb = b->to; 3899 // cout << endl;3900 3447 NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1); 3901 3448 NbNewPoints = NbEdgeCurve-1; … … 3907 3454 3908 3455 lstep = lcurve / NbEdgeCurve; 3909 // cout <<"lstep " << lstep << " lcurve "3910 // << lcurve << " NbEdgeCurve " << NbEdgeCurve << " " <<NbVerticesOnGeomEdge0<<" " <<NbVerticesOnGeomEdge<<" step =" <<step<< endl;3911 3456 } 3912 3457 // end of curve -- … … 3979 3524 int j= ii; 3980 3525 while (!(*on)[j].Required()) { 3981 // cout << i << " " << ii << " j= " << j << " curve = "3982 // << on->CurveNumber << " " << Gh.Number(on) << " on " << j3983 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1])3984 // << " -> " ;3985 3526 Adj(on,j); // next geom edge 3986 3527 j=1-j; 3987 // cout << Gh.Number(on) << " " << j << " e[ON] = " << e[Gh.Number(on)]3988 // << " s0 " << Gh.Number( (*on)[0]) << " s1 " << Gh.Number( (*on)[1]) << endl;3989 3528 if (e[Gh.Number(on)]) break; // optimisation 3990 3529 e[Gh.Number(on)] = ei; … … 3993 3532 3994 3533 int kk=0; 3995 for ( i=0;i<Gh.nbe ; i++) 3996 if (!e[i])3997 if(kk++<10) {3998 cerr << " Bug -- the geometrical edge " << i << " is on no edge curve = " << Gh.edges[i].CurveNumber3999 << " s0 " << Gh.Number( Gh.edges[i][0]) << " s1 " << Gh.Number( Gh.edges[i][1]) << endl;4000 3534 for ( i=0;i<Gh.nbe ; i++){ 3535 if (!e[i]){ 3536 kk++; 3537 if(kk<10) printf("BUG: the geometrical edge %i is on no edge curve\n",i); 3538 } 3539 } 4001 3540 if(kk) throw ErrorException(__FUNCT__,exprintf("See above")); 4002 3541 … … 4217 3756 for (Int4 icount=2; icount<nbvb; icount++) { 4218 3757 Vertex *vi = ordre[icount]; 4219 // cout << " Add vertex " << Number(vi) << endl;4220 3758 Icoor2 dete[3]; 4221 3759 Triangle *tcvi = FindTriangleContening(vi->i,dete); … … 4322 3860 triangles = savetriangles; 4323 3861 subdomains = savesubdomains; 4324 // cout << triangles << " <> " << OutSidesTriangles << endl;4325 /* k=0;4326 for (i=0;i<nbt;i++)4327 for (int j=0;j<3;j++)4328 if (!triangles[i].TriangleAdj(j))4329 k++;4330 */4331 3862 if (k) { 4332 3863 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); 4333 3864 } 4334 3865 FindSubDomain(); 4335 // cout << " NbTOld = " << NbTold << " == " << nbt - NbOutT << " " << nbt << endl;4336 3866 4337 3867 delete edge4; … … 4347 3877 if (!edges[i].adj[j]) 4348 3878 if(!edges[i][j].on->IsRequiredVertex()) { 4349 cerr << " Erreur adj et sommet requis edges [" << i << "][ " << j << "]= " 4350 << Number(edges[i][j]) << " : " << " on = " << Gh.Number(edges[i].on) ; 4351 if (edges[i][j].on->OnGeomVertex()) 4352 cerr << " vertex " << Gh.Number(edges[i][j].on->gv); 4353 else if (edges[i][j].on->OnGeomEdge()) 4354 cerr << "Edges " << Gh.Number(edges[i][j].on->ge); 4355 else 4356 cerr << " = " << edges[i][j].on ; 4357 cerr << endl; 3879 throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)")); 4358 3880 } 4359 3881 } … … 4372 3894 jp = aa[jp]&3; 4373 3895 do { 4374 // cout << *t << " " << j << "\n\t try swap " ;4375 3896 while (t->swap(j,koption)) 4376 3897 { … … 4382 3903 t= tp->at[jp]; // set unchange t qnd j for previous triangles 4383 3904 j= NextEdge[tp->aa[jp]&3]; 4384 // cout << "\n\t s " << *t << " " << j << endl;4385 3905 } 4386 3906 // end on this Triangle … … 4416 3936 for ( k=0;k<NbVerticesOnGeomEdge;k++ ) 4417 3937 tstart[ Number(VerticesOnGeomEdge[k].mv)]=&vide; 4418 if(verbosity>2) 4419 cout << " -- SmoothingVertex: nb Iteration = " << nbiter << " Omega = " << omega << endl; 3938 if(verbosity>2) printf(" SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega); 4420 3939 for (k=0;k<nbiter;k++) 4421 3940 { … … 4429 3948 if (tstart[i] != &vide) // not a boundary vertex 4430 3949 NbSwap += vertices[i].Optim(1); 4431 if (verbosity>3) 4432 cout << " Move max = " << sqrt(delta) << " iteration = " 4433 << k << " Nb of Swap = " << NbSwap << endl; 3950 if (verbosity>3) printf(" move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap); 4434 3951 } 4435 3952 … … 4441 3958 void Triangles::MakeQuadTree() { 4442 3959 long int verbosity=0; 4443 if(verbosity>8)4444 cout << " MakeQuadTree" << endl;4445 3960 if ( !quadtree ) quadtree = new QuadTree(this); 4446 3961 … … 4454 3969 D2xD2 Br(D2xD2(Beq,Heq).t()); 4455 3970 D2xD2 B1r(Br.inv()); 4456 /* D2xD2 BB = Br.t()*Br;4457 cout << " BB = " << BB << " " << Br*B1r << endl;4458 MetricAnIso MMM(BB.x.x,BB.x.y,BB.y.y);4459 MatVVP2x2 VMM(MMM);4460 cout << " " << VMM.lambda1 << " " << VMM.lambda2 << endl;4461 */4462 3971 double gammamn=1e100,hmin=1e100; 4463 3972 double gammamx=0,hmax=0; … … 4484 3993 MatVVP2x2 VMK(MK); 4485 3994 alpha2 = Max(alpha2,Max(VMK.lambda1/VMK.lambda2,VMK.lambda2/VMK.lambda1)); 4486 // cout << B_K << " * " << B1r << " == " << BK << " " << B_K*B_K.inv() << endl;4487 3995 Real8 betaK=0; 4488 3996 … … 4499 4007 MetricAnIso M1(BMB.x.x,BMB.x.y,BMB.y.y); 4500 4008 MatVVP2x2 VM1(M1); 4501 //cout << B_K <<" " << M << " " << he << " " << BMB << " " << VM1.lambda1 << " " << VM1.lambda2<< endl;4502 4009 gammamn=Min3(gammamn,VM1.lambda1,VM1.lambda2); 4503 4010 gammamx=Max3(gammamx,VM1.lambda1,VM1.lambda2); … … 4505 4012 betaK *= area3;// 1/2 (somme sqrt(det))* area(K) 4506 4013 Marea+= betaK; 4507 // cout << betaK << " " << area3 << " " << beta << " " << beta0 << " " << area3*3*3*3 <<endl;4508 4014 beta=min(beta,betaK); 4509 4015 beta0=max(beta0,betaK); … … 4512 4018 gammamn=sqrt(gammamn); 4513 4019 gammamx=sqrt(gammamx); 4514 cout << " -- adaptmesh Regulary: Nb triangles " << nt << " , h min " << hmin << " , h max " << hmax << endl; 4515 cout << " area = " << area << " , M area = " << Marea << " , M area/( |Khat| nt) " << Marea/(aireKh*nt) << endl; 4516 cout << " infiny-regulaty: min " << gammamn << " max " << gammamx << endl; 4517 cout << " anisomax "<< sqrt(alpha2) << ", beta max = " << 1./sqrt(beta/aireKh) 4518 << " min "<< 1./sqrt(beta0/aireKh) << endl; 4020 printf(" Adaptmesh info:\n"); 4021 printf(" number of triangles = %i\n",nt); 4022 printf(" hmin = %g, hmax=%g\n",hmin,hmax); 4023 printf(" area = %g, M area = %g, M area/( |Khat| nt) = %g\n",area,Marea, Marea/(aireKh*nt)); 4024 printf(" infinite-regularity(?): min = %g, max = %g\n",gammamn,gammamx); 4025 printf(" anisomax = %g, beta max = %g, min = %g\n",pow(alpha2,0.5),1./pow(beta/aireKh,0.5), 1./pow(beta0/aireKh,0.5)); 4519 4026 } 4520 4027 /*}}}1*/ … … 4549 4056 } 4550 4057 } 4551 cout << " -- Histogram of the unit mesh, nb of edges" << nbedges << endl <<endl; 4552 4553 cout << " length of edge in | % of edge | Nb of edges " << endl; 4554 cout << " ------------------- | ---------- | ----------- " << endl; 4555 for (i=0;i<=kmax;i++) 4556 { 4557 cout << " " ; 4558 cout.width(10); 4559 if (i==0) cout << " 0 " ; 4560 else cout << exp(lmin+i/delta) ; 4561 cout.width(); cout << "," ; 4562 cout.width(10); 4563 if (i==kmax) cout << " +infty " ; 4564 else cout << exp(lmin+(i+1)/delta) ; 4565 cout.width();cout << " | " ; 4566 4567 cout.precision(4); 4568 cout.width(6); 4569 cout << ((long) ((10000.0 * histo[i])/ nbedges))/100.0 ; 4570 cout.width(); 4571 cout.precision(); 4572 cout << " | " << histo[i] <<endl; 4058 printf(" --- Histogram of the unit mesh, nb of edges = %i\n",nbedges); 4059 printf(" length of edge in | %% of edge | Nb of edges \n"); 4060 printf(" --------------------+-------------+-------------\n"); 4061 for (i=0;i<=kmax;i++){ 4062 if (i==0) printf(" %10i",0); 4063 else printf(" %10g",exp(lmin+i/delta)); 4064 if (i==kmax) printf(" +inf "); 4065 else printf(" %10g",exp(lmin+(i+1)/delta)); 4066 printf("| %10g |\n",((long) ((10000.0 * histo[i])/ nbedges))/100.0); 4067 printf(" %g\n",histo[i]); 4573 4068 } 4574 cout << " ------------------- | ---------- | ----------- " << endl <<endl; 4575 4069 printf(" --------------------+-------------+-------------\n"); 4576 4070 } 4577 4071 /*}}}1*/ … … 4605 4099 if( k==0) return 0; 4606 4100 CurrentTh = this; 4607 cout << " Nb of Cracked Edges = " << k << endl;4101 printf(" number of Cracked Edges = %i\n",k); 4608 4102 NbCrackedEdges =k; 4609 4103 CrackedEdges = new CrackedEdge[k]; … … 4689 4183 4690 4184 // set the ref 4691 cout << " set the ref " << endl ;4692 4185 NbCrackedVertices = nbcrakev; 4693 4186 // int nbvo = nbv; 4694 4187 nbv = vlast - vertices; 4695 4188 int nbnewv = nbv - nbv; // nb of new vrtices 4696 if (nbcrakev && verbosity > 1 ) 4697 cout << " Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nbnewv<< endl; 4189 if (nbcrakev && verbosity > 1 ) printf(" number of Cracked vertices = %i, number of created vertices = %i\n",nbcrakev,nbnewv); 4698 4190 // all the new vertices are on geometry 4699 // BOFBO-- A VOIR4700 4191 if (nbnewv) 4701 4192 { // … … 4738 4229 throw ErrorException(__FUNCT__,exprintf("Not enougth vertices: to crack the mesh we need %i vertices",nbv)); 4739 4230 } 4740 cout << " NbCrackedVertices " << NbCrackedVertices << endl;4741 4231 CurrentTh = CurrentThOld; 4742 4232 return NbCrackedVertices; … … 4835 4325 void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){ 4836 4326 long int verbosity=0; 4837 4838 if(verbosity>1)4839 cout << " -- IntersectGeomMetric geometric err=" << err << (iso ? " iso " : " aniso " ) << endl;4840 4327 Real8 ss[2]={0.00001,0.99999}; 4841 4328 Real8 errC = 2*sqrt(2*err); … … 4854 4341 Vertex V; 4855 4342 VertexOnGeom GV; 4856 // cerr << Number(edges[i]) << " " << ss[j] << endl;4857 4343 Gh.ProjectOnCurve(edges[i],ss[j],V,GV); 4858 4344 { … … 4870 4356 } 4871 4357 MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); 4872 //cerr << " : " ;4873 4358 Metric MVp(Vp); 4874 // cerr << " : " << MVp << endl;4875 4359 edges[i][j].m.IntersectWith(MVp); 4876 //cerr << " . " << endl;4877 4360 } 4878 4361 … … 4886 4369 4887 4370 double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100)); 4888 if (verbosity > 1) 4889 cout << " -- BoundAnisotropy by " << anisomax << endl; 4371 if (verbosity > 1) printf(" BoundAnisotropy by %g\n",anisomax); 4890 4372 Real8 h1=1.e30,h2=1e-30,rx=0; 4891 4373 Real8 coef = 1./(anisomax*anisomax); … … 4913 4395 } 4914 4396 4915 if (verbosity>2) 4916 { 4917 cout << " input : Hmin = " << sqrt(1/h2) << " Hmax = " << sqrt(1/h1) 4918 << " factor of anisotropy max = " << sqrt(rx) << endl; 4919 cout << " output: Hmin = " << sqrt(1/hn2) << " Hmax = " << sqrt(1/hn1) 4920 << " factor of anisotropy max = " << sqrt(rnx) << endl; 4397 if (verbosity>2){ 4398 printf(" input: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5)); 4399 printf(" output: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 4921 4400 } 4922 4401 } … … 4960 4439 Real8 coef2 = 1/(coef*coef); 4961 4440 4962 if(verbosity>1) 4963 { 4964 cout << " -- Construction of Metric: Nb of field. " << n << " nbt = " 4965 << nbt << " nbv= " << nbv 4966 << " coef = " << coef << endl 4967 << " hmin = " << hmin << " hmax=" << hmax 4968 << " anisomax = " << anisomax << " Nb Jacobi " << NbJacobi << " Power = " << power ; 4969 if (RelativeMetric) 4970 cout << " RelativeErr with CutOff= " << CutOff << endl; 4971 else 4972 cout << " Absolute Err" <<endl; 4441 if(verbosity>1) { 4442 printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); 4443 printf(" coef = %g\n",coef); 4444 printf(" hmin = %g hmax = %g\n",hmin,hmax); 4445 printf(" anisomax = %g nb Jacobi = %i, power = %i\n",anisomax,NbJacobi,power); 4446 if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); 4447 else printf(" Absolute error\n"); 4973 4448 } 4974 double *ss=(double*)s; //, *ssiii = ss;4449 double *ss=(double*)s; 4975 4450 4976 4451 double sA,sB,sC; … … 5070 4545 Real8 cnorm = DoNormalisation ? coef2/sdelta : coef2; 5071 4546 5072 if(verbosity>2) 5073 cout << " Solution " << nusol << " Min = " << smin << " Max = " 5074 << smax << " Delta =" << sdelta << " cnorm = " << cnorm << " Nb of fields =" << nbfield << endl; 5075 5076 5077 if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) 5078 { 5079 if (verbosity>2) 5080 cout << " Solution " << nusol << " is constant. We skip. " 5081 << " Min = " << smin << " Max = " << smax << endl; 4547 if(verbosity>2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield); 4548 4549 if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) { 4550 if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); 5082 4551 continue; 5083 4552 } … … 5144 4613 taa[1][j] = lC*lA; 5145 4614 taa[2][j] = lA*lB; 5146 //Real8 xx = V.x-V.y;5147 //Real8 yy = V.x + V.y;5148 //cout << " iv " << ss[iV*n] << " == " << (8*xx*xx+yy*yy)5149 // << " l = " << lA << " " << lB << " " << lC5150 // << " = " << lA+lB+lC << " " << V << " == " << A*lA+B*lB+C*lC << endl;5151 5152 4615 lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; 5153 4616 … … 5184 4647 Real8 coef = 1.0/(3*dd*det33); 5185 4648 Real8 coef2 = 2*coef; 5186 // cout << " H = " << Hxx << " " << Hyy << " " << Hxy/2 << " coef2 = " << coef2 << endl;5187 4649 Hxx *= coef2; 5188 4650 Hyy *= coef2; 5189 4651 Hxy *= coef2; 5190 //cout << i << " H = " << 3*Hxx/dd << " " << 3*Hyy/dd << " " << 3*Hxy/(dd*2) << " nbb = " << nbb << endl;5191 4652 if(nbb==0) 5192 4653 { … … 5245 4706 Metric M(dxdx[iv], dxdy[iv], dydy[iv]); 5246 4707 MatVVP2x2 Vp(M); 5247 //cout <<iv << " M = " << M << " aniso= " << Vp.Aniso() ;5248 4708 Vp.Abs(); 5249 4709 M = Vp; … … 5251 4711 dxdy[iv] = M.a21; 5252 4712 dydy[iv] = M.a22; 5253 // cout << " (abs) iv M = " << M << " aniso= " << Vp.Aniso() <<endl;5254 4713 } 5255 4714 else kk++; … … 5311 4770 for ( iv=0,k=0 ; iv<nbv; iv++,k+=n ) 5312 4771 { // for all vertices 5313 //{5314 //Metric M(dxdx[iv], dxdy[iv], dydy[iv]);5315 // MatVVP2x2 Vp(M);5316 //cout << " iv M="<< M << " Vp = " << Vp << " aniso " << Vp.Aniso() << endl;5317 //}5318 4772 MetricIso Miso; 5319 // new code to compute ci ---5320 4773 Real8 ci ; 5321 4774 if (RelativeMetric) … … 5329 4782 else ci = cnorm; 5330 4783 5331 // old5332 // Real8 ci = RelativeMetric ? coef2/(Max(Abs(ss[k]),rCutOff)) : cnorm ;5333 // modif F Hecht 1010995334 4784 Metric Miv(dxdx[iv]*ci, dxdy[iv]*ci, dydy[iv]*ci); 5335 4785 MatVVP2x2 Vp(Miv); … … 5358 4808 vertices[iv].m.IntersectWith(MVp); 5359 4809 }// for all vertices 5360 if (verbosity>2) 5361 { 5362 cout << " Field " << nufield << " of solution " << nusol << endl; 5363 cout << " before bounding : Hmin = " << sqrt(1/h2) << " Hmax = " 5364 << sqrt(1/h1) << " factor of anisotropy max = " << sqrt(rx) << endl; 5365 cout << " after bounding : Hmin = " << sqrt(1/hn2) << " Hmax = " 5366 << sqrt(1/hn1) << " factor of anisotropy max = " << sqrt(rnx) << endl; 4810 if (verbosity>2) { 4811 printf(" Field %i of solution %i\n",nufield,nusol); 4812 printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); 4813 printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); 5367 4814 } 5368 4815 } // end of for all field … … 5386 4833 5387 4834 const Real8 maxsubdiv2 = maxsubdiv*maxsubdiv; 5388 if(verbosity>1) 5389 cout << " -- Limit the subdivision of a edges in the new mesh by " << maxsubdiv << endl ; 4835 if(verbosity>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv); 5390 4836 // for all the edges 5391 4837 // if the len of the edge is to long … … 5434 4880 } 5435 4881 } 5436 if(verbosity>3) 5437 cout << " Nb of metric change = " << nbchange 5438 << " Max of the subdivision of a edges before change = " << sqrt(lmax) << endl; 5439 4882 if(verbosity>3){ 4883 printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5)); 4884 } 5440 4885 } 5441 4886 /*}}}1*/ … … 5445 4890 5446 4891 if(raisonmax<1.1) return; 5447 if(verbosity > 1) 5448 cout << " -- Triangles::SmoothMetric raisonmax = " << raisonmax << " " <<nbv <<endl; 4892 if(verbosity > 1) printf(" Triangles::SmoothMetric raisonmax = %g\n",raisonmax); 5449 4893 ReMakeTriangleContainingTheVertex(); 5450 4894 Int4 i,j,kch,kk,ip; … … 5464 4908 for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) { 5465 4909 // pour tous les triangles autour du sommet s 5466 // cout << kk << " i = " << i << " " << ip << endl;5467 4910 register Triangle* t= vertices[i].t; 5468 4911 if (!t){ … … 5491 4934 { 5492 4935 Real8 dh=(hj-hi)/ll; 5493 //cout << " dh = " << dh << endl;5494 4936 if (dh>logseuil) { 5495 4937 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); … … 5502 4944 { 5503 4945 Real8 li = vi.m(Aij); 5504 //Real8 lj = vj.m(Aij);5505 // if ( i == 2 || j == 2)5506 // cout << " inter " << i << " " << j << " " << ((1 +logseuil*li)) << endl;5507 4946 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) 5508 //if( vj.m.IntersectWith(vi.m*(lj/li/(1 +logseuil*lj))) )5509 4947 if(first_np_or_next_t1[j]<0) // if the metrix change 5510 4948 kch++,first_np_or_next_t1[j]=Head1,Head1=j; … … 5517 4955 Head1 = -1; 5518 4956 Exchange(first_np_or_next_t0,first_np_or_next_t1); 5519 if(verbosity>5) 5520 cout << " Iteration " << kk << " Nb de vertices with change " << kch << endl; 5521 } 5522 if(verbosity>2 && verbosity < 5) 5523 cout << " Nb of Loop " << kch << endl; 4957 } 4958 if(verbosity>2) printf(" number of iterations = %i\n",kch); 5524 4959 delete [] first_np_or_next_t0; 5525 4960 delete [] first_np_or_next_t1; … … 5599 5034 strcat(strcat(identity,";"),countbuf); 5600 5035 strcat(identity,buf); 5601 // cout << "New MAILLAGE "<< identity << endl;5602 5036 } 5603 5037 5604 5038 quadtree=0; 5605 // edgescomponante=0;5606 5039 edges=0; 5607 5040 VerticesOnGeomVertex=0; … … 5609 5042 NbVerticesOnGeomVertex=0; 5610 5043 NbVerticesOnGeomEdge=0; 5611 // nbMaxIntersectionTriangles=0;5612 // lIntTria;5613 5044 subdomains=0; 5614 5045 NbSubDomains=0; 5615 // Meshbegin = vertices;5616 // Meshend = vertices + nbvx;5617 if (verbosity>98)5618 cout << "Triangles::PreInit() " << nbvx << " " << nbtx5619 << " " << vertices5620 << " " << ordre << " " << triangles <<endl;5621 5046 } 5622 5047 /*}}}1*/ … … 5724 5149 I2 I=vI, J=vJ, IJ= J-I; 5725 5150 IJ_IA = (IJ ,(A-I)); 5726 // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <<endl;5727 5151 if (IJ_IA<0) { 5728 5152 if (dir>0) {a=1;b=0;return edge;}// change of signe => I … … 5740 5164 a= IJ_AJ/IJ2; 5741 5165 b= IJ_IA/IJ2; 5742 // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl;5743 5166 return edge; 5744 5167 } … … 5786 5209 if ( (d2=(double) ACAC) < dd2) 5787 5210 { 5788 // cout << " A " << d2 << " " << dd2;5789 5211 er = ta; 5790 5212 l0 = ACAC; … … 5798 5220 if ( (d2=(double) BCBC) < dd2) 5799 5221 { 5800 // cout << " B " << d2 << " " << dd2;5801 5222 dd2 = d2; 5802 5223 er = Adj(ta); // other direction … … 5822 5243 s = 0; 5823 5244 cas = -1; 5824 // cout << " ABAC " << ABAC << " ABAC " << ABAC5825 // << " AB2 " << AB2 << endl;5826 5245 b = ((double) ABAC/(double) AB2); 5827 5246 a = 1 - b; … … 5911 5330 if ( o < oo) 5912 5331 o=oo,pi=BigPrimeNumber[i];} 5913 // cout << " AGoodNumberPrimeWith " << n << " " <<pi << " "<< o << endl;5914 5332 return pi; 5915 5333 } … … 6050 5468 } 6051 5469 det2 = det(*v2,a,b); 6052 // cout << " No Change try the next" << endl;6053 5470 } 6054 5471 … … 6068 5485 } 6069 5486 Icoor2 detss = 0,l=0,ks; 6070 // cout << "Real ForcingEdge " << *va << *vb << detss << endl;6071 5487 while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 6072 5488 if(l++ > 10000000) {
Note:
See TracChangeset
for help on using the changeset viewer.
