![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3279
- Timestamp:
- 03/15/10 16:08:59 (15 years ago)
- Location:
- issm/trunk/src
- Files:
-
- 5 edited
-
c/Bamgx/Bamgx.cpp (modified) (3 diffs)
-
c/Bamgx/objects/Geometry.cpp (modified) (9 diffs)
-
c/Bamgx/objects/Triangles.cpp (modified) (14 diffs)
-
c/objects/BamgOpts.cpp (modified) (1 diff)
-
mex/Bamg/Bamg.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Bamgx.cpp
r3278 r3279 53 53 hmax=bamgopts->hmax; 54 54 verbosity=bamgopts->verbose; 55 56 // some verification57 if ( maxsubdiv > 10 || maxsubdiv <= 1.0){58 throw ErrorException(__FUNCT__,exprintf("maxsubdiv (%g) should be in ]1,10]",maxsubdiv));59 }60 55 61 56 // no metric -> no smoothing … … 115 110 if (verbosity>1) printf(" Processing initial mesh and geometry...\n"); 116 111 Triangles BTh(bamggeom_in,bamgmesh_in,bamgopts); 117 hmin =Max(hmin,BTh.MinimalHmin());118 hmax =Min(hmax,BTh.MaximalHmax());112 hmin=Max(hmin,BTh.MinimalHmin()); 113 hmax=Min(hmax,BTh.MaximalHmax()); 119 114 120 115 //Make Quadtree from background mesh … … 138 133 } 139 134 } 135 140 136 } 141 137 -
issm/trunk/src/c/Bamgx/objects/Geometry.cpp
r3278 r3279 3 3 #include <cmath> 4 4 #include <ctime> 5 #include <assert.h> 5 6 6 7 #include "BamgObjects.h" … … 227 228 228 229 double Hmin = HUGE_VAL;// the infinie value 229 long hvertices=0;230 230 int i,j,k,n,i1,i2; 231 231 … … 290 290 edges = new GeometricalEdge[nbe]; 291 291 292 //if hvertices==0, initialize len (length of each edge) 293 if (!hvertices) { 294 len = new double[nbv]; 295 for(i=0;i<nbv;i++) len[i]=0; 296 } 292 //initialize len (length of each edge) 293 len = new double[nbv]; 294 for(i=0;i<nbv;i++) len[i]=0; 297 295 298 296 for (i=0;i<nbe;i++){ … … 316 314 edges[i].Adj[0] = edges[i].Adj[1] = 0; 317 315 edges[i].flag = 0; 318 if (!hvertices) { 319 vertices[i1].color++;320 vertices[i2].color++;321 len[i1] += l12;322 len[i2] += l12;323 }316 317 //prepare metric 318 vertices[i1].color++; 319 vertices[i2].color++; 320 len[i1] += l12; 321 len[i2] += l12; 324 322 } 325 323 326 324 // definition the default of the given mesh size 327 if (!hvertices){ 328 for (i=0;i<nbv;i++) 329 if (vertices[i].color > 0) 330 vertices[i].m=Metric(len[i] /(double) vertices[i].color); 331 else 332 vertices[i].m=Metric(Hmin); 333 delete [] len; 334 } 325 for (i=0;i<nbv;i++) 326 if (vertices[i].color > 0) 327 vertices[i].m=Metric(len[i] /(double) vertices[i].color); 328 else 329 vertices[i].m=Metric(Hmin); 330 delete [] len; 331 335 332 } 336 333 else{ … … 353 350 if(verbose>5) printf(" processing hVertices\n"); 354 351 for (i=0;i< nbv;i++){ 355 vertices[i].m=Metric((double)bamggeom->hVertices[i]); 352 if (!isnan(bamggeom->hVertices[i])){ 353 vertices[i].m=Metric((double)bamggeom->hVertices[i]); 354 } 356 355 } 357 356 } … … 363 362 if(bamggeom->MetricVertices){ 364 363 if(verbose>5) printf(" processing MetricVertices\n"); 365 hvertices=1;366 364 for (i=0;i< nbv;i++) { 367 365 vertices[i].m = Metric((double)bamggeom->MetricVertices[i*3+0],(double)bamggeom->MetricVertices[i*3+1],(double)bamggeom->MetricVertices[i*3+2]); … … 376 374 if(verbose>5) printf(" processing h1h2VpVertices\n"); 377 375 double h1,h2,v1,v2; 378 hvertices =1;379 376 for (i=0;i< nbv;i++) { 380 377 h1=(double)bamggeom->MetricVertices[i*4+0]; … … 438 435 //RequiredVertices 439 436 if(bamggeom->RequiredVertices){ 440 if(verbose>5) printf(" processing RequiredVertices ");437 if(verbose>5) printf(" processing RequiredVertices\n"); 441 438 n=bamggeom->NumRequiredVertices; 442 439 for (i=0;i<n;i++) { … … 453 450 //RequiredEdges 454 451 if(bamggeom->RequiredEdges){ 455 if(verbose>5) printf(" processing RequiredEdges ");452 if(verbose>5) printf(" processing RequiredEdges\n"); 456 453 n=bamggeom->NumRequiredEdges; 457 454 for (i=0;i<n;i++) { -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r3278 r3279 261 261 long i1,i2,i3,iref; 262 262 long i,j; 263 long hvertices =0;264 263 Metric M1(1); 265 int Version=1,dim=2;266 264 int verbose=0; 267 265 … … 311 309 long i1,i2,i3,iref; 312 310 long i,j; 313 long hvertices =0;314 311 long ifgeom=0; 315 312 Metric M1(1); 316 int Version=1,dim=2;317 313 318 314 verbose=bamgopts->verbose; … … 387 383 } 388 384 389 //hVertices390 if(bamgmesh->hVertices){391 if(verbose>5) printf(" processing hVertices\n");392 hvertices=1;393 for (i=0;i< nbv;i++){394 if (!isnan(bamgmesh->hVertices[i])){395 vertices[i].m=Metric((double)bamgmesh->hVertices[i]);396 }397 }398 }399 else{400 if(verbose>5) printf(" no hVertices found\n");401 }402 403 385 //VerticesOnGeometricEdge 404 386 if(bamgmesh->VerticesOnGeometricEdge){ … … 443 425 nbe=bamgmesh->NumEdges; 444 426 edges= new Edge[nbe]; 445 446 if (!hvertices) { 447 len= new double[nbv]; 448 for(i=0;i<nbv;i++) 449 len[i]=0; 450 } 427 //initialize length of each edge (used to provided metric) 428 len= new double[nbv]; 429 for(i=0;i<nbv;i++) len[i]=0; 451 430 452 431 for (i=0;i<nbe;i++){ 453 432 i1=(int)bamgmesh->Edges[i*3+0]-1; //-1 for C indexing 454 433 i2=(int)bamgmesh->Edges[i*3+1]-1; //-1 for C indexing 434 edges[i].ref=(long)bamgmesh->Edges[i*3+2]; 455 435 edges[i].v[0]= vertices +i1; 456 436 edges[i].v[1]= vertices +i2; … … 460 440 double l12=sqrt((x12,x12)); 461 441 462 if (!hvertices){ 463 vertices[i1].color++; 464 vertices[i2].color++; 465 len[i1]+=l12; 466 len[i2]+=l12; 467 } 442 //prepare metric 443 vertices[i1].color++; 444 vertices[i2].color++; 445 len[i1]+=l12; 446 len[i2]+=l12; 468 447 Hmin = Min(Hmin,l12); 469 448 } 470 449 471 450 // definition the default of the given mesh size 472 if (!hvertices){ 473 for (i=0;i<nbv;i++) 474 if (vertices[i].color > 0) 475 vertices[i].m=Metric(len[i] /(double) vertices[i].color); 476 else 477 vertices[i].m=Metric(Hmin); 478 delete [] len; 479 } 451 for (i=0;i<nbv;i++){ 452 if (vertices[i].color>0) 453 vertices[i].m=Metric(len[i]/(double)vertices[i].color); 454 else 455 vertices[i].m=Metric(Hmin); 456 } 457 delete [] len; 480 458 481 459 // construction of edges[].adj 482 460 for (i=0;i<nbv;i++){ 483 vertices[i].color = (vertices[i].color ==2) ? -1 :-2;461 vertices[i].color=(vertices[i].color ==2) ?-1:-2; 484 462 } 485 463 for (i=0;i<nbe;i++){ … … 537 515 else{ 538 516 if(verbose>5) printf(" no EdgesOnGeometricEdge found\n"); 517 } 518 519 //hVertices 520 if(bamgmesh->hVertices){ 521 if(verbose>5) printf(" processing hVertices\n"); 522 for (i=0;i< nbv;i++){ 523 if (!isnan(bamgmesh->hVertices[i])){ 524 vertices[i].m=Metric((double)bamgmesh->hVertices[i]); 525 } 526 } 527 } 528 else{ 529 if(verbose>5) printf(" no hVertices found\n"); 539 530 } 540 531 … … 1586 1577 1587 1578 /*Options*/ 1588 const int dim = 2;1589 1579 double* s=NULL; 1590 1580 long nbsol; … … 1790 1780 1791 1781 /*Options*/ 1792 const int dim = 2;1793 1782 double* s=NULL; 1794 1783 long nbsol; … … 2622 2611 } 2623 2612 /*}}}1*/ 2624 /*FUNCTION Triangles::ReconstructExistingMesh{{{1*/2625 void Triangles::ReconstructExistingMesh() {2626 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/2627 2628 /*This routine reconstruct an existing mesh to make it CONVEX:2629 * -all the holes are filled2630 * -concave boundaries are filled2631 * A convex mesh is required for a lot of operations. This is why every mesh2632 * goes through this process.2633 * This routine also generates mesh properties such as adjencies,...2634 */2635 2636 /*Intermediary*/2637 int verbosity=0;2638 2639 // generation of the integer coordinate2640 2641 // find extrema coordinates of vertices pmin,pmax2642 long i;2643 if(verbosity>2) printf(" Reconstruct mesh of %i vertices\n",nbv);2644 2645 //initialize ordre2646 assert(ordre);2647 for (i=0;i<nbv;i++) ordre[i]=0;2648 2649 //Initialize NbSubDomains2650 NbSubDomains =0;2651 2652 /* generation of triangles adjacency*/2653 2654 //First add existing edges2655 long kk =0;2656 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv);2657 for (i=0;i<nbe;i++){2658 kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1])));2659 }2660 if (kk != nbe){2661 throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe));2662 }2663 2664 //Add edges of all triangles in existing mesh2665 long* st = new long[nbt*3];2666 for (i=0;i<nbt*3;i++) st[i]=-1;2667 for (i=0;i<nbt;i++){2668 for (int j=0;j<3;j++){2669 2670 //Add current triangle edge to edge42671 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]));2672 2673 long invisible=triangles[i].Hidden(j);2674 2675 //If the edge has not been added to st, add it2676 if(st[k]==-1) st[k]=3*i+j;2677 2678 //If the edge already exists, add adjacency2679 else if(st[k]>=0) {2680 assert(!triangles[i].TriangleAdj(j));2681 assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3)));2682 2683 triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3));2684 if (invisible) triangles[i].SetHidden(j);2685 if (k<nbe) triangles[i].SetLocked(j);2686 2687 //Make st[k] negative so that it will throw an error message if it is found again2688 st[k]=-2-st[k];2689 }2690 2691 //An edge belongs to 2 triangles2692 else {2693 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles",2694 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])));2695 }2696 }2697 }2698 2699 //Display info if required2700 if(verbosity>5) {2701 printf(" info of Mesh:\n");2702 printf(" - number of vertices = %i \n",nbv);2703 printf(" - number of triangles = %i \n",nbt);2704 printf(" - number of given edges = %i \n",nbe);2705 printf(" - number of all edges = %i \n" ,edge4->nb());2706 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv);2707 }2708 2709 //check the consistency of edge[].adj and the geometrical required vertex2710 long k=0;2711 for (i=0;i<edge4->nb();i++){2712 if (st[i]>=0){ // edge alone2713 if (i<nbe){2714 long i0=edge4->i(i);2715 ordre[i0] = vertices+i0;2716 long i1=edge4->j(i);2717 ordre[i1] = vertices+i1;2718 }2719 else {2720 k=k+1;2721 if (k<10) {2722 //print only 10 edges2723 printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i));2724 }2725 else if (k==10){2726 printf("Other lost boundary edges not shown...\n");2727 }2728 }2729 }2730 }2731 if(k) {2732 throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k));2733 }2734 2735 /* mesh generation with boundary points*/2736 long nbvb=0;2737 for (i=0;i<nbv;i++){2738 vertices[i].t=0;2739 vertices[i].vint=0;2740 if (ordre[i]) ordre[nbvb++]=ordre[i];2741 }2742 2743 Triangle* savetriangles=triangles;2744 long savenbt=nbt;2745 long savenbtx=nbtx;2746 SubDomain* savesubdomains=subdomains;2747 subdomains=0;2748 2749 long Nbtriafillhole=2*nbvb;2750 Triangle* triafillhole=new Triangle[Nbtriafillhole];2751 triangles = triafillhole;2752 2753 nbt=2;2754 nbtx= Nbtriafillhole;2755 2756 //Find a vertex that is not aligned with vertices 0 and 12757 for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;)2758 if (++i>=nbvb) {2759 throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned"));2760 }2761 //Move this vertex (i) to the 2d position in ordre2762 Exchange(ordre[2], ordre[i]);2763 2764 /*Reconstruct mesh beginning with 2 triangles*/2765 Vertex * v0=ordre[0], *v1=ordre[1];2766 2767 triangles[0](0) = NULL; // Infinite vertex2768 triangles[0](1) = v0;2769 triangles[0](2) = v1;2770 2771 triangles[1](0) = NULL;// Infinite vertex2772 triangles[1](2) = v0;2773 triangles[1](1) = v1;2774 const int e0 = OppositeEdge[0];2775 const int e1 = NextEdge[e0];2776 const int e2 = PreviousEdge[e0];2777 triangles[0].SetAdj2(e0, &triangles[1] ,e0);2778 triangles[0].SetAdj2(e1, &triangles[1] ,e2);2779 triangles[0].SetAdj2(e2, &triangles[1] ,e1);2780 2781 triangles[0].det = -1; // boundary triangles2782 triangles[1].det = -1; // boundary triangles2783 2784 triangles[0].SetTriangleContainingTheVertex();2785 triangles[1].SetTriangleContainingTheVertex();2786 2787 triangles[0].link=&triangles[1];2788 triangles[1].link=&triangles[0];2789 2790 if (!quadtree) delete quadtree; //ReInitialise;2791 quadtree = new QuadTree(this,0);2792 quadtree->Add(*v0);2793 quadtree->Add(*v1);2794 2795 // vertices are added one by one2796 long NbSwap=0;2797 for (int icount=2; icount<nbvb; icount++) {2798 Vertex *vi = ordre[icount];2799 Icoor2 dete[3];2800 Triangle *tcvi = FindTriangleContaining(vi->i,dete);2801 quadtree->Add(*vi);2802 AddVertex(*vi,tcvi,dete);2803 NbSwap += vi->Optim(1,1);2804 }2805 2806 //enforce the boundary2807 TriangleAdjacent ta(0,0);2808 long nbloss = 0,knbe=0;2809 for ( i = 0; i < nbe; i++){2810 if (st[i] >=0){ //edge alone => on border2811 Vertex &a=edges[i][0], &b=edges[i][1];2812 if (a.t && b.t){2813 knbe++;2814 if (ForceEdge(a,b,ta)<0) nbloss++;2815 }2816 }2817 }2818 if(nbloss) {2819 throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe));2820 }2821 2822 FindSubDomain(1);2823 // remove all the hole2824 // remove all the good sub domain2825 long krm =0;2826 for (i=0;i<nbt;i++){2827 if (triangles[i].link){ // remove triangles2828 krm++;2829 for (int j=0;j<3;j++){2830 TriangleAdjacent ta = triangles[i].Adj(j);2831 Triangle &tta = *(Triangle*)ta;2832 //if edge between remove and not remove2833 if(! tta.link){2834 // change the link of ta;2835 int ja = ta;2836 Vertex *v0= ta.EdgeVertex(0);2837 Vertex *v1= ta.EdgeVertex(1);2838 long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv);2839 2840 assert(st[k]>=0);2841 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3));2842 ta.SetLock();2843 st[k]=-2-st[k];2844 }2845 }2846 }2847 }2848 long NbTfillHoll =0;2849 for (i=0;i<nbt;i++){2850 if (triangles[i].link) {2851 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL);2852 triangles[i].color=-1;2853 }2854 else{2855 triangles[i].color= savenbt+ NbTfillHoll++;2856 }2857 }2858 assert(savenbt+NbTfillHoll<=savenbtx);2859 2860 // copy of the outside triangles in saveTriangles2861 for (i=0;i<nbt;i++){2862 if(triangles[i].color>=0) {2863 savetriangles[savenbt]=triangles[i];2864 savetriangles[savenbt].link=0;2865 savenbt++;2866 }2867 }2868 // gestion of the adj2869 k =0;2870 Triangle * tmax = triangles + nbt;2871 for (i=0;i<savenbt;i++) {2872 Triangle & ti = savetriangles[i];2873 for (int j=0;j<3;j++){2874 Triangle * ta = ti.TriangleAdj(j);2875 int aa = ti.NuEdgeTriangleAdj(j);2876 int lck = ti.Locked(j);2877 if (!ta) k++; // bug2878 else if ( ta >= triangles && ta < tmax){2879 ta= savetriangles + ta->color;2880 ti.SetAdj2(j,ta,aa);2881 if(lck) ti.SetLocked(j);2882 }2883 }2884 }2885 2886 // restore triangles;2887 nbt=savenbt;2888 nbtx=savenbtx;2889 delete [] triangles;2890 delete [] subdomains;2891 triangles = savetriangles;2892 subdomains = savesubdomains;2893 if (k) {2894 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k));2895 }2896 FindSubDomain();2897 2898 delete edge4;2899 delete [] st;2900 for (i=0;i<nbv;i++) quadtree->Add(vertices[i]);2901 2902 SetVertexFieldOn();2903 2904 /*Check requirements consistency*/2905 for (i=0;i<nbe;i++){2906 /*If the current mesh edge is on Geometry*/2907 if(edges[i].onGeometry){2908 for(int j=0;j<2;j++){2909 /*Go through the edges adjacent to current edge (if on the same curve)*/2910 if (!edges[i].adj[j]){2911 /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/2912 /*Check that the 2 vertices are on geometry AND required*/2913 if(!edges[i][j].onGeometry->IsRequiredVertex()){2914 printf("ReconstructExistingMesh error message: problem with the edge %i: [%i][%i]\n",Number(edges[i][j])+1,i,j);2915 printf(" This edge is on geometry (Edge of geometry %i)",Gh.Number(edges[i].onGeometry));2916 if (edges[i][j].onGeometry->OnGeomVertex())2917 printf(" Vertex %i\n",Gh.Number(edges[i][j].onGeometry->gv));2918 else if (edges[i][j].onGeometry->OnGeomEdge())2919 printf(" Edges %i\n",Gh.Number(edges[i][j].onGeometry->ge));2920 else2921 printf(" = %p\n",edges[i][j].onGeometry);2922 throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)"));2923 }2924 }2925 }2926 }2927 }2928 }2929 /*}}}1*/2930 2613 /*FUNCTION Triangles::GeomToTriangles0{{{1*/ 2931 2614 void Triangles::GeomToTriangles0(long inbvx,BamgOpts* bamgopts){ … … 3893 3576 long it,nbchange=0; 3894 3577 double lmax=0; 3895 for (it=0;it<nbt;it++) 3896 { 3578 for (it=0;it<nbt;it++){ 3897 3579 Triangle &t=triangles[it]; 3898 for (int j=0;j<3;j++) 3899 { 3580 for (int j=0;j<3;j++){ 3900 3581 Triangle &tt = *t.TriangleAdj(j); 3901 if ( ! &tt || it < Number(tt) && ( tt.link || t.link)) 3902 { 3582 if ( ! &tt || it < Number(tt) && ( tt.link || t.link)){ 3903 3583 Vertex &v0 = t[VerticesOfTriangularEdge[j][0]]; 3904 3584 Vertex &v1 = t[VerticesOfTriangularEdge[j][1]]; … … 3907 3587 double l = M(AB,AB); 3908 3588 lmax = Max(lmax,l); 3909 if(l> maxsubdiv2) 3910 {R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M3589 if(l> maxsubdiv2){ 3590 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3911 3591 double lc = M(AC,AC); 3912 3592 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; … … 3916 3596 v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y); 3917 3597 nbchange++; 3918 3598 } 3919 3599 M = v1; 3920 3600 l = M(AB,AB); 3921 3601 lmax = Max(lmax,l); 3922 if(l> maxsubdiv2) 3923 {R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M3602 if(l> maxsubdiv2){ 3603 R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M 3924 3604 double lc = M(AC,AC); 3925 3605 D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; … … 3929 3609 v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y); 3930 3610 nbchange++; 3931 } 3932 3933 3934 } 3935 } 3936 } 3611 } 3612 } 3613 } 3614 } 3937 3615 if(verbosity>3){ 3938 3616 printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5)); … … 4285 3963 } 4286 3964 /*}}}1*/ 3965 /*FUNCTION Triangles::ReconstructExistingMesh{{{1*/ 3966 void Triangles::ReconstructExistingMesh(){ 3967 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FillHoleInMesh)*/ 3968 3969 /*This routine reconstruct an existing mesh to make it CONVEX: 3970 * -all the holes are filled 3971 * -concave boundaries are filled 3972 * A convex mesh is required for a lot of operations. This is why every mesh 3973 * goes through this process. 3974 * This routine also generates mesh properties such as adjencies,... 3975 */ 3976 3977 /*Intermediary*/ 3978 int verbosity=0; 3979 3980 // generation of the integer coordinate 3981 3982 // find extrema coordinates of vertices pmin,pmax 3983 long i; 3984 if(verbosity>2) printf(" Reconstruct mesh of %i vertices\n",nbv); 3985 3986 //initialize ordre 3987 assert(ordre); 3988 for (i=0;i<nbv;i++) ordre[i]=0; 3989 3990 //Initialize NbSubDomains 3991 NbSubDomains =0; 3992 3993 /* generation of triangles adjacency*/ 3994 3995 //First add existing edges 3996 long kk =0; 3997 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 3998 for (i=0;i<nbe;i++){ 3999 kk=kk+(i==edge4->SortAndAdd(Number(edges[i][0]),Number(edges[i][1]))); 4000 } 4001 if (kk != nbe){ 4002 throw ErrorException(__FUNCT__,exprintf("There are %i double edges in the mesh",kk-nbe)); 4003 } 4004 4005 //Add edges of all triangles in existing mesh 4006 long* st = new long[nbt*3]; 4007 for (i=0;i<nbt*3;i++) st[i]=-1; 4008 for (i=0;i<nbt;i++){ 4009 for (int j=0;j<3;j++){ 4010 4011 //Add current triangle edge to edge4 4012 long k =edge4->SortAndAdd(Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]])); 4013 4014 long invisible=triangles[i].Hidden(j); 4015 4016 //If the edge has not been added to st, add it 4017 if(st[k]==-1) st[k]=3*i+j; 4018 4019 //If the edge already exists, add adjacency 4020 else if(st[k]>=0) { 4021 assert(!triangles[i].TriangleAdj(j)); 4022 assert(!triangles[st[k]/3].TriangleAdj((int) (st[k]%3))); 4023 4024 triangles[i].SetAdj2(j,triangles+st[k]/3,(int)(st[k]%3)); 4025 if (invisible) triangles[i].SetHidden(j); 4026 if (k<nbe) triangles[i].SetLocked(j); 4027 4028 //Make st[k] negative so that it will throw an error message if it is found again 4029 st[k]=-2-st[k]; 4030 } 4031 4032 //An edge belongs to 2 triangles 4033 else { 4034 throw ErrorException(__FUNCT__,exprintf("The edge (%i , %i) belongs to more than 2 triangles", 4035 Number(triangles[i][VerticesOfTriangularEdge[j][0]]),Number(triangles[i][VerticesOfTriangularEdge[j][1]]))); 4036 } 4037 } 4038 } 4039 4040 //Display info if required 4041 if(verbosity>5) { 4042 printf(" info of Mesh:\n"); 4043 printf(" - number of vertices = %i \n",nbv); 4044 printf(" - number of triangles = %i \n",nbt); 4045 printf(" - number of given edges = %i \n",nbe); 4046 printf(" - number of all edges = %i \n" ,edge4->nb()); 4047 printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); 4048 } 4049 4050 //check the consistency of edge[].adj and the geometrical required vertex 4051 long k=0; 4052 for (i=0;i<edge4->nb();i++){ 4053 if (st[i]>=0){ // edge alone 4054 if (i<nbe){ 4055 long i0=edge4->i(i); 4056 ordre[i0] = vertices+i0; 4057 long i1=edge4->j(i); 4058 ordre[i1] = vertices+i1; 4059 } 4060 else { 4061 k=k+1; 4062 if (k<10) { 4063 //print only 10 edges 4064 printf("Lost boundary edges %i : %i %i\n",i,edge4->i(i),edge4->j(i)); 4065 } 4066 else if (k==10){ 4067 printf("Other lost boundary edges not shown...\n"); 4068 } 4069 } 4070 } 4071 } 4072 if(k) { 4073 throw ErrorException(__FUNCT__,exprintf("%i boundary edges (from the geometry) are not defined as mesh edges",k)); 4074 } 4075 4076 /* mesh generation with boundary points*/ 4077 long nbvb=0; 4078 for (i=0;i<nbv;i++){ 4079 vertices[i].t=0; 4080 vertices[i].vint=0; 4081 if (ordre[i]) ordre[nbvb++]=ordre[i]; 4082 } 4083 4084 Triangle* savetriangles=triangles; 4085 long savenbt=nbt; 4086 long savenbtx=nbtx; 4087 SubDomain* savesubdomains=subdomains; 4088 subdomains=0; 4089 4090 long Nbtriafillhole=2*nbvb; 4091 Triangle* triafillhole=new Triangle[Nbtriafillhole]; 4092 triangles = triafillhole; 4093 4094 nbt=2; 4095 nbtx= Nbtriafillhole; 4096 4097 //Find a vertex that is not aligned with vertices 0 and 1 4098 for (i=2;det(ordre[0]->i,ordre[1]->i,ordre[i]->i)==0;) 4099 if (++i>=nbvb) { 4100 throw ErrorException(__FUNCT__,exprintf("ReconstructExistingMesh: All the vertices are aligned")); 4101 } 4102 //Move this vertex (i) to the 2d position in ordre 4103 Exchange(ordre[2], ordre[i]); 4104 4105 /*Reconstruct mesh beginning with 2 triangles*/ 4106 Vertex * v0=ordre[0], *v1=ordre[1]; 4107 4108 triangles[0](0) = NULL; // Infinite vertex 4109 triangles[0](1) = v0; 4110 triangles[0](2) = v1; 4111 4112 triangles[1](0) = NULL;// Infinite vertex 4113 triangles[1](2) = v0; 4114 triangles[1](1) = v1; 4115 const int e0 = OppositeEdge[0]; 4116 const int e1 = NextEdge[e0]; 4117 const int e2 = PreviousEdge[e0]; 4118 triangles[0].SetAdj2(e0, &triangles[1] ,e0); 4119 triangles[0].SetAdj2(e1, &triangles[1] ,e2); 4120 triangles[0].SetAdj2(e2, &triangles[1] ,e1); 4121 4122 triangles[0].det = -1; // boundary triangles 4123 triangles[1].det = -1; // boundary triangles 4124 4125 triangles[0].SetTriangleContainingTheVertex(); 4126 triangles[1].SetTriangleContainingTheVertex(); 4127 4128 triangles[0].link=&triangles[1]; 4129 triangles[1].link=&triangles[0]; 4130 4131 if (!quadtree) delete quadtree; //ReInitialise; 4132 quadtree = new QuadTree(this,0); 4133 quadtree->Add(*v0); 4134 quadtree->Add(*v1); 4135 4136 // vertices are added one by one 4137 long NbSwap=0; 4138 for (int icount=2; icount<nbvb; icount++) { 4139 Vertex *vi = ordre[icount]; 4140 Icoor2 dete[3]; 4141 Triangle *tcvi = FindTriangleContaining(vi->i,dete); 4142 quadtree->Add(*vi); 4143 AddVertex(*vi,tcvi,dete); 4144 NbSwap += vi->Optim(1,1); 4145 } 4146 4147 //enforce the boundary 4148 TriangleAdjacent ta(0,0); 4149 long nbloss = 0,knbe=0; 4150 for ( i = 0; i < nbe; i++){ 4151 if (st[i] >=0){ //edge alone => on border 4152 Vertex &a=edges[i][0], &b=edges[i][1]; 4153 if (a.t && b.t){ 4154 knbe++; 4155 if (ForceEdge(a,b,ta)<0) nbloss++; 4156 } 4157 } 4158 } 4159 if(nbloss) { 4160 throw ErrorException(__FUNCT__,exprintf("we lost %i existing edges other %i",nbloss,knbe)); 4161 } 4162 4163 FindSubDomain(1); 4164 // remove all the hole 4165 // remove all the good sub domain 4166 long krm =0; 4167 for (i=0;i<nbt;i++){ 4168 if (triangles[i].link){ // remove triangles 4169 krm++; 4170 for (int j=0;j<3;j++){ 4171 TriangleAdjacent ta = triangles[i].Adj(j); 4172 Triangle &tta = *(Triangle*)ta; 4173 //if edge between remove and not remove 4174 if(! tta.link){ 4175 // change the link of ta; 4176 int ja = ta; 4177 Vertex *v0= ta.EdgeVertex(0); 4178 Vertex *v1= ta.EdgeVertex(1); 4179 long k =edge4->SortAndAdd(v0?Number(v0):nbv,v1? Number(v1):nbv); 4180 4181 assert(st[k]>=0); 4182 tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); 4183 ta.SetLock(); 4184 st[k]=-2-st[k]; 4185 } 4186 } 4187 } 4188 } 4189 long NbTfillHoll =0; 4190 for (i=0;i<nbt;i++){ 4191 if (triangles[i].link) { 4192 triangles[i]=Triangle((Vertex *) NULL,(Vertex *) NULL,(Vertex *) NULL); 4193 triangles[i].color=-1; 4194 } 4195 else{ 4196 triangles[i].color= savenbt+ NbTfillHoll++; 4197 } 4198 } 4199 assert(savenbt+NbTfillHoll<=savenbtx); 4200 4201 // copy of the outside triangles in saveTriangles 4202 for (i=0;i<nbt;i++){ 4203 if(triangles[i].color>=0) { 4204 savetriangles[savenbt]=triangles[i]; 4205 savetriangles[savenbt].link=0; 4206 savenbt++; 4207 } 4208 } 4209 // gestion of the adj 4210 k =0; 4211 Triangle * tmax = triangles + nbt; 4212 for (i=0;i<savenbt;i++) { 4213 Triangle & ti = savetriangles[i]; 4214 for (int j=0;j<3;j++){ 4215 Triangle * ta = ti.TriangleAdj(j); 4216 int aa = ti.NuEdgeTriangleAdj(j); 4217 int lck = ti.Locked(j); 4218 if (!ta) k++; // bug 4219 else if ( ta >= triangles && ta < tmax){ 4220 ta= savetriangles + ta->color; 4221 ti.SetAdj2(j,ta,aa); 4222 if(lck) ti.SetLocked(j); 4223 } 4224 } 4225 } 4226 4227 // restore triangles; 4228 nbt=savenbt; 4229 nbtx=savenbtx; 4230 delete [] triangles; 4231 delete [] subdomains; 4232 triangles = savetriangles; 4233 subdomains = savesubdomains; 4234 if (k) { 4235 throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); 4236 } 4237 FindSubDomain(); 4238 4239 delete edge4; 4240 delete [] st; 4241 for (i=0;i<nbv;i++) quadtree->Add(vertices[i]); 4242 4243 SetVertexFieldOn(); 4244 4245 /*Check requirements consistency*/ 4246 for (i=0;i<nbe;i++){ 4247 /*If the current mesh edge is on Geometry*/ 4248 if(edges[i].onGeometry){ 4249 for(int j=0;j<2;j++){ 4250 /*Go through the edges adjacent to current edge (if on the same curve)*/ 4251 if (!edges[i].adj[j]){ 4252 /*The edge is on Geometry and does not have 2 adjacent edges... (not on a closed curve)*/ 4253 /*Check that the 2 vertices are on geometry AND required*/ 4254 if(!edges[i][j].onGeometry->IsRequiredVertex()){ 4255 printf("ReconstructExistingMesh error message: problem with the edge number %i: [%i %i]\n",i+1,Number(edges[i][0])+1,Number(edges[i][1])+1); 4256 printf("This edge is on geometrical edge number %i\n",Gh.Number(edges[i].onGeometry)+1); 4257 if (edges[i][j].onGeometry->OnGeomVertex()) 4258 printf("the vertex number %i of this edge is a geometric Vertex number %i\n",Gh.Number(edges[i][j].onGeometry->gv)+1); 4259 else if (edges[i][j].onGeometry->OnGeomEdge()) 4260 printf("the vertex number %i of this edge is a geometric Edge number %i\n",Gh.Number(edges[i][j].onGeometry->ge)+1); 4261 else 4262 printf("Its pointer is %p\n",edges[i][j].onGeometry); 4263 4264 printf("This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n"); 4265 throw ErrorException(__FUNCT__,exprintf("See above (might be cryptic...)")); 4266 } 4267 } 4268 } 4269 } 4270 } 4271 } 4272 /*}}}1*/ 4287 4273 /*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/ 4288 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) 4274 void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress){ 4289 4275 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ReNumberingTheTriangleBySubDomain)*/ 4290 4276 -
issm/trunk/src/c/objects/BamgOpts.cpp
r3277 r3279 37 37 38 38 if (bamgopts->coef==0) throw ErrorException(__FUNCT__,exprintf("'coef' should be positive")); 39 if (bamgopts->maxsubdiv< 1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >=1"));39 if (bamgopts->maxsubdiv<=1) throw ErrorException(__FUNCT__,exprintf("'maxsubdiv' should be >1")); 40 40 if (bamgopts->Hessiantype!=0 && bamgopts->Hessiantype!=1) throw ErrorException(__FUNCT__,exprintf("'Hessiantype' supported options are 0 and 1")); 41 41 if (bamgopts->Metrictype!=0 && bamgopts->Metrictype!=1 && bamgopts->Metrictype!=2) throw ErrorException(__FUNCT__,exprintf("'Metrictype' supported options are 0, 1 and 2")); -
issm/trunk/src/mex/Bamg/Bamg.cpp
r3277 r3279 30 30 FetchData(&bamggeom_in.NumEdges,mxGetField(BAMGGEOMETRY,0,"NumEdges")); 31 31 FetchData(&bamggeom_in.Edges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Edges")); 32 FetchData(&bamggeom_in.NumCorners,mxGetField(BAMGGEOMETRY,0,"NumCorners")); 33 FetchData(&bamggeom_in.Corners,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"Corners")); 34 FetchData(&bamggeom_in.NumRequiredVertices,mxGetField(BAMGGEOMETRY,0,"NumRequiredVertices")); 35 FetchData(&bamggeom_in.RequiredVertices,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"RequiredVertices")); 36 FetchData(&bamggeom_in.NumRequiredEdges,mxGetField(BAMGGEOMETRY,0,"NumRequiredEdges")); 37 FetchData(&bamggeom_in.RequiredEdges,NULL,NULL,mxGetField(BAMGGEOMETRY,0,"RequiredEdges")); 32 38 FetchData(&bamggeom_in.hVertices,&lines,&cols,mxGetField(BAMGGEOMETRY,0,"hVertices")); 33 39 FetchData(&bamggeom_in.NumCrackedEdges,mxGetField(BAMGGEOMETRY,0,"NumCrackedEdges"));
Note:
See TracChangeset
for help on using the changeset viewer.
