![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22085
- Timestamp:
- 09/15/17 18:44:44 (8 years ago)
- Location:
- issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg
- Files:
-
- 18 edited
-
AdjacentTriangle.cpp (modified) (1 diff)
-
AdjacentTriangle.h (modified) (1 diff)
-
BamgOpts.cpp (modified) (2 diffs)
-
BamgOpts.h (modified) (2 diffs)
-
BamgVertex.cpp (modified) (3 diffs)
-
BamgVertex.h (modified) (1 diff)
-
GeomEdge.cpp (modified) (1 diff)
-
GeomEdge.h (modified) (1 diff)
-
Geometry.cpp (modified) (5 diffs)
-
Geometry.h (modified) (1 diff)
-
ListofIntersectionTriangles.cpp (modified) (3 diffs)
-
ListofIntersectionTriangles.h (modified) (1 diff)
-
Mesh.cpp (modified) (36 diffs)
-
Mesh.h (modified) (8 diffs)
-
Triangle.cpp (modified) (3 diffs)
-
Triangle.h (modified) (2 diffs)
-
det.h (modified) (1 diff)
-
typedefs.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.cpp ¶
r18064 r22085 34 34 } 35 35 /*}}}*/ 36 Icoor2& AdjacentTriangle::det() const {/*{{{*/36 long long & AdjacentTriangle::det() const {/*{{{*/ 37 37 return t->det; 38 38 } -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/AdjacentTriangle.h ¶
r16233 r22085 37 37 AdjacentTriangle Adj() const; 38 38 BamgVertex* EdgeVertex(const int & i) const; 39 Icoor2& det() const;39 long long& det() const; 40 40 }; 41 41 } -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.cpp ¶
r18064 r22085 5 5 BamgOpts::BamgOpts(){/*{{{*/ 6 6 7 this->anisomax=0; 8 this->cutoff=0; 9 this->coeff=0; 10 this->errg=0; 11 this->gradation=0; 12 this->Hessiantype=0; 13 this->MaxCornerAngle=0; 14 this->maxnbv=0; 15 this->maxsubdiv=0; 16 this->Metrictype=0; 17 this->nbjacobi=0; 18 this->nbsmooth=0; 19 this->omega=0; 20 this->power=0; 21 this->random=0; 22 this->verbose=0; 7 this->anisomax = 0; 8 this->cutoff = 0; 9 this->coeff = 0; 10 this->errg = 0; 11 this->gradation = 0; 12 this->Hessiantype = 0; 13 this->maxnbv = 0; 14 this->maxsubdiv = 0; 15 this->Metrictype = 0; 16 this->nbjacobi = 0; 17 this->nbsmooth = 0; 18 this->omega = 0; 19 this->power = 0; 20 this->verbose = 0; 23 21 24 this->Crack=0; 25 this->geometricalmetric=0; 26 this->KeepVertices=0; 27 this->splitcorners=0; 22 this->Crack = 0; 23 this->KeepVertices = 0; 24 this->splitcorners = 0; 28 25 29 this->hmin =0;30 this->hmax =0;26 this->hmin = 0; 27 this->hmax = 0; 31 28 this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0; 32 29 this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=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"); -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgOpts.h ¶
r16967 r22085 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; -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.cpp ¶
r21759 r22085 122 122 } 123 123 /*}}}*/ 124 double BamgVertex::Smoothing(Mesh &Th, constMesh &BTh,Triangle* &tstart ,double omega){/*{{{*/124 double BamgVertex::Smoothing(Mesh &Th,Mesh &BTh,Triangle* &tstart ,double omega){/*{{{*/ 125 125 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/ 126 126 … … 157 157 double delta=Norme2_2(Xmove); 158 158 159 Icoor2deta[3];159 long long deta[3]; 160 160 I2 IBTh = BTh.R2ToI2(PNew); 161 161 … … 194 194 tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]); 195 195 if (loop) { 196 BamgVertex *v0,*v1,*v2,*v3;197 196 if (tria->det<0) ok =1; 198 197 else if ( (double)tria->det < detold/2 ) ok=1; -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/BamgVertex.h ¶
r21759 r22085 42 42 /*methods (No constructor and no destructors...)*/ 43 43 BamgVertex(); 44 double Smoothing(Mesh & , constMesh & ,Triangle * & ,double =1);44 double Smoothing(Mesh & ,Mesh & ,Triangle * & ,double =1); 45 45 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); 46 46 void Echo(); -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.cpp ¶
r18064 r22085 63 63 return type &16; 64 64 }/*}}}*/ 65 double GeomEdge::R1tg(double theta,R2 & t) const{/*{{{*/66 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/67 // 1/R of radius of cuvature68 69 R2 A=v[0]->r,B=v[1]->r;70 double dca,dcb,dcta,dctb;71 double ddca,ddcb,ddcta,ddctb;72 double tt = theta*theta;73 74 //check theta75 _assert_(theta>=0 && theta<=1);76 77 if (TgA()){78 if (TgB()){79 // Tangent A and B provided:80 // interpolation d'hermite81 dcb = 6*theta*(1-theta);82 ddcb = 6*(1-2*theta);83 dca = -dcb;84 ddca = -ddcb;85 dcta = (3*theta - 4)*theta + 1;86 ddcta=6*theta-4;87 dctb = 3*tt - 2*theta;88 ddctb = 6*theta-2;89 }90 else {91 //Tangent A provided but tangent B not provided92 // 1-t*t, t-t*t, t*t93 double t = theta;94 dcb = 2*t;95 ddcb = 2;96 dca = -dcb;97 ddca = -2;98 dcta = 1-dcb;99 ddcta = -ddcb;100 dctb=0;101 ddctb=0;102 }103 }104 else{105 if (TgB()){106 //Tangent B provided but tangent A not provided107 double t = 1-theta;108 dca = -2*t;109 ddca = 2;110 dcb = -dca;111 ddcb = -2;112 dctb = 1+dca;113 ddctb= ddca;114 dcta =0;115 ddcta =0;116 }117 else {118 //Neither thangent A nor tangent B provided119 // lagrange P1120 t=B-A;121 return 0;122 }123 }124 R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb;125 R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;126 double d2=(d,d);127 double sd2 = sqrt(d2);128 t=d;129 if(d2>1.0e-20){130 t/=sd2;131 return Abs(Det(d,dd))/(d2*sd2);132 }133 else return 0;134 }135 /*}}}*/136 65 int GeomEdge::Required() {/*{{{*/ 137 66 return type &64; -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/GeomEdge.h ¶
r16231 r22085 27 27 //Methods 28 28 R2 F(double theta) const ; // parametrization of the curve edge 29 double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente30 29 int Cracked() const; 31 30 int TgA() const; -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.cpp ¶
r21759 r22085 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(); … … 634 621 } 635 622 } 636 if(ord != 2) {637 vertices[i].SetCorner();638 }639 vertices[i].SetCorner();640 623 641 624 /*close the list around the vertex to have a circular loop*/ -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Geometry.h ¶
r21759 r22085 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 -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.cpp ¶
r18064 r22085 188 188 } 189 189 /*}}}*/ 190 void ListofIntersectionTriangles::SplitEdge( constMesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/190 void ListofIntersectionTriangles::SplitEdge(Mesh & Bh, const R2 &A,const R2 &B,int nbegin) {/*{{{*/ 191 191 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ListofIntersectionTriangles)*/ 192 192 193 193 Triangle *tbegin, *t; 194 194 195 Icoor2deta[3], deti,detj;195 long long deta[3], deti,detj; 196 196 double ba[3]; 197 197 int ifirst=-1,ilast; … … 199 199 int ocut,i,j,k=-1; 200 200 // int OnAVertices =0; 201 Icoor2dt[3];201 long long dt[3]; 202 202 I2 a = Bh.R2ToI2(A) ,b= Bh.R2ToI2(B);// compute the Icoor a,b 203 203 I2 vi,vj; … … 356 356 k = OppositeVertex[ocut]; 357 357 358 Icoor2detbij = bamg::det((*t)[i],(*t)[j],b);358 long long detbij = bamg::det((*t)[i],(*t)[j],b); 359 359 360 360 if (detbij >= 0) { //we find the triangle contening b -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/ListofIntersectionTriangles.h ¶
r16233 r22085 63 63 int NewItem(Triangle *tt,double d0,double d1,double d2); 64 64 int NewItem(R2 ,const Metric &); 65 void SplitEdge( constMesh &,const R2 &,const R2 &,int nbegin=0);65 void SplitEdge(Mesh &,const R2 &,const R2 &,int nbegin=0); 66 66 double Length(); 67 67 long NewPoints(BamgVertex *,long &nbv,long maxnbv); -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.cpp ¶
r21759 r22085 133 133 nbt++; 134 134 } 135 if (nbt==0 && nbv==0) { 135 if (nbt==0 && nbv==0){ 136 delete [] refv; 136 137 _error_("All triangles have been removed"); 137 138 } … … 698 699 for (i=0;i<nbsubdomains;i++){ 699 700 bamgmesh->SubDomains[i*4+0]=3; 700 bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)] ;701 bamgmesh->SubDomains[i*4+1]=reft[GetId(subdomains[i].head)]+1;//MATLAB indexing 701 702 bamgmesh->SubDomains[i*4+2]=1; 702 703 bamgmesh->SubDomains[i*4+3]=subdomains[i].ReferenceNumber; … … 980 981 981 982 /*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 983 void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ 1034 984 // Hessiantype = 0 => H is computed using double L2 projection … … 1049 999 } 1050 1000 /*}}}*/ 1051 void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3){/*{{{*/1001 void Mesh::AddVertex( BamgVertex &s,Triangle* t, long long* det3){/*{{{*/ 1052 1002 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Add)*/ 1053 1003 // ------------------------------- … … 1068 1018 /*Intermediaries*/ 1069 1019 Triangle* tt[3]; //the three triangles 1070 Icoor2det3local[3]; //three determinants (integer)1020 long long det3local[3]; //three determinants (integer) 1071 1021 int nbzerodet =0; //number of zeros in det3 1072 1022 int izerodet=-1; //egde containing the vertex s … … 1079 1029 1080 1030 //determinant of t 1081 Icoor2detOld=t->det;1031 long long detOld=t->det; 1082 1032 1083 1033 /* infvertexindex = index of the infinite vertex (NULL) … … 1231 1181 int i,j,k,kk,it,jt; 1232 1182 int verbose=0; 1233 double cutoffradian=10*Pi/180;1234 1183 1235 1184 /*Recover options*/ 1236 if 1185 if(bamgopts){ 1237 1186 verbose=bamgopts->verbose; 1238 cutoffradian=bamgopts->MaxCornerAngle*Pi/180;1239 1187 } 1240 1188 … … 1243 1191 1244 1192 //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; 1193 if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); 1251 1194 1252 1195 /*Construction of the edges*/ … … 1447 1390 1448 1391 //check that nbsubdomains is empty 1449 if (nbsubdomains){ 1450 _error_("nbsubdomains should be 0"); 1451 } 1392 if(nbsubdomains) _error_("nbsubdomains should be 0"); 1452 1393 nbsubdomains=0; 1453 1394 … … 1653 1594 long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]); 1654 1595 k = edge4->SortAndFind(i0,i1); 1655 _assert_(k>=0); 1596 if(k<0){ 1597 delete [] colorV; 1598 _error_("This should not happen"); 1599 } 1656 1600 subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1; 1657 1601 subdomains[i].edge = edges+k; … … 2113 2057 else if (xy==1) dd=dxdy; 2114 2058 else if (xy==2) dd=dydy; 2115 else _error_("not supported yet"); 2059 else{ 2060 delete [] detT; 2061 delete [] Mmass; 2062 delete [] workT; 2063 delete [] workV; 2064 delete [] Mmassxx; 2065 delete [] OnBoundary; 2066 _error_("not supported yet"); 2067 } 2116 2068 // do leat 2 iteration for boundary problem 2117 2069 for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){ … … 2162 2114 delete [] dxdy; 2163 2115 delete [] dydy; 2164 delete [] 2116 delete [] workT; 2165 2117 delete [] workV; 2166 2118 delete [] Mmassxx; 2167 delete [] 2119 delete [] OnBoundary; 2168 2120 2169 2121 } … … 2453 2405 if (nbt == nbtout || !NbSubDomTot) { 2454 2406 delete [] HeapArete; 2407 delete [] HeapTriangle; 2455 2408 _error_("The boundary is not close: all triangles are outside"); 2456 2409 } … … 2649 2602 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2650 2603 2651 /* initialize random seed: */2652 srand(19999999);2653 2654 2604 /*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; 2605 this->NbRef = 0; 2606 this->quadtree = NULL; 2607 this->nbv = 0; 2608 this->nbt = 0; 2609 this->nbe = 0; 2610 this->edges = NULL; 2611 this->nbsubdomains = 0; 2612 this->subdomains = NULL; 2613 this->maxnbv = maxnbv_in; 2614 this->maxnbt = 2 *maxnbv_in-2; 2615 this->NbVertexOnBThVertex = 0; 2616 this->VertexOnBThVertex = NULL; 2617 this->NbVertexOnBThEdge = 0; 2618 this->VertexOnBThEdge = NULL; 2619 this->NbCrackedVertices = 0; 2620 this->CrackedVertices = NULL; 2621 this->NbCrackedEdges = 0; 2622 this->CrackedEdges = NULL; 2623 this->NbVerticesOnGeomVertex = 0; 2624 this->VerticesOnGeomVertex = NULL; 2625 this->NbVerticesOnGeomEdge = 0; 2626 this->VerticesOnGeomEdge = NULL; 2627 2628 /*Initialize random seed*/ 2629 this->randomseed = 1; 2677 2630 2678 2631 /*Allocate if maxnbv_in>0*/ 2679 2632 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);2633 this->vertices=new BamgVertex[this->maxnbv]; 2634 this->orderedvertices=new BamgVertex* [this->maxnbv]; 2635 this->triangles=new Triangle[this->maxnbt]; 2636 _assert_(this->vertices); 2637 _assert_(this->orderedvertices); 2638 _assert_(this->triangles); 2686 2639 } 2687 2640 else{ 2688 vertices=NULL;2689 orderedvertices=NULL;2690 t riangles=NULL;2691 maxnbt=0;2641 this->vertices = NULL; 2642 this->orderedvertices = NULL; 2643 this->triangles = NULL; 2644 this->maxnbt = 0; 2692 2645 } 2693 2646 } 2694 2647 /*}}}*/ 2695 void Mesh::Insert( bool random){/*{{{*/2648 void Mesh::Insert(){/*{{{*/ 2696 2649 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2697 2650 … … 2740 2693 //Get Prime number 2741 2694 const long PrimeNumber= BigPrimeNumber(nbv); 2742 int temp = rand(); if(!random) temp = 756804110; 2743 int k0=temp%nbv; 2695 int k0=this->RandomNumber(nbv); 2744 2696 2745 2697 //Build orderedvertices … … 2805 2757 2806 2758 //Find the triangle in which newvertex is located 2807 Icoor2det3[3];2759 long long det3[3]; 2808 2760 Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates) 2809 2761 … … 2825 2777 } 2826 2778 /*}}}*/ 2827 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap ,bool random) {/*{{{*/2779 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ 2828 2780 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2829 2781 … … 2832 2784 long i; 2833 2785 long NbSwap=0; 2834 Icoor2det3[3];2786 long long det3[3]; 2835 2787 2836 2788 //number of new points … … 2845 2797 /*construction of a random order*/ 2846 2798 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; 2799 long k3 = long(this->RandomNumber(nbvnew)); 2850 2800 //loop over the new points 2851 2801 for (int is3=0; is3<nbvnew; is3++){ … … 2945 2895 } 2946 2896 } 2947 if(kk) _error_("See above"); 2897 if(kk){ 2898 delete [] e; 2899 _error_("See above"); 2900 } 2948 2901 2949 2902 return e; … … 3013 2966 } 3014 2967 /*}}}*/ 3015 Metric Mesh::MetricAt(const R2 & A) const{ /*{{{*/2968 Metric Mesh::MetricAt(const R2 & A){ /*{{{*/ 3016 2969 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/MetricAt)*/ 3017 2970 3018 2971 I2 a = R2ToI2(A); 3019 Icoor2deta[3];2972 long long deta[3]; 3020 2973 Triangle * t =TriangleFindFromCoord(a,deta); 3021 2974 if (t->det <0) { // outside … … 3086 3039 } 3087 3040 } 3088 if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n");3041 //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); 3089 3042 Bh.CreateSingleVertexToTriangleConnectivity(); 3090 InsertNewPoints(nbvold,NbTSwap ,bamgopts->random);3043 InsertNewPoints(nbvold,NbTSwap); 3091 3044 } 3092 3045 else Bh.CreateSingleVertexToTriangleConnectivity(); … … 3146 3099 }// for triangle 3147 3100 3148 if (!InsertNewPoints(nbvold,NbTSwap ,bamgopts->random)) break;3101 if (!InsertNewPoints(nbvold,NbTSwap)) break; 3149 3102 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3150 3103 Headt = nbt; // empty list … … 3476 3429 for (int icount=2; icount<nbvb; icount++) { 3477 3430 BamgVertex *vi = orderedvertices[icount]; 3478 Icoor2det3[3];3431 long long det3[3]; 3479 3432 Triangle *tcvi = TriangleFindFromCoord(vi->i,det3); 3480 3433 quadtree->Add(*vi); … … 3893 3846 long iv = nbvold; 3894 3847 long NbSwap = 0; 3895 Icoor2det3[3];3848 long long det3[3]; 3896 3849 for (int i=nbvold;i<nbv;i++) {// for all the new point 3897 3850 BamgVertex & vi = vertices[i]; … … 3934 3887 } 3935 3888 /*}}}*/ 3936 Triangle * Mesh::TriangleFindFromCoord(const I2 & B, Icoor2 det3[3], Triangle *tstart) const{/*{{{*/3889 Triangle * Mesh::TriangleFindFromCoord(const I2 & B,long long det3[3], Triangle *tstart){/*{{{*/ 3937 3890 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/ 3938 3891 … … 3963 3916 } 3964 3917 3965 Icoor2detop ;3918 long long detop ; 3966 3919 3967 3920 /*initialize number of test triangle*/ … … 4009 3962 4010 3963 if (k==0) break; 4011 if (k==2 && BinaryRand()) Exchange(ii[0],ii[1]);3964 if (k==2 && this->RandomNumber(1)) Exchange(ii[0],ii[1]); 4012 3965 _assert_(k<3); 4013 3966 AdjacentTriangle t1 = t->Adj(jj=ii[0]); … … 4095 4048 4096 4049 /*Insert Vertices*/ 4097 Insert( true);4050 Insert(); 4098 4051 } 4099 4052 /*}}}*/ … … 4396 4349 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4397 4350 if (verbose>3) _printf_(" Inserting boundary points\n"); 4398 Insert( bamgopts->random);4351 Insert(); 4399 4352 4400 4353 //Force the boundary … … 4460 4413 printf("\n"); 4461 4414 if(NbVerticesOnGeomVertex >= maxnbv){ 4415 delete [] bcurve; 4462 4416 _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv); 4463 4417 } … … 4722 4676 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4723 4677 if (verbose>3) _printf_(" Inserting boundary points\n"); 4724 Insert( bamgopts->random);4678 Insert(); 4725 4679 4726 4680 //Force the boundary … … 4735 4689 NewPoints(BTh,bamgopts,KeepVertices) ; 4736 4690 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); 4691 } 4692 /*}}}*/ 4693 int Mesh::RandomNumber(int max){/*{{{*/ 4694 /* Generate a random number from 0 to max-1 using linear congruential 4695 * random number generator*/ 4696 4697 this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l; 4698 return this->randomseed / (714025l / max + 1); 4699 } /*}}}*/ 4700 int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4701 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4702 4703 int NbSwap =0; 4704 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4705 _error_("!a.t || !b.t"); 4706 } 4707 int k=0; 4708 taret=AdjacentTriangle(0,0); // erreur 4709 4710 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4711 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4712 // we turn around a in the direct direction 4713 4714 long long det2 = v2 ? det(*v2,a,b): -1 , det1; 4715 if(v2) // normal case 4716 det2 = det(*v2,a,b); 4717 else { // no chance infini vertex try the next 4718 tta= Previous(Adj(tta)); 4719 v2 = tta.EdgeVertex(0); 4720 vbegin =v2; 4721 if (!v2){ 4722 _error_("!v2"); 4723 } 4724 det2 = det(*v2,a,b); 4725 } 4726 4727 while (v2 != &b) { 4728 AdjacentTriangle tc = Previous(Adj(tta)); 4729 v1 = v2; 4730 v2 = tc.EdgeVertex(0); 4731 det1 = det2; 4732 det2 = v2 ? det(*v2,a,b): det2; 4733 4734 if((det1 < 0) && (det2 >0)) { 4735 // try to force the edge 4736 BamgVertex * va = &a, *vb = &b; 4737 tc = Previous(tc); 4738 if (!v1 || !v2){ 4739 _error_("!v1 || !v2"); 4740 } 4741 long long detss = 0,l=0; 4742 while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4743 if(l++ > 10000000) { 4744 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4745 } 4746 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4747 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4748 tc.SetLock(); 4749 a.Optim(1,0); 4750 b.Optim(1,0); 4751 taret = tc; 4752 return NbSwap; 4753 } 4754 else 4755 { 4756 taret = tc; 4757 return -2; // error boundary is crossing 4758 } 4759 } 4760 tta = tc; 4761 k++; 4762 if (k>=2000){ 4763 _error_("k>=2000"); 4764 } 4765 if ( vbegin == v2 ) return -1;// error 4766 } 4767 4768 tta.SetLock(); 4769 taret=tta; 4770 a.Optim(1,0); 4771 b.Optim(1,0); 4772 return NbSwap; 4773 } 4774 /*}}}*/ 4775 int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ 4776 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ 4777 // l'arete ta coupe l'arete pva pvb 4778 // de cas apres le swap sa coupe toujours 4779 // on cherche l'arete suivante 4780 // on suppose que detsa >0 et detsb <0 4781 // attention la routine echange pva et pvb 4782 4783 if(tt1.Locked()) return 0; // frontiere croise 4784 4785 AdjacentTriangle tt2 = Adj(tt1); 4786 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 4787 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 4788 if ( a1<0 || a1>=3 ){ 4789 _error_("a1<0 || a1>=3"); 4790 } 4791 4792 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 4793 BamgVertex & s1= (*t1)[OppositeVertex[a1]]; 4794 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 4795 4796 long long dets2 = det(*pva,*pvb,s2); 4797 long long det1=t1->det , det2=t2->det ; 4798 long long detT = det1+det2; 4799 if ((det1<=0 ) || (det2<=0)){ 4800 _error_("(det1<=0 ) || (det2<=0)"); 4801 } 4802 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 4803 _error_("(detsa>=0) || (detsb<=0)"); 4804 } 4805 long long ndet1 = bamg::det(s1,sa,s2); 4806 long long ndet2 = detT - ndet1; 4807 4808 int ToSwap =0; //pas de swap 4809 if ((ndet1 >0) && (ndet2 >0)) 4810 { // on peut swaper 4811 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 4812 ToSwap =1; 4813 else // swap alleatoire 4814 if (this->RandomNumber(1)) 4815 ToSwap =2; 4816 } 4817 if (ToSwap) NbSwap++, 4818 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 4819 4820 int ret=1; 4821 4822 if (dets2 < 0) {// haut 4823 dets1 = ToSwap ? dets1 : detsa ; 4824 detsa = dets2; 4825 tt1 = Previous(tt2) ;} 4826 else if (dets2 > 0){// bas 4827 dets1 = ToSwap ? dets1 : detsb ; 4828 detsb = dets2; 4829 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 4830 if(!ToSwap) tt1 = Next(tt2); 4831 } 4832 else { // changement de direction 4833 ret = -1; 4834 Exchange(pva,pvb); 4835 Exchange(detsa,detsb); 4836 Exchange(dets1,dets2); 4837 Exchange(tt1,tt2); 4838 dets1=-dets1; 4839 dets2=-dets2; 4840 detsa=-detsa; 4841 detsb=-detsb; 4842 4843 if(ToSwap){ 4844 if (dets2 < 0) {// haut 4845 dets1 = (ToSwap ? dets1 : detsa) ; 4846 detsa = dets2; 4847 tt1 = Previous(tt2) ;} 4848 else if(dets2 > 0){// bas 4849 dets1 = (ToSwap ? dets1 : detsb) ; 4850 detsb = dets2; 4851 if(!ToSwap) tt1 = Next(tt2); 4852 } 4853 else {// on a fin ??? 4854 tt1 = Next(tt2); 4855 ret =0;} 4856 } 4857 4858 } 4859 return ret; 4737 4860 } 4738 4861 /*}}}*/ … … 4748 4871 } 4749 4872 int kkk=0; 4750 Icoor2IJ_IA,IJ_AJ;4873 long long IJ_IA,IJ_AJ; 4751 4874 AdjacentTriangle edge(t,OppositeEdge[k]); 4752 4875 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { … … 4778 4901 } 4779 4902 /*}}}*/ 4780 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4781 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4782 4783 int NbSwap =0; 4784 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4785 _error_("!a.t || !b.t"); 4786 } 4787 int k=0; 4788 taret=AdjacentTriangle(0,0); // erreur 4789 4790 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4791 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4792 // we turn around a in the direct direction 4793 4794 Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; 4795 if(v2) // normal case 4796 det2 = det(*v2,a,b); 4797 else { // no chance infini vertex try the next 4798 tta= Previous(Adj(tta)); 4799 v2 = tta.EdgeVertex(0); 4800 vbegin =v2; 4801 if (!v2){ 4802 _error_("!v2"); 4803 } 4804 det2 = det(*v2,a,b); 4805 } 4806 4807 while (v2 != &b) { 4808 AdjacentTriangle tc = Previous(Adj(tta)); 4809 v1 = v2; 4810 v2 = tc.EdgeVertex(0); 4811 det1 = det2; 4812 det2 = v2 ? det(*v2,a,b): det2; 4813 4814 if((det1 < 0) && (det2 >0)) { 4815 // try to force the edge 4816 BamgVertex * va = &a, *vb = &b; 4817 tc = Previous(tc); 4818 if (!v1 || !v2){ 4819 _error_("!v1 || !v2"); 4820 } 4821 Icoor2 detss = 0,l=0; 4822 while ((SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4823 if(l++ > 10000000) { 4824 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4825 } 4826 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4827 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4828 tc.SetLock(); 4829 a.Optim(1,0); 4830 b.Optim(1,0); 4831 taret = tc; 4832 return NbSwap; 4833 } 4834 else 4835 { 4836 taret = tc; 4837 return -2; // error boundary is crossing 4838 } 4839 } 4840 tta = tc; 4841 k++; 4842 if (k>=2000){ 4843 _error_("k>=2000"); 4844 } 4845 if ( vbegin == v2 ) return -1;// error 4846 } 4847 4848 tta.SetLock(); 4849 taret=tta; 4850 a.Optim(1,0); 4851 b.Optim(1,0); 4852 return NbSwap; 4853 } 4854 /*}}}*/ 4855 void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2){ /*{{{*/ 4903 void swap(Triangle *t1,short a1, Triangle *t2,short a2, BamgVertex *s1,BamgVertex *s2,long long det1,long long det2){ /*{{{*/ 4856 4904 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/swap)*/ 4857 4905 // -------------------------------------------------------------- … … 4897 4945 } // end swap 4898 4946 /*}}}*/ 4899 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & 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 Icoor2 dets2 = det(*pva,*pvb,s2);4921 Icoor2 det1=t1->det , det2=t2->det ;4922 Icoor2 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 Icoor2 ndet1 = bamg::det(s1,sa,s2);4930 Icoor2 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 4947 } -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Mesh.h ¶
r21759 r22085 38 38 39 39 R2 pmin,pmax; // extrema 40 double coefIcoor; // coef to integer Icoor1;40 double coefIcoor; // coef to integer 41 41 ListofIntersectionTriangles lIntTria; 42 int randomseed; //used for random number generation 42 43 43 44 long NbVerticesOnGeomVertex; … … 76 77 I2 R2ToI2(const R2 & P) const; 77 78 R2 I2ToR2(const I2 & P) const; 78 void AddVertex(BamgVertex & s,Triangle * t, Icoor2* =0) ;79 void Insert( bool random);79 void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; 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); 96 Metric MetricAt (const R2 &) const;97 Metric MetricAt (const R2 &); 97 98 GeomEdge* ProjectOnCurve( Edge & AB, BamgVertex & A, BamgVertex & B,double theta, BamgVertex & R,VertexOnEdge & BR,VertexOnGeom & GR); 98 99 long GetId(const Triangle & t) const; … … 103 104 long GetId(const Edge * t) const; 104 105 BamgVertex* NearestVertex(int i,int j) ; 105 Triangle* TriangleFindFromCoord(const I2 & , Icoor2 [3],Triangle *tstart=0) const;106 Triangle* TriangleFindFromCoord(const I2 & ,long long [3],Triangle *tstart=0); 106 107 void ReadMesh(int* index,double* x,double* y,int nods,int nels); 107 108 void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts); … … 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 … … 149 154 void swap(Triangle *t1,short a1, 150 155 Triangle *t2,short a2, 151 BamgVertex *s1,BamgVertex *s2,Icoor2 det1,Icoor2 det2); 152 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , 153 AdjacentTriangle & tt1,Icoor2 & dets1, 154 Icoor2 & detsa,Icoor2 & detsb, int & nbswap); 155 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; 156 BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); 157 156 158 inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ 157 159 return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); … … 166 168 int j=i;i=on->AdjVertexIndex[i];on=on->Adj[j]; 167 169 } 168 inline double qualite(const BamgVertex &va,const BamgVertex &vb,const BamgVertex &vc){169 double ret;170 I2 ia=va,ib=vb,ic=vc;171 I2 ab=ib-ia,bc=ic-ib,ac=ic-ia;172 Icoor2 deta=Det(ab,ac);173 if (deta <=0) ret = -1;174 else {175 double a = sqrt((double) (ac,ac)),176 b = sqrt((double) (bc,bc)),177 c = sqrt((double) (ab,ab)),178 p = a+b+c;179 double h= Max(Max(a,b),c),ro=deta/p;180 ret = ro/h;181 }182 return ret;183 }184 185 170 } 186 171 #endif -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.cpp ¶
r21759 r22085 238 238 BamgVertex *s2=t2->vertices[OppositeVertex[a2]]; 239 239 240 Icoor2det1=t1->det , det2=t2->det ;241 Icoor2detT = det1+det2;242 Icoor2detA = Abs(det1) + Abs(det2);243 Icoor2detMin = Min(det1,det2);240 long long det1=t1->det , det2=t2->det ; 241 long long detT = det1+det2; 242 long long detA = Abs(det1) + Abs(det2); 243 long long detMin = Min(det1,det2); 244 244 245 245 int OnSwap = 0; … … 256 256 OnSwap = (Abs(det1) + Abs(det2)) < detA; 257 257 258 Icoor2 detMinNew=Min(det1,det2); 259 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test 258 long long detMinNew=Min(det1,det2); 260 259 if (! OnSwap &&(detMinNew>0)) { 261 260 OnSwap = detMin ==0; … … 265 264 if(kopt) { 266 265 // critere de Delaunay pure isotrope 267 Icoor2xb1 = sb->i.x - s1->i.x,266 long long xb1 = sb->i.x - s1->i.x, 268 267 x21 = s2->i.x - s1->i.x, 269 268 yb1 = sb->i.y - s1->i.y, -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/Triangle.h ¶
r21759 r22085 21 21 22 22 public: 23 Icoor2det; //Integer determinant (twice its area)23 long long det; //Integer determinant (twice its area) 24 24 union { 25 25 Triangle *link; … … 60 60 61 61 //Inline methods 62 double qualite() ;63 62 void Set(const Triangle &,const Mesh &,Mesh &); 64 63 int In(BamgVertex *v) const { return vertices[0]==v || vertices[1]==v || vertices[2]==v ;} -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/det.h ¶
r15061 r22085 6 6 namespace bamg { 7 7 8 Icoor2inline det(const I2 &a,const I2 & b,const I2 &c){9 Icoor2bax = b.x - a.x ,bay = b.y - a.y;10 Icoor2cax = c.x - a.x ,cay = c.y - a.y;8 long long inline det(const I2 &a,const I2 & b,const I2 &c){ 9 long long bax = b.x - a.x ,bay = b.y - a.y; 10 long long cax = c.x - a.x ,cay = c.y - a.y; 11 11 return bax*cay - bay*cax; 12 12 } -
TabularUnified issm/branches/trunk-larour-NatGeoScience2016/src/c/bamg/typedefs.h ¶
r15101 r22085 7 7 8 8 /*Integer coordinates types*/ 9 typedef int Icoor1;10 typedef long long Icoor2;11 9 12 10 /*I2 and R2*/ 13 typedef P2< Icoor1,Icoor2> I2;11 typedef P2<int,long long> I2; 14 12 typedef P2<double,double> R2; 15 13 }
Note:
See TracChangeset
for help on using the changeset viewer.
