Index: /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp =================================================================== --- /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21837) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.cpp (revision 21838) @@ -11,4 +11,5 @@ this->gradation = 0; this->Hessiantype = 0; + this->MaxCornerAngle = 0; this->maxnbv = 0; this->maxsubdiv = 0; @@ -18,7 +19,9 @@ this->omega = 0; this->power = 0; + this->random = 0; this->verbose = 0; this->Crack = 0; + this->geometricalmetric = 0; this->KeepVertices = 0; this->splitcorners = 0; @@ -65,4 +68,5 @@ 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 21837) +++ /issm/trunk-jpl/src/c/bamg/BamgOpts.h (revision 21838) @@ -17,4 +17,5 @@ double gradation; int Hessiantype; + double MaxCornerAngle; int maxnbv; double maxsubdiv; @@ -24,8 +25,10 @@ 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 21837) +++ /issm/trunk-jpl/src/c/bamg/Geometry.cpp (revision 21838) @@ -187,4 +187,10 @@ } + //MaxCornerAngle + if (bamgopts->MaxCornerAngle){ + if(verbose>5) _printf_(" processing MaxCornerAngle\n"); + MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180; + } + //TangentAtEdges if (bamggeom->TangentAtEdges){ @@ -395,4 +401,5 @@ _printf_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")\n"); _printf_(" coefIcoor: " << coefIcoor << "\n"); + _printf_(" MaxCornerAngle: " << MaxCornerAngle << "\n"); return; @@ -412,4 +419,5 @@ nbcurves=0; subdomains=NULL; + MaxCornerAngle = 10*Pi/180; //default is 10 degres } /*}}}*/ @@ -606,12 +614,17 @@ } - //all vertices provided in geometry are corners (ord = number of edges holding i) - vertices[i].SetCorner() ; - if(ord==2){ + // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) + 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 - /* if the edge type/referencenumber a changing then is SetRequired();*/ + 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(); if (edges[i1].type != edges[i2].type || edges[i1].Required()){ vertices[i].SetRequired(); @@ -620,4 +633,7 @@ 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 21837) +++ /issm/trunk-jpl/src/c/bamg/Geometry.h (revision 21838) @@ -32,4 +32,5 @@ 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 21837) +++ /issm/trunk-jpl/src/c/bamg/Mesh.cpp (revision 21838) @@ -980,4 +980,55 @@ /*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 @@ -1180,8 +1231,10 @@ 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; } @@ -1190,5 +1243,10 @@ //check that the mesh is not empty - if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh 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; /*Construction of the edges*/ @@ -1389,5 +1447,7 @@ //check that nbsubdomains is empty - if(nbsubdomains) _error_("nbsubdomains should be 0"); + if (nbsubdomains){ + _error_("nbsubdomains should be 0"); + } nbsubdomains=0; @@ -2589,49 +2649,49 @@ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/PreInit)*/ + /* initialize random seed: */ + srand(19999999); + /*Initialize fields*/ - 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; + 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; /*Allocate if maxnbv_in>0*/ if(maxnbv_in){ - 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); + vertices=new BamgVertex[maxnbv]; + _assert_(vertices); + orderedvertices=new BamgVertex* [maxnbv]; + _assert_(orderedvertices); + triangles=new Triangle[maxnbt]; + _assert_(triangles); } else{ - this->vertices = NULL; - this->orderedvertices = NULL; - this->triangles = NULL; - this->maxnbt = 0; + vertices=NULL; + orderedvertices=NULL; + triangles=NULL; + maxnbt=0; } } /*}}}*/ - void Mesh::Insert(){/*{{{*/ + void Mesh::Insert(bool random) {/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/Insert)*/ @@ -2680,5 +2740,6 @@ //Get Prime number const long PrimeNumber= BigPrimeNumber(nbv); - int k0=this->RandomNumber(nbv); + int temp = rand(); if(!random) temp = 756804110; + int k0=temp%nbv; //Build orderedvertices @@ -2764,5 +2825,5 @@ } /*}}}*/ - long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ + long Mesh::InsertNewPoints(long nbvold,long & NbTSwap,bool random) {/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ @@ -2784,5 +2845,7 @@ /*construction of a random order*/ const long PrimeNumber= BigPrimeNumber(nbv) ; - long k3 = long(this->RandomNumber(nbvnew)); + //remainder of the division of rand() by nbvnew + int temp = rand(); if(!random) temp = 1098566905; + long k3 = temp%nbvnew; //loop over the new points for (int is3=0; is3random); } else Bh.CreateSingleVertexToTriangleConnectivity(); @@ -3083,5 +3146,5 @@ }// for triangle - if (!InsertNewPoints(nbvold,NbTSwap)) break; + if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break; for (i=nbtold;i3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); if (verbose>3) _printf_(" Inserting boundary points\n"); - Insert(); + Insert(bamgopts->random); //Force the boundary @@ -4659,5 +4722,5 @@ if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); if (verbose>3) _printf_(" Inserting boundary points\n"); - Insert(); + Insert(bamgopts->random); //Force the boundary @@ -4672,173 +4735,4 @@ 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; } /*}}}*/ @@ -4882,4 +4776,79 @@ 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; } /*}}}*/ @@ -4928,3 +4897,90 @@ } // 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 21837) +++ /issm/trunk-jpl/src/c/bamg/Mesh.h (revision 21838) @@ -40,5 +40,4 @@ double coefIcoor; // coef to integer ListofIntersectionTriangles lIntTria; - int randomseed; //used for random number generation long NbVerticesOnGeomVertex; @@ -78,5 +77,5 @@ R2 I2ToR2(const I2 & P) const; void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; - void Insert(); + void Insert(bool random); void Echo(void); void ForceBoundary(); @@ -92,5 +91,5 @@ void MaxSubDivision(double maxsubdiv); void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); - long InsertNewPoints(long nbvold,long & NbTSwap); + long InsertNewPoints(long nbvold,long & NbTSwap,bool random); void TrianglesRenumberBySubDomain(bool justcompress=false); void SmoothingVertex(int =3,double=0.3); @@ -114,6 +113,6 @@ void BuildMetric0(BamgOpts* bamgopts); void BuildMetric1(BamgOpts* bamgopts); + void AddGeometryMetric(BamgOpts* bamgopts); void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); - int RandomNumber(int max); void ReconstructExistingMesh(); @@ -144,8 +143,4 @@ 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); }; @@ -155,5 +150,8 @@ 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 21837) +++ /issm/trunk-jpl/src/c/bamg/Triangle.cpp (revision 21838) @@ -257,4 +257,5 @@ 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 21837) +++ /issm/trunk-jpl/src/c/classes/AmrBamg.cpp (revision 21838) @@ -32,4 +32,5 @@ this->options->gradation = 1.5; this->options->Hessiantype = 0; + this->options->MaxCornerAngle = 1.e-12; this->options->maxnbv = 1e6; this->options->maxsubdiv = 10; @@ -39,6 +40,8 @@ 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 21837) +++ /issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp (revision 21838) @@ -149,4 +149,7 @@ } + //Add geometry metric if provided + if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts); + //Smoothe metric BTh.SmoothMetric(bamgopts->gradation);