Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21864) @@ -11,5 +11,4 @@ this->gradation = 0; this->Hessiantype = 0; - this->MaxCornerAngle = 0; this->maxnbv = 0; this->maxsubdiv = 0; @@ -19,9 +18,7 @@ this->omega = 0; this->power = 0; - this->random = 0; this->verbose = 0; this->Crack = 0; - this->geometricalmetric = 0; this->KeepVertices = 0; this->splitcorners = 0; @@ -68,5 +65,4 @@ if (this->Crack!=0 && this->Crack!=1) _error_("'Crack' supported options are 0 and 1"); if (this->KeepVertices!=0 && this->KeepVertices!=1) _error_("'KeepVertices' supported options are 0 and 1"); - if (this->geometricalmetric!=0 && this->geometricalmetric!=1) _error_("'geometricalmetric' supported options are 0 and 1"); if (this->hmin<=0) _error_("'hmin' option should be >0"); Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.h (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.h (revision 21864) @@ -17,5 +17,4 @@ double gradation; int Hessiantype; - double MaxCornerAngle; int maxnbv; double maxsubdiv; @@ -25,10 +24,8 @@ double omega; double power; - bool random; int verbose; /*Flags*/ int Crack; - int geometricalmetric; int KeepVertices; int splitcorners; Index: /issm/trunk-jpl/src/c/bamg/Geometry.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 21864) @@ -187,10 +187,4 @@ } - //MaxCornerAngle - if (bamgopts->MaxCornerAngle){ - if(verbose>5) _printf_(" processing MaxCornerAngle\n"); - MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180; - } - //TangentAtEdges if (bamggeom->TangentAtEdges){ @@ -401,5 +395,4 @@ _printf_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")\n"); _printf_(" coefIcoor: " << coefIcoor << "\n"); - _printf_(" MaxCornerAngle: " << MaxCornerAngle << "\n"); return; @@ -419,5 +412,4 @@ nbcurves=0; subdomains=NULL; - MaxCornerAngle = 10*Pi/180; //default is 10 degres } /*}}}*/ @@ -614,17 +606,12 @@ } - // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) - if(ord==2) { + //all vertices provided in geometry are corners (ord = number of edges holding i) + vertices[i].SetCorner() ; + if(ord==2){ long n1 = head_v[i]; long n2 = next_p[n1]; long i1 = n1/2, i2 = n2/2; // edge number - long j1 = n1%2, j2 = n2%2; // vertex in the edge - float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; - float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; - float da12 = Abs(angle2-angle1); - if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) { - vertices[i].SetCorner() ; - } - // if the edge type/referencenumber a changing then is SetRequired(); + long j1 = n1%2, j2 = n2%2; // vertex in the edge + /* if the edge type/referencenumber a changing then is SetRequired();*/ if (edges[i1].type != edges[i2].type || edges[i1].Required()){ vertices[i].SetRequired(); @@ -633,7 +620,4 @@ vertices[i].SetRequired(); } - } - if(ord != 2) { - vertices[i].SetCorner(); } Index: /issm/trunk-jpl/src/c/bamg/Geometry.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Geometry.h (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/Geometry.h (revision 21864) @@ -32,5 +32,4 @@ R2 pmin,pmax; // domain extrema coordinates double coefIcoor; // coef to integer coordinates; - double MaxCornerAngle; //Constructor/Destructors Index: /issm/trunk-jpl/src/c/bamg/Mesh.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21864) @@ -980,55 +980,4 @@ /*Methods*/ - void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Metric.cpp/IntersectGeomMetric)*/ - - /*Get options*/ - double anisomax = bamgopts->anisomax; - double errg = bamgopts->errg; - - double ss[2]={0.00001,0.99999}; - double errC = 2*sqrt(2*errg); - double hmax = Gh.MaximalHmax(); - double hmin = Gh.MinimalHmin(); - - //check that hmax is positive - _assert_(hmax>0); - - //errC cannot be higher than 1 - if(errC>1) errC=1; - - //Set all vertices to "on" - SetVertexFieldOn(); - - //loop over all the vertices on edges - for(int i=0;iR1tg(s,tg); - double ht=hmax; - // err relative to the length of the edge - if (R1>1.0e-20) { - ht = Min(Max(errC/R1,hmin),hmax); - } - double hn=Min(hmax,ht*anisomax); - - if (ht<=0 || hn<=0){ - _error_("ht<=0 || hn<=0"); - } - EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg); - Metric MVp(Vp); - edges[i][j].m.IntersectWith(MVp); - } - } - // the problem is for the vertex on vertex - } - /*}}}*/ void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ // Hessiantype = 0 => H is computed using double L2 projection @@ -1231,10 +1180,8 @@ int i,j,k,kk,it,jt; int verbose=0; - double cutoffradian=10*Pi/180; /*Recover options*/ - if (bamgopts){ + if(bamgopts){ verbose=bamgopts->verbose; - cutoffradian=bamgopts->MaxCornerAngle*Pi/180; } @@ -1243,10 +1190,5 @@ //check that the mesh is not empty - if (nbt<=0 || nbv <=0 ) { - _error_("nbt or nbv is negative (Mesh empty?)"); - } - - //Gh is the geometry of the mesh (this), initialize MaxCornerAngle - if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian; + if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); /*Construction of the edges*/ @@ -1447,7 +1389,5 @@ //check that nbsubdomains is empty - if (nbsubdomains){ - _error_("nbsubdomains should be 0"); - } + if(nbsubdomains) _error_("nbsubdomains should be 0"); nbsubdomains=0; @@ -2649,49 +2589,49 @@ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/PreInit)*/ - /* initialize random seed: */ - srand(19999999); - /*Initialize fields*/ - NbRef=0; - quadtree=NULL; - nbv=0; - nbt=0; - nbe=0; - edges=NULL; - nbsubdomains=0; - subdomains=NULL; - maxnbv=maxnbv_in; - maxnbt=2 *maxnbv_in-2; - NbVertexOnBThVertex=0; - VertexOnBThVertex=NULL; - NbVertexOnBThEdge=0; - VertexOnBThEdge=NULL; - NbCrackedVertices=0; - CrackedVertices =NULL; - NbCrackedEdges =0; - CrackedEdges =NULL; - NbVerticesOnGeomVertex=0; - VerticesOnGeomVertex=NULL; - NbVerticesOnGeomEdge=0; - VerticesOnGeomEdge=NULL; + this->NbRef = 0; + this->quadtree = NULL; + this->nbv = 0; + this->nbt = 0; + this->nbe = 0; + this->edges = NULL; + this->nbsubdomains = 0; + this->subdomains = NULL; + this->maxnbv = maxnbv_in; + this->maxnbt = 2 *maxnbv_in-2; + this->NbVertexOnBThVertex = 0; + this->VertexOnBThVertex = NULL; + this->NbVertexOnBThEdge = 0; + this->VertexOnBThEdge = NULL; + this->NbCrackedVertices = 0; + this->CrackedVertices = NULL; + this->NbCrackedEdges = 0; + this->CrackedEdges = NULL; + this->NbVerticesOnGeomVertex = 0; + this->VerticesOnGeomVertex = NULL; + this->NbVerticesOnGeomEdge = 0; + this->VerticesOnGeomEdge = NULL; + + /*Initialize random seed*/ + this->randomseed = 1; /*Allocate if maxnbv_in>0*/ if(maxnbv_in){ - vertices=new BamgVertex[maxnbv]; - _assert_(vertices); - orderedvertices=new BamgVertex* [maxnbv]; - _assert_(orderedvertices); - triangles=new Triangle[maxnbt]; - _assert_(triangles); + this->vertices=new BamgVertex[this->maxnbv]; + this->orderedvertices=new BamgVertex* [this->maxnbv]; + this->triangles=new Triangle[this->maxnbt]; + _assert_(this->vertices); + _assert_(this->orderedvertices); + _assert_(this->triangles); } else{ - vertices=NULL; - orderedvertices=NULL; - triangles=NULL; - maxnbt=0; + this->vertices = NULL; + this->orderedvertices = NULL; + this->triangles = NULL; + this->maxnbt = 0; } } /*}}}*/ - void Mesh::Insert(bool random) {/*{{{*/ + void Mesh::Insert(){/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/Insert)*/ @@ -2740,6 +2680,5 @@ //Get Prime number const long PrimeNumber= BigPrimeNumber(nbv); - int temp = rand(); if(!random) temp = 756804110; - int k0=temp%nbv; + int k0=this->RandomNumber(nbv); //Build orderedvertices @@ -2825,5 +2764,5 @@ } /*}}}*/ - long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {/*{{{*/ + long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ @@ -2845,7 +2784,5 @@ /*construction of a random order*/ const long PrimeNumber= BigPrimeNumber(nbv) ; - //remainder of the division of rand() by nbvnew - int temp = rand(); if(!random) temp = 1098566905; - long k3 = temp%nbvnew; + long k3 = long(this->RandomNumber(nbvnew)); //loop over the new points for (int is3=0; is3random); + InsertNewPoints(nbvold,NbTSwap); } else Bh.CreateSingleVertexToTriangleConnectivity(); @@ -3146,5 +3083,5 @@ }// for triangle - if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break; + if (!InsertNewPoints(nbvold,NbTSwap)) break; for (i=nbtold;i3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); if (verbose>3) _printf_(" Inserting boundary points\n"); - Insert(bamgopts->random); + Insert(); //Force the boundary @@ -4722,5 +4659,5 @@ if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); if (verbose>3) _printf_(" Inserting boundary points\n"); - Insert(bamgopts->random); + Insert(); //Force the boundary @@ -4735,4 +4672,173 @@ NewPoints(BTh,bamgopts,KeepVertices) ; if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); + } + /*}}}*/ + int Mesh::RandomNumber(int max){/*{{{*/ + /* Generate a random number from 0 to max-1 using linear congruential + * random number generator*/ + + this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l; + return this->randomseed / (714025l / max + 1); + } /*}}}*/ + int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ + /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ + + int NbSwap =0; + if (!a.t || !b.t){ // the 2 vertex is in a mesh + _error_("!a.t || !b.t"); + } + int k=0; + taret=AdjacentTriangle(0,0); // erreur + + AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); + BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; + // we turn around a in the direct direction + + long long det2 = v2 ? det(*v2,a,b): -1 , det1; + if(v2) // normal case + det2 = det(*v2,a,b); + else { // no chance infini vertex try the next + tta= Previous(Adj(tta)); + v2 = tta.EdgeVertex(0); + vbegin =v2; + if (!v2){ + _error_("!v2"); + } + det2 = det(*v2,a,b); + } + + while (v2 != &b) { + AdjacentTriangle tc = Previous(Adj(tta)); + v1 = v2; + v2 = tc.EdgeVertex(0); + det1 = det2; + det2 = v2 ? det(*v2,a,b): det2; + + if((det1 < 0) && (det2 >0)) { + // try to force the edge + BamgVertex * va = &a, *vb = &b; + tc = Previous(tc); + if (!v1 || !v2){ + _error_("!v1 || !v2"); + } + long long detss = 0,l=0; + while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) + if(l++ > 10000000) { + _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); + } + BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); + if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ + tc.SetLock(); + a.Optim(1,0); + b.Optim(1,0); + taret = tc; + return NbSwap; + } + else + { + taret = tc; + return -2; // error boundary is crossing + } + } + tta = tc; + k++; + if (k>=2000){ + _error_("k>=2000"); + } + if ( vbegin == v2 ) return -1;// error + } + + tta.SetLock(); + taret=tta; + a.Optim(1,0); + b.Optim(1,0); + return NbSwap; + } + /*}}}*/ + int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ + /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ + // l'arete ta coupe l'arete pva pvb + // de cas apres le swap sa coupe toujours + // on cherche l'arete suivante + // on suppose que detsa >0 et detsb <0 + // attention la routine echange pva et pvb + + if(tt1.Locked()) return 0; // frontiere croise + + AdjacentTriangle tt2 = Adj(tt1); + Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent + short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles + if ( a1<0 || a1>=3 ){ + _error_("a1<0 || a1>=3"); + } + + BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; + BamgVertex & s1= (*t1)[OppositeVertex[a1]]; + BamgVertex & s2= (*t2)[OppositeVertex[a2]]; + + long long dets2 = det(*pva,*pvb,s2); + long long det1=t1->det , det2=t2->det ; + long long detT = det1+det2; + if ((det1<=0 ) || (det2<=0)){ + _error_("(det1<=0 ) || (det2<=0)"); + } + if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb + _error_("(detsa>=0) || (detsb<=0)"); + } + long long ndet1 = bamg::det(s1,sa,s2); + long long ndet2 = detT - ndet1; + + int ToSwap =0; //pas de swap + if ((ndet1 >0) && (ndet2 >0)) + { // on peut swaper + if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) + ToSwap =1; + else // swap alleatoire + if (BinaryRand()) + ToSwap =2; + } + if (ToSwap) NbSwap++, + bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); + + int ret=1; + + if (dets2 < 0) {// haut + dets1 = ToSwap ? dets1 : detsa ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = ToSwap ? dets1 : detsb ; + detsb = dets2; + //xxxx tt1 = ToSwap ? tt1 : Next(tt2); + if(!ToSwap) tt1 = Next(tt2); + } + else { // changement de direction + ret = -1; + Exchange(pva,pvb); + Exchange(detsa,detsb); + Exchange(dets1,dets2); + Exchange(tt1,tt2); + dets1=-dets1; + dets2=-dets2; + detsa=-detsa; + detsb=-detsb; + + if(ToSwap){ + if (dets2 < 0) {// haut + dets1 = (ToSwap ? dets1 : detsa) ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if(dets2 > 0){// bas + dets1 = (ToSwap ? dets1 : detsb) ; + detsb = dets2; + if(!ToSwap) tt1 = Next(tt2); + } + else {// on a fin ??? + tt1 = Next(tt2); + ret =0;} + } + + } + return ret; } /*}}}*/ @@ -4776,79 +4882,4 @@ return edge; } - } - /*}}}*/ - int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ - - int NbSwap =0; - if (!a.t || !b.t){ // the 2 vertex is in a mesh - _error_("!a.t || !b.t"); - } - int k=0; - taret=AdjacentTriangle(0,0); // erreur - - AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); - BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; - // we turn around a in the direct direction - - long long det2 = v2 ? det(*v2,a,b): -1 , det1; - if(v2) // normal case - det2 = det(*v2,a,b); - else { // no chance infini vertex try the next - tta= Previous(Adj(tta)); - v2 = tta.EdgeVertex(0); - vbegin =v2; - if (!v2){ - _error_("!v2"); - } - det2 = det(*v2,a,b); - } - - while (v2 != &b) { - AdjacentTriangle tc = Previous(Adj(tta)); - v1 = v2; - v2 = tc.EdgeVertex(0); - det1 = det2; - det2 = v2 ? det(*v2,a,b): det2; - - if((det1 < 0) && (det2 >0)) { - // try to force the edge - BamgVertex * va = &a, *vb = &b; - tc = Previous(tc); - if (!v1 || !v2){ - _error_("!v1 || !v2"); - } - long long detss = 0,l=0; - while ((SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) - if(l++ > 10000000) { - _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); - } - BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); - if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ - tc.SetLock(); - a.Optim(1,0); - b.Optim(1,0); - taret = tc; - return NbSwap; - } - else - { - taret = tc; - return -2; // error boundary is crossing - } - } - tta = tc; - k++; - if (k>=2000){ - _error_("k>=2000"); - } - if ( vbegin == v2 ) return -1;// error - } - - tta.SetLock(); - taret=tta; - a.Optim(1,0); - b.Optim(1,0); - return NbSwap; } /*}}}*/ @@ -4897,90 +4928,3 @@ } // end swap /*}}}*/ - int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ - // l'arete ta coupe l'arete pva pvb - // de cas apres le swap sa coupe toujours - // on cherche l'arete suivante - // on suppose que detsa >0 et detsb <0 - // attention la routine echange pva et pvb - - if(tt1.Locked()) return 0; // frontiere croise - - AdjacentTriangle tt2 = Adj(tt1); - Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent - short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles - if ( a1<0 || a1>=3 ){ - _error_("a1<0 || a1>=3"); - } - - BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; - BamgVertex & s1= (*t1)[OppositeVertex[a1]]; - BamgVertex & s2= (*t2)[OppositeVertex[a2]]; - - long long dets2 = det(*pva,*pvb,s2); - long long det1=t1->det , det2=t2->det ; - long long detT = det1+det2; - if ((det1<=0 ) || (det2<=0)){ - _error_("(det1<=0 ) || (det2<=0)"); - } - if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb - _error_("(detsa>=0) || (detsb<=0)"); - } - long long ndet1 = bamg::det(s1,sa,s2); - long long ndet2 = detT - ndet1; - - int ToSwap =0; //pas de swap - if ((ndet1 >0) && (ndet2 >0)) - { // on peut swaper - if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) - ToSwap =1; - else // swap alleatoire - if (BinaryRand()) - ToSwap =2; - } - if (ToSwap) NbSwap++, - bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); - - int ret=1; - - if (dets2 < 0) {// haut - dets1 = ToSwap ? dets1 : detsa ; - detsa = dets2; - tt1 = Previous(tt2) ;} - else if (dets2 > 0){// bas - dets1 = ToSwap ? dets1 : detsb ; - detsb = dets2; - //xxxx tt1 = ToSwap ? tt1 : Next(tt2); - if(!ToSwap) tt1 = Next(tt2); - } - else { // changement de direction - ret = -1; - Exchange(pva,pvb); - Exchange(detsa,detsb); - Exchange(dets1,dets2); - Exchange(tt1,tt2); - dets1=-dets1; - dets2=-dets2; - detsa=-detsa; - detsb=-detsb; - - if(ToSwap){ - if (dets2 < 0) {// haut - dets1 = (ToSwap ? dets1 : detsa) ; - detsa = dets2; - tt1 = Previous(tt2) ;} - else if(dets2 > 0){// bas - dets1 = (ToSwap ? dets1 : detsb) ; - detsb = dets2; - if(!ToSwap) tt1 = Next(tt2); - } - else {// on a fin ??? - tt1 = Next(tt2); - ret =0;} - } - - } - return ret; - } - /*}}}*/ } Index: /issm/trunk-jpl/src/c/bamg/Mesh.h =================================================================== --- /issm/trunk-jpl/src/c/bamg/Mesh.h (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/Mesh.h (revision 21864) @@ -40,4 +40,5 @@ double coefIcoor; // coef to integer ListofIntersectionTriangles lIntTria; + int randomseed; //used for random number generation long NbVerticesOnGeomVertex; @@ -77,5 +78,5 @@ R2 I2ToR2(const I2 & P) const; void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; - void Insert(bool random); + void Insert(); void Echo(void); void ForceBoundary(); @@ -91,5 +92,5 @@ void MaxSubDivision(double maxsubdiv); void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); - long InsertNewPoints(long nbvold,long & NbTSwap,bool random); + long InsertNewPoints(long nbvold,long & NbTSwap); void TrianglesRenumberBySubDomain(bool justcompress=false); void SmoothingVertex(int =3,double=0.3); @@ -113,6 +114,6 @@ void BuildMetric0(BamgOpts* bamgopts); void BuildMetric1(BamgOpts* bamgopts); - void AddGeometryMetric(BamgOpts* bamgopts); void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); + int RandomNumber(int max); void ReconstructExistingMesh(); @@ -143,4 +144,8 @@ void Triangulate(double* x,double* y,int nods); void Init(long); + int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; + int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , + AdjacentTriangle & tt1,long long & dets1, + long long & detsa,long long & detsb, int & nbswap); }; @@ -150,8 +155,5 @@ Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); - int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , - AdjacentTriangle & tt1,long long & dets1, - long long & detsa,long long & detsb, int & nbswap); - int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; + inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); Index: /issm/trunk-jpl/src/c/bamg/Triangle.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/Triangle.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/bamg/Triangle.cpp (revision 21864) @@ -257,5 +257,4 @@ long long detMinNew=Min(det1,det2); - // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test if (! OnSwap &&(detMinNew>0)) { OnSwap = detMin ==0; Index: /issm/trunk-jpl/src/c/classes/AmrBamg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21864) @@ -32,5 +32,4 @@ this->options->gradation = 1.5; this->options->Hessiantype = 0; - this->options->MaxCornerAngle = 1.e-12; this->options->maxnbv = 1e6; this->options->maxsubdiv = 10; @@ -40,8 +39,6 @@ this->options->omega = 1.8; this->options->power = 1; - this->options->random = 0; this->options->verbose = 0; this->options->Crack = 0; - this->options->geometricalmetric = 0; this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ this->options->splitcorners = 1; Index: /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 21863) +++ /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 21864) @@ -149,7 +149,4 @@ } - //Add geometry metric if provided - if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts); - //Smoothe metric BTh.SmoothMetric(bamgopts->gradation); Index: /issm/trunk-jpl/src/m/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.m (revision 21863) +++ /issm/trunk-jpl/src/m/mesh/bamg.m (revision 21864) @@ -22,5 +22,4 @@ % 1 -> use Green formula % - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) -% - MaxCornerAngle : maximum angle of corners in degree (default is 10) % - maxnbv : maximum number of vertices used to allocate memory (default is 10^6) % - maxsubdiv : maximum subdivision of exisiting elements (default is 10) @@ -34,5 +33,4 @@ % - power : power applied to the metric (default is 1) % - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) -% - geometricalmetric : take the geometry into account to generate the metric (default is 0) % - verbose : level of verbosity (default is 1) % @@ -330,5 +328,4 @@ bamg_options.hVertices=getfieldvalue(options,'hVertices',[]); bamg_options.KeepVertices=getfieldvalue(options,'KeepVertices',1); -bamg_options.MaxCornerAngle=getfieldvalue(options,'MaxCornerAngle',10.); bamg_options.maxnbv=getfieldvalue(options,'maxnbv',10^6); bamg_options.maxsubdiv=getfieldvalue(options,'maxsubdiv',10.); @@ -340,6 +337,4 @@ bamg_options.power=getfieldvalue(options,'power',1.); bamg_options.splitcorners=getfieldvalue(options,'splitcorners',1); -bamg_options.geometricalmetric=getfieldvalue(options,'geometricalmetric',0); -bamg_options.random=getfieldvalue(options,'rand',true); bamg_options.verbose=getfieldvalue(options,'verbose',1); %}}} Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 21863) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 21864) @@ -37,5 +37,4 @@ 1 -> use Green formula - KeepVertices : try to keep initial vertices when adaptation is done on an existing mesh (default 1) - - MaxCornerAngle : maximum angle of corners in degree (default is 10) - maxnbv : maximum number of vertices used to allocate memory (default is 10^6) - maxsubdiv : maximum subdivision of exisiting elements (default is 10) @@ -49,5 +48,4 @@ - power : power applied to the metric (default is 1) - splitcorners : split triangles whuch have 3 vertices on the outline (default is 1) - - geometricalmetric : take the geometry into account to generate the metric (default is 0) - verbose : level of verbosity (default is 1) @@ -303,5 +301,4 @@ bamg_options['hVertices']=options.getfieldvalue('hVertices',np.empty((0,1))) bamg_options['KeepVertices']=options.getfieldvalue('KeepVertices',1) - bamg_options['MaxCornerAngle']=options.getfieldvalue('MaxCornerAngle',10.) bamg_options['maxnbv']=options.getfieldvalue('maxnbv',10**6) bamg_options['maxsubdiv']=options.getfieldvalue('maxsubdiv',10.) @@ -313,6 +310,4 @@ bamg_options['power']=options.getfieldvalue('power',1.) bamg_options['splitcorners']=options.getfieldvalue('splitcorners',1) - bamg_options['geometricalmetric']=options.getfieldvalue('geometricalmetric',0) - bamg_options['random']=options.getfieldvalue('rand',True) bamg_options['verbose']=options.getfieldvalue('verbose',1) #}}} Index: /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp (revision 21863) +++ /issm/trunk-jpl/src/wrappers/matlab/io/FetchMatlabData.cpp (revision 21864) @@ -537,5 +537,4 @@ FetchData(&bamgopts->gradation,mxGetField(dataref,0,"gradation")); FetchData(&bamgopts->Hessiantype,mxGetField(dataref,0,"Hessiantype")); - FetchData(&bamgopts->MaxCornerAngle,mxGetField(dataref,0,"MaxCornerAngle")); FetchData(&bamgopts->maxnbv,mxGetField(dataref,0,"maxnbv")); FetchData(&bamgopts->maxsubdiv,mxGetField(dataref,0,"maxsubdiv")); @@ -545,9 +544,7 @@ FetchData(&bamgopts->omega,mxGetField(dataref,0,"omega")); FetchData(&bamgopts->power,mxGetField(dataref,0,"power")); - FetchData(&bamgopts->random,mxGetField(dataref,0,"random")); FetchData(&bamgopts->verbose,mxGetField(dataref,0,"verbose")); FetchData(&bamgopts->Crack,mxGetField(dataref,0,"Crack")); - FetchData(&bamgopts->geometricalmetric,mxGetField(dataref,0,"geometricalmetric")); FetchData(&bamgopts->KeepVertices,mxGetField(dataref,0,"KeepVertices")); FetchData(&bamgopts->splitcorners,mxGetField(dataref,0,"splitcorners")); Index: /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp (revision 21863) +++ /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp (revision 21864) @@ -692,5 +692,4 @@ FetchData(&bamgopts->gradation,PyDict_GetItemString(py_dict,"gradation")); FetchData(&bamgopts->Hessiantype,PyDict_GetItemString(py_dict,"Hessiantype")); - FetchData(&bamgopts->MaxCornerAngle,PyDict_GetItemString(py_dict,"MaxCornerAngle")); FetchData(&bamgopts->maxnbv,PyDict_GetItemString(py_dict,"maxnbv")); FetchData(&bamgopts->maxsubdiv,PyDict_GetItemString(py_dict,"maxsubdiv")); @@ -700,9 +699,7 @@ FetchData(&bamgopts->omega,PyDict_GetItemString(py_dict,"omega")); FetchData(&bamgopts->power,PyDict_GetItemString(py_dict,"power")); - FetchData(&bamgopts->random,PyDict_GetItemString(py_dict,"random")); FetchData(&bamgopts->verbose,PyDict_GetItemString(py_dict,"verbose")); FetchData(&bamgopts->Crack,PyDict_GetItemString(py_dict,"Crack")); - FetchData(&bamgopts->geometricalmetric,PyDict_GetItemString(py_dict,"geometricalmetric")); FetchData(&bamgopts->KeepVertices,PyDict_GetItemString(py_dict,"KeepVertices")); FetchData(&bamgopts->splitcorners,PyDict_GetItemString(py_dict,"splitcorners"));