![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21864
- Timestamp:
- 07/25/17 13:03:02 (8 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 13 edited
-
c/bamg/BamgOpts.cpp (modified) (3 diffs)
-
c/bamg/BamgOpts.h (modified) (2 diffs)
-
c/bamg/Geometry.cpp (modified) (5 diffs)
-
c/bamg/Geometry.h (modified) (1 diff)
-
c/bamg/Mesh.cpp (modified) (16 diffs)
-
c/bamg/Mesh.h (modified) (6 diffs)
-
c/bamg/Triangle.cpp (modified) (1 diff)
-
c/classes/AmrBamg.cpp (modified) (2 diffs)
-
c/modules/Bamgx/Bamgx.cpp (modified) (1 diff)
-
m/mesh/bamg.m (modified) (4 diffs)
-
m/mesh/bamg.py (modified) (4 diffs)
-
wrappers/matlab/io/FetchMatlabData.cpp (modified) (2 diffs)
-
wrappers/python/io/FetchPythonData.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/bamg/BamgOpts.cpp
r21838 r21864 11 11 this->gradation = 0; 12 12 this->Hessiantype = 0; 13 this->MaxCornerAngle = 0;14 13 this->maxnbv = 0; 15 14 this->maxsubdiv = 0; … … 19 18 this->omega = 0; 20 19 this->power = 0; 21 this->random = 0;22 20 this->verbose = 0; 23 21 24 22 this->Crack = 0; 25 this->geometricalmetric = 0;26 23 this->KeepVertices = 0; 27 24 this->splitcorners = 0; … … 68 65 if (this->Crack!=0 && this->Crack!=1) _error_("'Crack' supported options are 0 and 1"); 69 66 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");71 67 72 68 if (this->hmin<=0) _error_("'hmin' option should be >0"); -
issm/trunk-jpl/src/c/bamg/BamgOpts.h
r21838 r21864 17 17 double gradation; 18 18 int Hessiantype; 19 double MaxCornerAngle;20 19 int maxnbv; 21 20 double maxsubdiv; … … 25 24 double omega; 26 25 double power; 27 bool random;28 26 int verbose; 29 27 30 28 /*Flags*/ 31 29 int Crack; 32 int geometricalmetric;33 30 int KeepVertices; 34 31 int splitcorners; -
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r21838 r21864 187 187 } 188 188 189 //MaxCornerAngle190 if (bamgopts->MaxCornerAngle){191 if(verbose>5) _printf_(" processing MaxCornerAngle\n");192 MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;193 }194 195 189 //TangentAtEdges 196 190 if (bamggeom->TangentAtEdges){ … … 401 395 _printf_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")\n"); 402 396 _printf_(" coefIcoor: " << coefIcoor << "\n"); 403 _printf_(" MaxCornerAngle: " << MaxCornerAngle << "\n");404 397 405 398 return; … … 419 412 nbcurves=0; 420 413 subdomains=NULL; 421 MaxCornerAngle = 10*Pi/180; //default is 10 degres422 414 } 423 415 /*}}}*/ … … 614 606 } 615 607 616 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) 617 if(ord==2) { 608 //all vertices provided in geometry are corners (ord = number of edges holding i) 609 vertices[i].SetCorner() ; 610 if(ord==2){ 618 611 long n1 = head_v[i]; 619 612 long n2 = next_p[n1]; 620 613 long i1 = n1/2, i2 = n2/2; // edge number 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(); 614 long j1 = n1%2, j2 = n2%2; // vertex in the edge 615 /* if the edge type/referencenumber a changing then is SetRequired();*/ 629 616 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 630 617 vertices[i].SetRequired(); … … 633 620 vertices[i].SetRequired(); 634 621 } 635 }636 if(ord != 2) {637 vertices[i].SetCorner();638 622 } 639 623 -
issm/trunk-jpl/src/c/bamg/Geometry.h
r21838 r21864 32 32 R2 pmin,pmax; // domain extrema coordinates 33 33 double coefIcoor; // coef to integer coordinates; 34 double MaxCornerAngle;35 34 36 35 //Constructor/Destructors -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r21838 r21864 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 positive995 _assert_(hmax>0);996 997 //errC cannot be higher than 1998 if(errC>1) errC=1;999 1000 //Set all vertices to "on"1001 SetVertexFieldOn();1002 1003 //loop over all the vertices on edges1004 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 edge1017 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 vertex1031 }1032 /*}}}*/1033 982 void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ 1034 983 // Hessiantype = 0 => H is computed using double L2 projection … … 1231 1180 int i,j,k,kk,it,jt; 1232 1181 int verbose=0; 1233 double cutoffradian=10*Pi/180;1234 1182 1235 1183 /*Recover options*/ 1236 if 1184 if(bamgopts){ 1237 1185 verbose=bamgopts->verbose; 1238 cutoffradian=bamgopts->MaxCornerAngle*Pi/180;1239 1186 } 1240 1187 … … 1243 1190 1244 1191 //check that the mesh is not 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; 1192 if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); 1251 1193 1252 1194 /*Construction of the edges*/ … … 1447 1389 1448 1390 //check that nbsubdomains is empty 1449 if (nbsubdomains){ 1450 _error_("nbsubdomains should be 0"); 1451 } 1391 if(nbsubdomains) _error_("nbsubdomains should be 0"); 1452 1392 nbsubdomains=0; 1453 1393 … … 2649 2589 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2650 2590 2651 /* initialize random seed: */2652 srand(19999999);2653 2654 2591 /*Initialize fields*/ 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; 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; 2677 2617 2678 2618 /*Allocate if maxnbv_in>0*/ 2679 2619 if(maxnbv_in){ 2680 vertices=new BamgVertex[maxnbv];2681 _assert_(vertices);2682 orderedvertices=new BamgVertex* [maxnbv];2683 _assert_(orderedvertices);2684 triangles=new Triangle[maxnbt];2685 _assert_(t riangles);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_(this->triangles); 2686 2626 } 2687 2627 else{ 2688 vertices=NULL;2689 orderedvertices=NULL;2690 t riangles=NULL;2691 maxnbt=0;2628 this->vertices = NULL; 2629 this->orderedvertices = NULL; 2630 this->triangles = NULL; 2631 this->maxnbt = 0; 2692 2632 } 2693 2633 } 2694 2634 /*}}}*/ 2695 void Mesh::Insert( bool random){/*{{{*/2635 void Mesh::Insert(){/*{{{*/ 2696 2636 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2697 2637 … … 2740 2680 //Get Prime number 2741 2681 const long PrimeNumber= BigPrimeNumber(nbv); 2742 int temp = rand(); if(!random) temp = 756804110; 2743 int k0=temp%nbv; 2682 int k0=this->RandomNumber(nbv); 2744 2683 2745 2684 //Build orderedvertices … … 2825 2764 } 2826 2765 /*}}}*/ 2827 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap ,bool random) {/*{{{*/2766 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ 2828 2767 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2829 2768 … … 2845 2784 /*construction of a random order*/ 2846 2785 const long PrimeNumber= BigPrimeNumber(nbv) ; 2847 //remainder of the division of rand() by nbvnew 2848 int temp = rand(); if(!random) temp = 1098566905; 2849 long k3 = temp%nbvnew; 2786 long k3 = long(this->RandomNumber(nbvnew)); 2850 2787 //loop over the new points 2851 2788 for (int is3=0; is3<nbvnew; is3++){ … … 3088 3025 //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); 3089 3026 Bh.CreateSingleVertexToTriangleConnectivity(); 3090 InsertNewPoints(nbvold,NbTSwap ,bamgopts->random);3027 InsertNewPoints(nbvold,NbTSwap); 3091 3028 } 3092 3029 else Bh.CreateSingleVertexToTriangleConnectivity(); … … 3146 3083 }// for triangle 3147 3084 3148 if (!InsertNewPoints(nbvold,NbTSwap ,bamgopts->random)) break;3085 if (!InsertNewPoints(nbvold,NbTSwap)) break; 3149 3086 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3150 3087 Headt = nbt; // empty list … … 4095 4032 4096 4033 /*Insert Vertices*/ 4097 Insert( true);4034 Insert(); 4098 4035 } 4099 4036 /*}}}*/ … … 4396 4333 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4397 4334 if (verbose>3) _printf_(" Inserting boundary points\n"); 4398 Insert( bamgopts->random);4335 Insert(); 4399 4336 4400 4337 //Force the boundary … … 4722 4659 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4723 4660 if (verbose>3) _printf_(" Inserting boundary points\n"); 4724 Insert( bamgopts->random);4661 Insert(); 4725 4662 4726 4663 //Force the boundary … … 4735 4672 NewPoints(BTh,bamgopts,KeepVertices) ; 4736 4673 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 congruential 4678 * 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 mesh 4688 _error_("!a.t || !b.t"); 4689 } 4690 int k=0; 4691 taret=AdjacentTriangle(0,0); // erreur 4692 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 direction 4696 4697 long long det2 = v2 ? det(*v2,a,b): -1 , det1; 4698 if(v2) // normal case 4699 det2 = det(*v2,a,b); 4700 else { // no chance infini vertex try the next 4701 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 edge 4719 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 else 4738 { 4739 taret = tc; 4740 return -2; // error boundary is crossing 4741 } 4742 } 4743 tta = tc; 4744 k++; 4745 if (k>=2000){ 4746 _error_("k>=2000"); 4747 } 4748 if ( vbegin == v2 ) return -1;// error 4749 } 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 pvb 4761 // de cas apres le swap sa coupe toujours 4762 // on cherche l'arete suivante 4763 // on suppose que detsa >0 et detsb <0 4764 // attention la routine echange pva et pvb 4765 4766 if(tt1.Locked()) return 0; // frontiere croise 4767 4768 AdjacentTriangle tt2 = Adj(tt1); 4769 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 4770 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 4771 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,bb 4786 _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 swap 4792 if ((ndet1 >0) && (ndet2 >0)) 4793 { // on peut swaper 4794 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 4795 ToSwap =1; 4796 else // swap alleatoire 4797 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) {// haut 4806 dets1 = ToSwap ? dets1 : detsa ; 4807 detsa = dets2; 4808 tt1 = Previous(tt2) ;} 4809 else if (dets2 > 0){// bas 4810 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 direction 4816 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) {// haut 4828 dets1 = (ToSwap ? dets1 : detsa) ; 4829 detsa = dets2; 4830 tt1 = Previous(tt2) ;} 4831 else if(dets2 > 0){// bas 4832 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; 4737 4843 } 4738 4844 /*}}}*/ … … 4776 4882 return edge; 4777 4883 } 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 mesh4785 _error_("!a.t || !b.t");4786 }4787 int k=0;4788 taret=AdjacentTriangle(0,0); // erreur4789 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 direction4793 4794 long long det2 = v2 ? det(*v2,a,b): -1 , det1;4795 if(v2) // normal case4796 det2 = det(*v2,a,b);4797 else { // no chance infini vertex try the next4798 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 edge4816 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 else4835 {4836 taret = tc;4837 return -2; // error boundary is crossing4838 }4839 }4840 tta = tc;4841 k++;4842 if (k>=2000){4843 _error_("k>=2000");4844 }4845 if ( vbegin == v2 ) return -1;// error4846 }4847 4848 tta.SetLock();4849 taret=tta;4850 a.Optim(1,0);4851 b.Optim(1,0);4852 return NbSwap;4853 4884 } 4854 4885 /*}}}*/ … … 4897 4928 } // end swap 4898 4929 /*}}}*/ 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 pvb4902 // de cas apres le swap sa coupe toujours4903 // on cherche l'arete suivante4904 // on suppose que detsa >0 et detsb <04905 // attention la routine echange pva et pvb4906 4907 if(tt1.Locked()) return 0; // frontiere croise4908 4909 AdjacentTriangle tt2 = Adj(tt1);4910 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent4911 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles4912 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,bb4927 _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 swap4933 if ((ndet1 >0) && (ndet2 >0))4934 { // on peut swaper4935 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))4936 ToSwap =1;4937 else // swap alleatoire4938 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) {// haut4947 dets1 = ToSwap ? dets1 : detsa ;4948 detsa = dets2;4949 tt1 = Previous(tt2) ;}4950 else if (dets2 > 0){// bas4951 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 direction4957 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) {// haut4969 dets1 = (ToSwap ? dets1 : detsa) ;4970 detsa = dets2;4971 tt1 = Previous(tt2) ;}4972 else if(dets2 > 0){// bas4973 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 /*}}}*/4986 4930 } -
issm/trunk-jpl/src/c/bamg/Mesh.h
r21838 r21864 40 40 double coefIcoor; // coef to integer 41 41 ListofIntersectionTriangles lIntTria; 42 int randomseed; //used for random number generation 42 43 43 44 long NbVerticesOnGeomVertex; … … 77 78 R2 I2ToR2(const I2 & P) const; 78 79 void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; 79 void Insert( bool random);80 void Insert(); 80 81 void Echo(void); 81 82 void ForceBoundary(); … … 91 92 void MaxSubDivision(double maxsubdiv); 92 93 void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); 93 long InsertNewPoints(long nbvold,long & NbTSwap ,bool random);94 long InsertNewPoints(long nbvold,long & NbTSwap); 94 95 void TrianglesRenumberBySubDomain(bool justcompress=false); 95 96 void SmoothingVertex(int =3,double=0.3); … … 113 114 void BuildMetric0(BamgOpts* bamgopts); 114 115 void BuildMetric1(BamgOpts* bamgopts); 115 void AddGeometryMetric(BamgOpts* bamgopts);116 116 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 117 int RandomNumber(int max); 117 118 void ReconstructExistingMesh(); 118 119 … … 143 144 void Triangulate(double* x,double* y,int nods); 144 145 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); 145 150 }; 146 151 … … 150 155 Triangle *t2,short a2, 151 156 BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); 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) ; 157 156 158 inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ 157 159 return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); -
issm/trunk-jpl/src/c/bamg/Triangle.cpp
r21838 r21864 257 257 258 258 long long detMinNew=Min(det1,det2); 259 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test260 259 if (! OnSwap &&(detMinNew>0)) { 261 260 OnSwap = detMin ==0; -
issm/trunk-jpl/src/c/classes/AmrBamg.cpp
r21838 r21864 32 32 this->options->gradation = 1.5; 33 33 this->options->Hessiantype = 0; 34 this->options->MaxCornerAngle = 1.e-12;35 34 this->options->maxnbv = 1e6; 36 35 this->options->maxsubdiv = 10; … … 40 39 this->options->omega = 1.8; 41 40 this->options->power = 1; 42 this->options->random = 0;43 41 this->options->verbose = 0; 44 42 this->options->Crack = 0; 45 this->options->geometricalmetric = 0;46 43 this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ 47 44 this->options->splitcorners = 1; -
issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
r21838 r21864 149 149 } 150 150 151 //Add geometry metric if provided152 if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts);153 154 151 //Smoothe metric 155 152 BTh.SmoothMetric(bamgopts->gradation); -
issm/trunk-jpl/src/m/mesh/bamg.m
r21453 r21864 22 22 % 1 -> use Green formula 23 23 % - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) 24 % - MaxCornerAngle : maximum angle of corners in degree (default is 10)25 24 % - maxnbv : maximum number of vertices used to allocate memory (default is 10^6) 26 25 % - maxsubdiv : maximum subdivision of exisiting elements (default is 10) … … 34 33 % - power : power applied to the metric (default is 1) 35 34 % - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) 36 % - geometricalmetric : take the geometry into account to generate the metric (default is 0)37 35 % - verbose : level of verbosity (default is 1) 38 36 % … … 330 328 bamg_options.hVertices=getfieldvalue(options,'hVertices',[]); 331 329 bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1); 332 bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10.);333 330 bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6); 334 331 bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10.); … … 340 337 bamg_options.power=getfieldvalue(options,'power',1.); 341 338 bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1); 342 bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0);343 bamg_options.random=getfieldvalue(options,'rand',true);344 339 bamg_options.verbose=getfieldvalue(options,'verbose',1); 345 340 %}}} -
issm/trunk-jpl/src/m/mesh/bamg.py
r21303 r21864 37 37 1 -> use Green formula 38 38 - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) 39 - MaxCornerAngle : maximum angle of corners in degree (default is 10)40 39 - maxnbv : maximum number of vertices used to allocate memory (default is 10^6) 41 40 - maxsubdiv : maximum subdivision of exisiting elements (default is 10) … … 49 48 - power : power applied to the metric (default is 1) 50 49 - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) 51 - geometricalmetric : take the geometry into account to generate the metric (default is 0)52 50 - verbose : level of verbosity (default is 1) 53 51 … … 303 301 bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1))) 304 302 bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1) 305 bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.)306 303 bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6) 307 304 bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.) … … 313 310 bamg_options['power']=options.getfieldvalue('power',1.) 314 311 bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1) 315 bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0)316 bamg_options['random']=options.getfieldvalue('rand',True)317 312 bamg_options['verbose']=options.getfieldvalue('verbose',1) 318 313 #}}} -
issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp
r16968 r21864 537 537 FetchData(&bamgopts->gradation,mxGetField(dataref,0,"gradation")); 538 538 FetchData(&bamgopts->Hessiantype,mxGetField(dataref,0,"Hessiantype")); 539 FetchData(&bamgopts->MaxCornerAngle,mxGetField(dataref,0,"MaxCornerAngle"));540 539 FetchData(&bamgopts->maxnbv,mxGetField(dataref,0,"maxnbv")); 541 540 FetchData(&bamgopts->maxsubdiv,mxGetField(dataref,0,"maxsubdiv")); … … 545 544 FetchData(&bamgopts->omega,mxGetField(dataref,0,"omega")); 546 545 FetchData(&bamgopts->power,mxGetField(dataref,0,"power")); 547 FetchData(&bamgopts->random,mxGetField(dataref,0,"random"));548 546 FetchData(&bamgopts->verbose,mxGetField(dataref,0,"verbose")); 549 547 550 548 FetchData(&bamgopts->Crack,mxGetField(dataref,0,"Crack")); 551 FetchData(&bamgopts->geometricalmetric,mxGetField(dataref,0,"geometricalmetric"));552 549 FetchData(&bamgopts->KeepVertices,mxGetField(dataref,0,"KeepVertices")); 553 550 FetchData(&bamgopts->splitcorners,mxGetField(dataref,0,"splitcorners")); -
issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp
r21201 r21864 692 692 FetchData(&bamgopts->gradation,PyDict_GetItemString(py_dict,"gradation")); 693 693 FetchData(&bamgopts->Hessiantype,PyDict_GetItemString(py_dict,"Hessiantype")); 694 FetchData(&bamgopts->MaxCornerAngle,PyDict_GetItemString(py_dict,"MaxCornerAngle"));695 694 FetchData(&bamgopts->maxnbv,PyDict_GetItemString(py_dict,"maxnbv")); 696 695 FetchData(&bamgopts->maxsubdiv,PyDict_GetItemString(py_dict,"maxsubdiv")); … … 700 699 FetchData(&bamgopts->omega,PyDict_GetItemString(py_dict,"omega")); 701 700 FetchData(&bamgopts->power,PyDict_GetItemString(py_dict,"power")); 702 FetchData(&bamgopts->random,PyDict_GetItemString(py_dict,"random"));703 701 FetchData(&bamgopts->verbose,PyDict_GetItemString(py_dict,"verbose")); 704 702 705 703 FetchData(&bamgopts->Crack,PyDict_GetItemString(py_dict,"Crack")); 706 FetchData(&bamgopts->geometricalmetric,PyDict_GetItemString(py_dict,"geometricalmetric"));707 704 FetchData(&bamgopts->KeepVertices,PyDict_GetItemString(py_dict,"KeepVertices")); 708 705 FetchData(&bamgopts->splitcorners,PyDict_GetItemString(py_dict,"splitcorners"));
Note:
See TracChangeset
for help on using the changeset viewer.
