Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.cpp (revision 22085) @@ -34,5 +34,5 @@ } /*}}}*/ - Icoor2 & AdjacentTriangle::det() const {/*{{{*/ + long long & AdjacentTriangle::det() const {/*{{{*/ return t->det; } Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.h (revision 22085) @@ -37,5 +37,5 @@ AdjacentTriangle Adj() const; BamgVertex* EdgeVertex(const int & i) const; - Icoor2& det() const; + long long& det() const; }; } Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.cpp (revision 22085) @@ -5,28 +5,25 @@ BamgOpts::BamgOpts(){/*{{{*/ - this->anisomax=0; - this->cutoff=0; - this->coeff=0; - this->errg=0; - this->gradation=0; - this->Hessiantype=0; - this->MaxCornerAngle=0; - this->maxnbv=0; - this->maxsubdiv=0; - this->Metrictype=0; - this->nbjacobi=0; - this->nbsmooth=0; - this->omega=0; - this->power=0; - this->random=0; - this->verbose=0; + this->anisomax = 0; + this->cutoff = 0; + this->coeff = 0; + this->errg = 0; + this->gradation = 0; + this->Hessiantype = 0; + this->maxnbv = 0; + this->maxsubdiv = 0; + this->Metrictype = 0; + this->nbjacobi = 0; + this->nbsmooth = 0; + this->omega = 0; + this->power = 0; + this->verbose = 0; - this->Crack=0; - this->geometricalmetric=0; - this->KeepVertices=0; - this->splitcorners=0; + this->Crack = 0; + this->KeepVertices = 0; + this->splitcorners = 0; - this->hmin=0; - this->hmax=0; + this->hmin = 0; + this->hmax = 0; this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0; this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=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/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.h (revision 22085) @@ -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/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.cpp (revision 22085) @@ -122,5 +122,5 @@ } /*}}}*/ - double BamgVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){/*{{{*/ + double BamgVertex::Smoothing(Mesh &Th,Mesh &BTh,Triangle* &tstart ,double omega){/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/Smoothing)*/ @@ -157,5 +157,5 @@ double delta=Norme2_2(Xmove); - Icoor2 deta[3]; + long long deta[3]; I2 IBTh = BTh.R2ToI2(PNew); @@ -194,5 +194,4 @@ tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]); if (loop) { - BamgVertex *v0,*v1,*v2,*v3; if (tria->det<0) ok =1; else if ( (double)tria->det < detold/2 ) ok=1; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.h (revision 22085) @@ -42,5 +42,5 @@ /*methods (No constructor and no destructors...)*/ BamgVertex(); - double Smoothing(Mesh & ,const Mesh & ,Triangle * & ,double =1); + double Smoothing(Mesh & ,Mesh & ,Triangle * & ,double =1); void MetricFromHessian(const double Hxx,const double Hyx, const double Hyy, const double smin,const double smax,const double s,const double err,BamgOpts* bamgopts); void Echo(); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.cpp (revision 22085) @@ -63,75 +63,4 @@ return type &16; }/*}}}*/ - double GeomEdge::R1tg(double theta,R2 & t) const{/*{{{*/ - /*Original code from Frederic Hecht (BAMG v1.01, MeshGeom.cpp/R1tg)*/ - // 1/R of radius of cuvature - - R2 A=v[0]->r,B=v[1]->r; - double dca,dcb,dcta,dctb; - double ddca,ddcb,ddcta,ddctb; - double tt = theta*theta; - - //check theta - _assert_(theta>=0 && theta<=1); - - if (TgA()){ - if (TgB()){ - // Tangent A and B provided: - // interpolation d'hermite - dcb = 6*theta*(1-theta); - ddcb = 6*(1-2*theta); - dca = -dcb; - ddca = -ddcb; - dcta = (3*theta - 4)*theta + 1; - ddcta=6*theta-4; - dctb = 3*tt - 2*theta; - ddctb = 6*theta-2; - } - else { - //Tangent A provided but tangent B not provided - // 1-t*t, t-t*t, t*t - double t = theta; - dcb = 2*t; - ddcb = 2; - dca = -dcb; - ddca = -2; - dcta = 1-dcb; - ddcta = -ddcb; - dctb=0; - ddctb=0; - } - } - else{ - if (TgB()){ - //Tangent B provided but tangent A not provided - double t = 1-theta; - dca = -2*t; - ddca = 2; - dcb = -dca; - ddcb = -2; - dctb = 1+dca; - ddctb= ddca; - dcta =0; - ddcta =0; - } - else { - //Neither thangent A nor tangent B provided - // lagrange P1 - t=B-A; - return 0; - } - } - R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb; - R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb; - double d2=(d,d); - double sd2 = sqrt(d2); - t=d; - if(d2>1.0e-20){ - t/=sd2; - return Abs(Det(d,dd))/(d2*sd2); - } - else return 0; - } - /*}}}*/ int GeomEdge::Required() {/*{{{*/ return type &64; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.h (revision 22085) @@ -27,5 +27,4 @@ //Methods R2 F(double theta) const ; // parametrization of the curve edge - double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente int Cracked() const; int TgA() const; Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.cpp (revision 22085) @@ -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(); @@ -634,8 +621,4 @@ } } - if(ord != 2) { - vertices[i].SetCorner(); - } - vertices[i].SetCorner(); /*close the list around the vertex to have a circular loop*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.h (revision 22085) @@ -32,5 +32,4 @@ R2 pmin,pmax; // domain extrema coordinates double coefIcoor; // coef to integer coordinates; - double MaxCornerAngle; //Constructor/Destructors Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.cpp (revision 22085) @@ -188,10 +188,10 @@ } /*}}}*/ - void ListofIntersectionTriangles::SplitEdge(const Mesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/ + void ListofIntersectionTriangles::SplitEdge(Mesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/ListofIntersectionTriangles)*/ Triangle *tbegin, *t; - Icoor2 deta[3], deti,detj; + long long deta[3], deti,detj; double ba[3]; int ifirst=-1,ilast; @@ -199,5 +199,5 @@ int ocut,i,j,k=-1; // int OnAVertices =0; - Icoor2 dt[3]; + long long dt[3]; I2 a = Bh.R2ToI2(A) ,b= Bh.R2ToI2(B);// compute the Icoor a,b I2 vi,vj; @@ -356,5 +356,5 @@ k = OppositeVertex[ocut]; - Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b); + long long detbij = bamg::det((*t)[i],(*t)[j],b); if (detbij >= 0) { //we find the triangle contening b Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.h (revision 22085) @@ -63,5 +63,5 @@ int NewItem(Triangle *tt,double d0,double d1,double d2); int NewItem(R2 ,const Metric &); - void SplitEdge(const Mesh &,const R2 &,const R2 &,int nbegin=0); + void SplitEdge(Mesh &,const R2 &,const R2 &,int nbegin=0); double Length(); long NewPoints(BamgVertex *,long &nbv,long maxnbv); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp (revision 22085) @@ -133,5 +133,6 @@ nbt++; } - if (nbt==0 && nbv==0) { + if (nbt==0 && nbv==0){ + delete [] refv; _error_("All triangles have been removed"); } @@ -698,5 +699,5 @@ for (i=0;iSubDomains[i*4+0]=3; - bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]; + bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing bamgmesh->SubDomains[i*4+2]=1; bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber; @@ -980,55 +981,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 @@ -1049,5 +999,5 @@ } /*}}}*/ - void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3){/*{{{*/ + void Mesh::AddVertex( BamgVertex &s,Triangle* t, long long* det3){/*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/Add)*/ // ------------------------------- @@ -1068,5 +1018,5 @@ /*Intermediaries*/ Triangle* tt[3]; //the three triangles - Icoor2 det3local[3]; //three determinants (integer) + long long det3local[3]; //three determinants (integer) int nbzerodet =0; //number of zeros in det3 int izerodet=-1; //egde containing the vertex s @@ -1079,5 +1029,5 @@ //determinant of t - Icoor2 detOld=t->det; + long long detOld=t->det; /* infvertexindex = index of the infinite vertex (NULL) @@ -1231,10 +1181,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 +1191,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 +1390,5 @@ //check that nbsubdomains is empty - if (nbsubdomains){ - _error_("nbsubdomains should be 0"); - } + if(nbsubdomains) _error_("nbsubdomains should be 0"); nbsubdomains=0; @@ -1653,5 +1594,8 @@ long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]); k = edge4->SortAndFind(i0,i1); - _assert_(k>=0); + if(k<0){ + delete [] colorV; + _error_("This should not happen"); + } subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; subdomains[i].edge = edges+k; @@ -2113,5 +2057,13 @@ else if (xy==1) dd=dxdy; else if (xy==2) dd=dydy; - else _error_("not supported yet"); + else{ + delete [] detT; + delete [] Mmass; + delete [] workT; + delete [] workV; + delete [] Mmassxx; + delete [] OnBoundary; + _error_("not supported yet"); + } // do leat 2 iteration for boundary problem for (int ijacobi=0;ijacobi (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 +2693,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 @@ -2805,5 +2757,5 @@ //Find the triangle in which newvertex is located - Icoor2 det3[3]; + long long det3[3]; Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates) @@ -2825,5 +2777,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)*/ @@ -2832,5 +2784,5 @@ long i; long NbSwap=0; - Icoor2 det3[3]; + long long det3[3]; //number of new points @@ -2845,7 +2797,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; is3 (BAMG v1.01, Mesh2.cpp/MetricAt)*/ I2 a = R2ToI2(A); - Icoor2 deta[3]; + long long deta[3]; Triangle * t =TriangleFindFromCoord(a,deta); if (t->det <0) { // outside @@ -3086,7 +3039,7 @@ } } - if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); + //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); Bh.CreateSingleVertexToTriangleConnectivity(); - InsertNewPoints(nbvold,NbTSwap,bamgopts->random); + InsertNewPoints(nbvold,NbTSwap); } else Bh.CreateSingleVertexToTriangleConnectivity(); @@ -3146,5 +3099,5 @@ }// for triangle - if (!InsertNewPoints(nbvold,NbTSwap,bamgopts->random)) break; + if (!InsertNewPoints(nbvold,NbTSwap)) break; for (i=nbtold;ii,det3); quadtree->Add(*vi); @@ -3893,5 +3846,5 @@ long iv = nbvold; long NbSwap = 0; - Icoor2 det3[3]; + long long det3[3]; for (int i=nbvold;i (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/ @@ -3963,5 +3916,5 @@ } - Icoor2 detop ; + long long detop ; /*initialize number of test triangle*/ @@ -4009,5 +3962,5 @@ if (k==0) break; - if (k==2 && BinaryRand()) Exchange(ii[0],ii[1]); + if (k==2 && this->RandomNumber(1)) Exchange(ii[0],ii[1]); _assert_(k<3); AdjacentTriangle t1 = t->Adj(jj=ii[0]); @@ -4095,5 +4048,5 @@ /*Insert Vertices*/ - Insert(true); + Insert(); } /*}}}*/ @@ -4396,5 +4349,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 @@ -4460,4 +4413,5 @@ printf("\n"); if(NbVerticesOnGeomVertex >= maxnbv){ + delete [] bcurve; _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv); } @@ -4722,5 +4676,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 +4689,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 (this->RandomNumber(1)) + 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; } /*}}}*/ @@ -4748,5 +4871,5 @@ } int kkk=0; - Icoor2 IJ_IA,IJ_AJ; + long long IJ_IA,IJ_AJ; AdjacentTriangle edge(t,OppositeEdge[k]); for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { @@ -4778,80 +4901,5 @@ } /*}}}*/ - 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 - - Icoor2 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"); - } - Icoor2 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; - } - /*}}}*/ - void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){ /*{{{*/ + void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,long long det1,long long det2){ /*{{{*/ /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/swap)*/ // -------------------------------------------------------------- @@ -4897,90 +4945,3 @@ } // end swap /*}}}*/ - int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & 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]]; - - Icoor2 dets2 = det(*pva,*pvb,s2); - Icoor2 det1=t1->det , det2=t2->det ; - Icoor2 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)"); - } - Icoor2 ndet1 = bamg::det(s1,sa,s2); - Icoor2 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/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.h (revision 22085) @@ -38,6 +38,7 @@ R2 pmin,pmax; // extrema - double coefIcoor; // coef to integer Icoor1; + double coefIcoor; // coef to integer ListofIntersectionTriangles lIntTria; + int randomseed; //used for random number generation long NbVerticesOnGeomVertex; @@ -76,6 +77,6 @@ I2 R2ToI2(const R2 & P) const; R2 I2ToR2(const I2 & P) const; - void AddVertex(BamgVertex & s,Triangle * t,Icoor2 * =0) ; - void Insert(bool random); + void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; + void Insert(); void Echo(void); void ForceBoundary(); @@ -91,8 +92,8 @@ 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); - Metric MetricAt (const R2 &) const; + Metric MetricAt (const R2 &); GeomEdge* ProjectOnCurve( Edge & AB, BamgVertex & A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR); long GetId(const Triangle & t) const; @@ -103,5 +104,5 @@ long GetId(const Edge * t) const; BamgVertex* NearestVertex(int i,int j) ; - Triangle* TriangleFindFromCoord(const I2 & ,Icoor2 [3],Triangle *tstart=0) const; + Triangle* TriangleFindFromCoord(const I2 & ,long long [3],Triangle *tstart=0); void ReadMesh(int* index,double* x,double* y,int nods,int nels); void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); @@ -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); }; @@ -149,9 +154,6 @@ void swap(Triangle *t1,short a1, Triangle *t2,short a2, - BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2); - int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , - AdjacentTriangle & tt1,Icoor2 & dets1, - Icoor2 & detsa,Icoor2 & detsb, int & nbswap); - int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; + BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); + inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); @@ -166,21 +168,4 @@ int j=i;i=on->AdjVertexIndex[i];on=on->Adj[j]; } - inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){ - double ret; - I2 ia=va,ib=vb,ic=vc; - I2 ab=ib-ia,bc=ic-ib,ac=ic-ia; - Icoor2 deta=Det(ab,ac); - if (deta <=0) ret = -1; - else { - double a = sqrt((double) (ac,ac)), - b = sqrt((double) (bc,bc)), - c = sqrt((double) (ab,ab)), - p = a+b+c; - double h= Max(Max(a,b),c),ro=deta/p; - ret = ro/h; - } - return ret; - } - } #endif Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.cpp (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.cpp (revision 22085) @@ -238,8 +238,8 @@ BamgVertex *s2=t2->vertices[OppositeVertex[a2]]; - Icoor2 det1=t1->det , det2=t2->det ; - Icoor2 detT = det1+det2; - Icoor2 detA = Abs(det1) + Abs(det2); - Icoor2 detMin = Min(det1,det2); + long long det1=t1->det , det2=t2->det ; + long long detT = det1+det2; + long long detA = Abs(det1) + Abs(det2); + long long detMin = Min(det1,det2); int OnSwap = 0; @@ -256,6 +256,5 @@ OnSwap = (Abs(det1) + Abs(det2)) < detA; - Icoor2 detMinNew=Min(det1,det2); - // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test + long long detMinNew=Min(det1,det2); if (! OnSwap &&(detMinNew>0)) { OnSwap = detMin ==0; @@ -265,5 +264,5 @@ if(kopt) { // critere de Delaunay pure isotrope - Icoor2 xb1 = sb->i.x - s1->i.x, + long long xb1 = sb->i.x - s1->i.x, x21 = s2->i.x - s1->i.x, yb1 = sb->i.y - s1->i.y, Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.h (revision 22085) @@ -21,5 +21,5 @@ public: - Icoor2 det; //Integer determinant (twice its area) + long long det; //Integer determinant (twice its area) union { Triangle *link; @@ -60,5 +60,4 @@ //Inline methods - double qualite() ; void Set(const Triangle &,const Mesh &,Mesh &); int In(BamgVertex *v) const { return vertices[0]==v || vertices[1]==v || vertices[2]==v ;} Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/det.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/det.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/det.h (revision 22085) @@ -6,7 +6,7 @@ namespace bamg { - Icoor2 inline det(const I2 &a,const I2 & b,const I2 &c){ - Icoor2 bax = b.x - a.x ,bay = b.y - a.y; - Icoor2 cax = c.x - a.x ,cay = c.y - a.y; + long long inline det(const I2 &a,const I2 & b,const I2 &c){ + long long bax = b.x - a.x ,bay = b.y - a.y; + long long cax = c.x - a.x ,cay = c.y - a.y; return bax*cay - bay*cax; } Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/typedefs.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/typedefs.h (revision 22084) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/typedefs.h (revision 22085) @@ -7,9 +7,7 @@ /*Integer coordinates types*/ - typedef int Icoor1; - typedef long long Icoor2; /*I2 and R2*/ - typedef P2 I2; + typedef P2 I2; typedef P2 R2; }