![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21838
- Timestamp:
- 07/24/17 10:37:43 (8 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 9 edited
-
bamg/BamgOpts.cpp (modified) (3 diffs)
-
bamg/BamgOpts.h (modified) (2 diffs)
-
bamg/Geometry.cpp (modified) (5 diffs)
-
bamg/Geometry.h (modified) (1 diff)
-
bamg/Mesh.cpp (modified) (16 diffs)
-
bamg/Mesh.h (modified) (6 diffs)
-
bamg/Triangle.cpp (modified) (1 diff)
-
classes/AmrBamg.cpp (modified) (2 diffs)
-
modules/Bamgx/Bamgx.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/bamg/BamgOpts.cpp
r21837 r21838 11 11 this->gradation = 0; 12 12 this->Hessiantype = 0; 13 this->MaxCornerAngle = 0; 13 14 this->maxnbv = 0; 14 15 this->maxsubdiv = 0; … … 18 19 this->omega = 0; 19 20 this->power = 0; 21 this->random = 0; 20 22 this->verbose = 0; 21 23 22 24 this->Crack = 0; 25 this->geometricalmetric = 0; 23 26 this->KeepVertices = 0; 24 27 this->splitcorners = 0; … … 65 68 if (this->Crack!=0 && this->Crack!=1) _error_("'Crack' supported options are 0 and 1"); 66 69 if (this->KeepVertices!=0 && this->KeepVertices!=1) _error_("'KeepVertices' supported options are 0 and 1"); 70 if (this->geometricalmetric!=0 && this->geometricalmetric!=1) _error_("'geometricalmetric' supported options are 0 and 1"); 67 71 68 72 if (this->hmin<=0) _error_("'hmin' option should be >0"); -
issm/trunk-jpl/src/c/bamg/BamgOpts.h
r21837 r21838 17 17 double gradation; 18 18 int Hessiantype; 19 double MaxCornerAngle; 19 20 int maxnbv; 20 21 double maxsubdiv; … … 24 25 double omega; 25 26 double power; 27 bool random; 26 28 int verbose; 27 29 28 30 /*Flags*/ 29 31 int Crack; 32 int geometricalmetric; 30 33 int KeepVertices; 31 34 int splitcorners; -
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r21837 r21838 187 187 } 188 188 189 //MaxCornerAngle 190 if (bamgopts->MaxCornerAngle){ 191 if(verbose>5) _printf_(" processing MaxCornerAngle\n"); 192 MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180; 193 } 194 189 195 //TangentAtEdges 190 196 if (bamggeom->TangentAtEdges){ … … 395 401 _printf_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")\n"); 396 402 _printf_(" coefIcoor: " << coefIcoor << "\n"); 403 _printf_(" MaxCornerAngle: " << MaxCornerAngle << "\n"); 397 404 398 405 return; … … 412 419 nbcurves=0; 413 420 subdomains=NULL; 421 MaxCornerAngle = 10*Pi/180; //default is 10 degres 414 422 } 415 423 /*}}}*/ … … 606 614 } 607 615 608 //all vertices provided in geometry are corners (ord = number of edges holding i) 609 vertices[i].SetCorner() ; 610 if(ord==2){ 616 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) 617 if(ord==2) { 611 618 long n1 = head_v[i]; 612 619 long n2 = next_p[n1]; 613 620 long i1 = n1/2, i2 = n2/2; // edge number 614 long j1 = n1%2, j2 = n2%2; // vertex in the edge 615 /* if the edge type/referencenumber a changing then is SetRequired();*/ 621 long j1 = n1%2, j2 = n2%2; // vertex in the edge 622 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; 623 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; 624 float da12 = Abs(angle2-angle1); 625 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) { 626 vertices[i].SetCorner() ; 627 } 628 // if the edge type/referencenumber a changing then is SetRequired(); 616 629 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 617 630 vertices[i].SetRequired(); … … 620 633 vertices[i].SetRequired(); 621 634 } 635 } 636 if(ord != 2) { 637 vertices[i].SetCorner(); 622 638 } 623 639 -
issm/trunk-jpl/src/c/bamg/Geometry.h
r21837 r21838 32 32 R2 pmin,pmax; // domain extrema coordinates 33 33 double coefIcoor; // coef to integer coordinates; 34 double MaxCornerAngle; 34 35 35 36 //Constructor/Destructors -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r21837 r21838 980 980 981 981 /*Methods*/ 982 void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/*{{{*/ 983 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectGeomMetric)*/ 984 985 /*Get options*/ 986 double anisomax = bamgopts->anisomax; 987 double errg = bamgopts->errg; 988 989 double ss[2]={0.00001,0.99999}; 990 double errC = 2*sqrt(2*errg); 991 double hmax = Gh.MaximalHmax(); 992 double hmin = Gh.MinimalHmin(); 993 994 //check that hmax is positive 995 _assert_(hmax>0); 996 997 //errC cannot be higher than 1 998 if(errC>1) errC=1; 999 1000 //Set all vertices to "on" 1001 SetVertexFieldOn(); 1002 1003 //loop over all the vertices on edges 1004 for(int i=0;i<nbe;i++){ 1005 for (int j=0;j<2;j++){ 1006 1007 BamgVertex V; 1008 VertexOnGeom GV; 1009 Gh.ProjectOnCurve(edges[i],ss[j],V,GV); 1010 1011 GeomEdge* eg = GV; 1012 double s = GV; 1013 R2 tg; 1014 double R1= eg->R1tg(s,tg); 1015 double ht=hmax; 1016 // err relative to the length of the edge 1017 if (R1>1.0e-20) { 1018 ht = Min(Max(errC/R1,hmin),hmax); 1019 } 1020 double hn=Min(hmax,ht*anisomax); 1021 1022 if (ht<=0 || hn<=0){ 1023 _error_("ht<=0 || hn<=0"); 1024 } 1025 EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg); 1026 Metric MVp(Vp); 1027 edges[i][j].m.IntersectWith(MVp); 1028 } 1029 } 1030 // the problem is for the vertex on vertex 1031 } 1032 /*}}}*/ 982 1033 void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ 983 1034 // Hessiantype = 0 => H is computed using double L2 projection … … 1180 1231 int i,j,k,kk,it,jt; 1181 1232 int verbose=0; 1233 double cutoffradian=10*Pi/180; 1182 1234 1183 1235 /*Recover options*/ 1184 if (bamgopts){1236 if (bamgopts){ 1185 1237 verbose=bamgopts->verbose; 1238 cutoffradian=bamgopts->MaxCornerAngle*Pi/180; 1186 1239 } 1187 1240 … … 1190 1243 1191 1244 //check that the mesh is not empty 1192 if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); 1245 if (nbt<=0 || nbv <=0 ) { 1246 _error_("nbt or nbv is negative (Mesh empty?)"); 1247 } 1248 1249 //Gh is the geometry of the mesh (this), initialize MaxCornerAngle 1250 if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian; 1193 1251 1194 1252 /*Construction of the edges*/ … … 1389 1447 1390 1448 //check that nbsubdomains is empty 1391 if(nbsubdomains) _error_("nbsubdomains should be 0"); 1449 if (nbsubdomains){ 1450 _error_("nbsubdomains should be 0"); 1451 } 1392 1452 nbsubdomains=0; 1393 1453 … … 2589 2649 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2590 2650 2651 /* initialize random seed: */ 2652 srand(19999999); 2653 2591 2654 /*Initialize fields*/ 2592 this->NbRef = 0; 2593 this->quadtree = NULL; 2594 this->nbv = 0; 2595 this->nbt = 0; 2596 this->nbe = 0; 2597 this->edges = NULL; 2598 this->nbsubdomains = 0; 2599 this->subdomains = NULL; 2600 this->maxnbv = maxnbv_in; 2601 this->maxnbt = 2 *maxnbv_in-2; 2602 this->NbVertexOnBThVertex = 0; 2603 this->VertexOnBThVertex = NULL; 2604 this->NbVertexOnBThEdge = 0; 2605 this->VertexOnBThEdge = NULL; 2606 this->NbCrackedVertices = 0; 2607 this->CrackedVertices = NULL; 2608 this->NbCrackedEdges = 0; 2609 this->CrackedEdges = NULL; 2610 this->NbVerticesOnGeomVertex = 0; 2611 this->VerticesOnGeomVertex = NULL; 2612 this->NbVerticesOnGeomEdge = 0; 2613 this->VerticesOnGeomEdge = NULL; 2614 2615 /*Initialize random seed*/ 2616 this->randomseed = 1; 2655 NbRef=0; 2656 quadtree=NULL; 2657 nbv=0; 2658 nbt=0; 2659 nbe=0; 2660 edges=NULL; 2661 nbsubdomains=0; 2662 subdomains=NULL; 2663 maxnbv=maxnbv_in; 2664 maxnbt=2 *maxnbv_in-2; 2665 NbVertexOnBThVertex=0; 2666 VertexOnBThVertex=NULL; 2667 NbVertexOnBThEdge=0; 2668 VertexOnBThEdge=NULL; 2669 NbCrackedVertices=0; 2670 CrackedVertices =NULL; 2671 NbCrackedEdges =0; 2672 CrackedEdges =NULL; 2673 NbVerticesOnGeomVertex=0; 2674 VerticesOnGeomVertex=NULL; 2675 NbVerticesOnGeomEdge=0; 2676 VerticesOnGeomEdge=NULL; 2617 2677 2618 2678 /*Allocate if maxnbv_in>0*/ 2619 2679 if(maxnbv_in){ 2620 this->vertices=new BamgVertex[this->maxnbv];2621 this->orderedvertices=new BamgVertex* [this->maxnbv];2622 this->triangles=new Triangle[this->maxnbt];2623 _assert_(this->vertices);2624 _assert_(this->orderedvertices);2625 _assert_(t his->triangles);2680 vertices=new BamgVertex[maxnbv]; 2681 _assert_(vertices); 2682 orderedvertices=new BamgVertex* [maxnbv]; 2683 _assert_(orderedvertices); 2684 triangles=new Triangle[maxnbt]; 2685 _assert_(triangles); 2626 2686 } 2627 2687 else{ 2628 this->vertices =NULL;2629 this->orderedvertices =NULL;2630 t his->triangles =NULL;2631 this->maxnbt =0;2688 vertices=NULL; 2689 orderedvertices=NULL; 2690 triangles=NULL; 2691 maxnbt=0; 2632 2692 } 2633 2693 } 2634 2694 /*}}}*/ 2635 void Mesh::Insert( ){/*{{{*/2695 void Mesh::Insert(bool random) {/*{{{*/ 2636 2696 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2637 2697 … … 2680 2740 //Get Prime number 2681 2741 const long PrimeNumber= BigPrimeNumber(nbv); 2682 int k0=this->RandomNumber(nbv); 2742 int temp = rand(); if(!random) temp = 756804110; 2743 int k0=temp%nbv; 2683 2744 2684 2745 //Build orderedvertices … … 2764 2825 } 2765 2826 /*}}}*/ 2766 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap ) {/*{{{*/2827 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {/*{{{*/ 2767 2828 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2768 2829 … … 2784 2845 /*construction of a random order*/ 2785 2846 const long PrimeNumber= BigPrimeNumber(nbv) ; 2786 long k3 = long(this->RandomNumber(nbvnew)); 2847 //remainder of the division of rand() by nbvnew 2848 int temp = rand(); if(!random) temp = 1098566905; 2849 long k3 = temp%nbvnew; 2787 2850 //loop over the new points 2788 2851 for (int is3=0; is3<nbvnew; is3++){ … … 3025 3088 //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); 3026 3089 Bh.CreateSingleVertexToTriangleConnectivity(); 3027 InsertNewPoints(nbvold,NbTSwap );3090 InsertNewPoints(nbvold,NbTSwap,bamgopts->random); 3028 3091 } 3029 3092 else Bh.CreateSingleVertexToTriangleConnectivity(); … … 3083 3146 }// for triangle 3084 3147 3085 if (!InsertNewPoints(nbvold,NbTSwap )) break;3148 if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break; 3086 3149 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3087 3150 Headt = nbt; // empty list … … 4032 4095 4033 4096 /*Insert Vertices*/ 4034 Insert( );4097 Insert(true); 4035 4098 } 4036 4099 /*}}}*/ … … 4333 4396 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4334 4397 if (verbose>3) _printf_(" Inserting boundary points\n"); 4335 Insert( );4398 Insert(bamgopts->random); 4336 4399 4337 4400 //Force the boundary … … 4659 4722 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4660 4723 if (verbose>3) _printf_(" Inserting boundary points\n"); 4661 Insert( );4724 Insert(bamgopts->random); 4662 4725 4663 4726 //Force the boundary … … 4672 4735 NewPoints(BTh,bamgopts,KeepVertices) ; 4673 4736 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); 4674 }4675 /*}}}*/4676 int Mesh::RandomNumber(int max){/*{{{*/4677 /* Generate a random number from 0 to max-1 using linear congruential4678 * random number generator*/4679 4680 this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l;4681 return this->randomseed / (714025l / max + 1);4682 } /*}}}*/4683 int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/4684 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/4685 4686 int NbSwap =0;4687 if (!a.t || !b.t){ // the 2 vertex is in a mesh4688 _error_("!a.t || !b.t");4689 }4690 int k=0;4691 taret=AdjacentTriangle(0,0); // erreur4692 4693 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);4694 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;4695 // we turn around a in the direct direction4696 4697 long long det2 = v2 ? det(*v2,a,b): -1 , det1;4698 if(v2) // normal case4699 det2 = det(*v2,a,b);4700 else { // no chance infini vertex try the next4701 tta= Previous(Adj(tta));4702 v2 = tta.EdgeVertex(0);4703 vbegin =v2;4704 if (!v2){4705 _error_("!v2");4706 }4707 det2 = det(*v2,a,b);4708 }4709 4710 while (v2 != &b) {4711 AdjacentTriangle tc = Previous(Adj(tta));4712 v1 = v2;4713 v2 = tc.EdgeVertex(0);4714 det1 = det2;4715 det2 = v2 ? det(*v2,a,b): det2;4716 4717 if((det1 < 0) && (det2 >0)) {4718 // try to force the edge4719 BamgVertex * va = &a, *vb = &b;4720 tc = Previous(tc);4721 if (!v1 || !v2){4722 _error_("!v1 || !v2");4723 }4724 long long detss = 0,l=0;4725 while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))4726 if(l++ > 10000000) {4727 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");4728 }4729 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);4730 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){4731 tc.SetLock();4732 a.Optim(1,0);4733 b.Optim(1,0);4734 taret = tc;4735 return NbSwap;4736 }4737 else4738 {4739 taret = tc;4740 return -2; // error boundary is crossing4741 }4742 }4743 tta = tc;4744 k++;4745 if (k>=2000){4746 _error_("k>=2000");4747 }4748 if ( vbegin == v2 ) return -1;// error4749 }4750 4751 tta.SetLock();4752 taret=tta;4753 a.Optim(1,0);4754 b.Optim(1,0);4755 return NbSwap;4756 }4757 /*}}}*/4758 int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/4759 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/4760 // l'arete ta coupe l'arete pva pvb4761 // de cas apres le swap sa coupe toujours4762 // on cherche l'arete suivante4763 // on suppose que detsa >0 et detsb <04764 // attention la routine echange pva et pvb4765 4766 if(tt1.Locked()) return 0; // frontiere croise4767 4768 AdjacentTriangle tt2 = Adj(tt1);4769 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent4770 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles4771 if ( a1<0 || a1>=3 ){4772 _error_("a1<0 || a1>=3");4773 }4774 4775 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];4776 BamgVertex & s1= (*t1)[OppositeVertex[a1]];4777 BamgVertex & s2= (*t2)[OppositeVertex[a2]];4778 4779 long long dets2 = det(*pva,*pvb,s2);4780 long long det1=t1->det , det2=t2->det ;4781 long long detT = det1+det2;4782 if ((det1<=0 ) || (det2<=0)){4783 _error_("(det1<=0 ) || (det2<=0)");4784 }4785 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb4786 _error_("(detsa>=0) || (detsb<=0)");4787 }4788 long long ndet1 = bamg::det(s1,sa,s2);4789 long long ndet2 = detT - ndet1;4790 4791 int ToSwap =0; //pas de swap4792 if ((ndet1 >0) && (ndet2 >0))4793 { // on peut swaper4794 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))4795 ToSwap =1;4796 else // swap alleatoire4797 if (BinaryRand())4798 ToSwap =2;4799 }4800 if (ToSwap) NbSwap++,4801 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);4802 4803 int ret=1;4804 4805 if (dets2 < 0) {// haut4806 dets1 = ToSwap ? dets1 : detsa ;4807 detsa = dets2;4808 tt1 = Previous(tt2) ;}4809 else if (dets2 > 0){// bas4810 dets1 = ToSwap ? dets1 : detsb ;4811 detsb = dets2;4812 //xxxx tt1 = ToSwap ? tt1 : Next(tt2);4813 if(!ToSwap) tt1 = Next(tt2);4814 }4815 else { // changement de direction4816 ret = -1;4817 Exchange(pva,pvb);4818 Exchange(detsa,detsb);4819 Exchange(dets1,dets2);4820 Exchange(tt1,tt2);4821 dets1=-dets1;4822 dets2=-dets2;4823 detsa=-detsa;4824 detsb=-detsb;4825 4826 if(ToSwap){4827 if (dets2 < 0) {// haut4828 dets1 = (ToSwap ? dets1 : detsa) ;4829 detsa = dets2;4830 tt1 = Previous(tt2) ;}4831 else if(dets2 > 0){// bas4832 dets1 = (ToSwap ? dets1 : detsb) ;4833 detsb = dets2;4834 if(!ToSwap) tt1 = Next(tt2);4835 }4836 else {// on a fin ???4837 tt1 = Next(tt2);4838 ret =0;}4839 }4840 4841 }4842 return ret;4843 4737 } 4844 4738 /*}}}*/ … … 4882 4776 return edge; 4883 4777 } 4778 } 4779 /*}}}*/ 4780 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4781 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4782 4783 int NbSwap =0; 4784 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4785 _error_("!a.t || !b.t"); 4786 } 4787 int k=0; 4788 taret=AdjacentTriangle(0,0); // erreur 4789 4790 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4791 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4792 // we turn around a in the direct direction 4793 4794 long long det2 = v2 ? det(*v2,a,b): -1 , det1; 4795 if(v2) // normal case 4796 det2 = det(*v2,a,b); 4797 else { // no chance infini vertex try the next 4798 tta= Previous(Adj(tta)); 4799 v2 = tta.EdgeVertex(0); 4800 vbegin =v2; 4801 if (!v2){ 4802 _error_("!v2"); 4803 } 4804 det2 = det(*v2,a,b); 4805 } 4806 4807 while (v2 != &b) { 4808 AdjacentTriangle tc = Previous(Adj(tta)); 4809 v1 = v2; 4810 v2 = tc.EdgeVertex(0); 4811 det1 = det2; 4812 det2 = v2 ? det(*v2,a,b): det2; 4813 4814 if((det1 < 0) && (det2 >0)) { 4815 // try to force the edge 4816 BamgVertex * va = &a, *vb = &b; 4817 tc = Previous(tc); 4818 if (!v1 || !v2){ 4819 _error_("!v1 || !v2"); 4820 } 4821 long long detss = 0,l=0; 4822 while ((SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4823 if(l++ > 10000000) { 4824 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4825 } 4826 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4827 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4828 tc.SetLock(); 4829 a.Optim(1,0); 4830 b.Optim(1,0); 4831 taret = tc; 4832 return NbSwap; 4833 } 4834 else 4835 { 4836 taret = tc; 4837 return -2; // error boundary is crossing 4838 } 4839 } 4840 tta = tc; 4841 k++; 4842 if (k>=2000){ 4843 _error_("k>=2000"); 4844 } 4845 if ( vbegin == v2 ) return -1;// error 4846 } 4847 4848 tta.SetLock(); 4849 taret=tta; 4850 a.Optim(1,0); 4851 b.Optim(1,0); 4852 return NbSwap; 4884 4853 } 4885 4854 /*}}}*/ … … 4928 4897 } // end swap 4929 4898 /*}}}*/ 4899 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ 4900 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ 4901 // l'arete ta coupe l'arete pva pvb 4902 // de cas apres le swap sa coupe toujours 4903 // on cherche l'arete suivante 4904 // on suppose que detsa >0 et detsb <0 4905 // attention la routine echange pva et pvb 4906 4907 if(tt1.Locked()) return 0; // frontiere croise 4908 4909 AdjacentTriangle tt2 = Adj(tt1); 4910 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 4911 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 4912 if ( a1<0 || a1>=3 ){ 4913 _error_("a1<0 || a1>=3"); 4914 } 4915 4916 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 4917 BamgVertex & s1= (*t1)[OppositeVertex[a1]]; 4918 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 4919 4920 long long dets2 = det(*pva,*pvb,s2); 4921 long long det1=t1->det , det2=t2->det ; 4922 long long detT = det1+det2; 4923 if ((det1<=0 ) || (det2<=0)){ 4924 _error_("(det1<=0 ) || (det2<=0)"); 4925 } 4926 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 4927 _error_("(detsa>=0) || (detsb<=0)"); 4928 } 4929 long long ndet1 = bamg::det(s1,sa,s2); 4930 long long ndet2 = detT - ndet1; 4931 4932 int ToSwap =0; //pas de swap 4933 if ((ndet1 >0) && (ndet2 >0)) 4934 { // on peut swaper 4935 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 4936 ToSwap =1; 4937 else // swap alleatoire 4938 if (BinaryRand()) 4939 ToSwap =2; 4940 } 4941 if (ToSwap) NbSwap++, 4942 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 4943 4944 int ret=1; 4945 4946 if (dets2 < 0) {// haut 4947 dets1 = ToSwap ? dets1 : detsa ; 4948 detsa = dets2; 4949 tt1 = Previous(tt2) ;} 4950 else if (dets2 > 0){// bas 4951 dets1 = ToSwap ? dets1 : detsb ; 4952 detsb = dets2; 4953 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 4954 if(!ToSwap) tt1 = Next(tt2); 4955 } 4956 else { // changement de direction 4957 ret = -1; 4958 Exchange(pva,pvb); 4959 Exchange(detsa,detsb); 4960 Exchange(dets1,dets2); 4961 Exchange(tt1,tt2); 4962 dets1=-dets1; 4963 dets2=-dets2; 4964 detsa=-detsa; 4965 detsb=-detsb; 4966 4967 if(ToSwap){ 4968 if (dets2 < 0) {// haut 4969 dets1 = (ToSwap ? dets1 : detsa) ; 4970 detsa = dets2; 4971 tt1 = Previous(tt2) ;} 4972 else if(dets2 > 0){// bas 4973 dets1 = (ToSwap ? dets1 : detsb) ; 4974 detsb = dets2; 4975 if(!ToSwap) tt1 = Next(tt2); 4976 } 4977 else {// on a fin ??? 4978 tt1 = Next(tt2); 4979 ret =0;} 4980 } 4981 4982 } 4983 return ret; 4984 } 4985 /*}}}*/ 4930 4986 } -
issm/trunk-jpl/src/c/bamg/Mesh.h
r21837 r21838 40 40 double coefIcoor; // coef to integer 41 41 ListofIntersectionTriangles lIntTria; 42 int randomseed; //used for random number generation43 42 44 43 long NbVerticesOnGeomVertex; … … 78 77 R2 I2ToR2(const I2 & P) const; 79 78 void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; 80 void Insert( );79 void Insert(bool random); 81 80 void Echo(void); 82 81 void ForceBoundary(); … … 92 91 void MaxSubDivision(double maxsubdiv); 93 92 void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); 94 long InsertNewPoints(long nbvold,long & NbTSwap );93 long InsertNewPoints(long nbvold,long & NbTSwap,bool random); 95 94 void TrianglesRenumberBySubDomain(bool justcompress=false); 96 95 void SmoothingVertex(int =3,double=0.3); … … 114 113 void BuildMetric0(BamgOpts* bamgopts); 115 114 void BuildMetric1(BamgOpts* bamgopts); 115 void AddGeometryMetric(BamgOpts* bamgopts); 116 116 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 117 int RandomNumber(int max);118 117 void ReconstructExistingMesh(); 119 118 … … 144 143 void Triangulate(double* x,double* y,int nods); 145 144 void Init(long); 146 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ;147 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,148 AdjacentTriangle & tt1,long long & dets1,149 long long & detsa,long long & detsb, int & nbswap);150 145 }; 151 146 … … 155 150 Triangle *t2,short a2, 156 151 BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); 157 152 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , 153 AdjacentTriangle & tt1,long long & dets1, 154 long long & detsa,long long & detsb, int & nbswap); 155 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; 158 156 inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ 159 157 return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); -
issm/trunk-jpl/src/c/bamg/Triangle.cpp
r21837 r21838 257 257 258 258 long long detMinNew=Min(det1,det2); 259 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test 259 260 if (! OnSwap &&(detMinNew>0)) { 260 261 OnSwap = detMin ==0; -
issm/trunk-jpl/src/c/classes/AmrBamg.cpp
r21837 r21838 32 32 this->options->gradation = 1.5; 33 33 this->options->Hessiantype = 0; 34 this->options->MaxCornerAngle = 1.e-12; 34 35 this->options->maxnbv = 1e6; 35 36 this->options->maxsubdiv = 10; … … 39 40 this->options->omega = 1.8; 40 41 this->options->power = 1; 42 this->options->random = 0; 41 43 this->options->verbose = 0; 42 44 this->options->Crack = 0; 45 this->options->geometricalmetric = 0; 43 46 this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ 44 47 this->options->splitcorners = 1; -
issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
r21837 r21838 149 149 } 150 150 151 //Add geometry metric if provided 152 if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts); 153 151 154 //Smoothe metric 152 155 BTh.SmoothMetric(bamgopts->gradation);
Note:
See TracChangeset
for help on using the changeset viewer.
