![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21837
- Timestamp:
- 07/24/17 10:34:45 (8 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 11 edited
-
m4/issm_options.m4 (modified) (1 diff)
-
src/c/Makefile.am (modified) (2 diffs)
-
src/c/bamg/BamgOpts.cpp (modified) (3 diffs)
-
src/c/bamg/BamgOpts.h (modified) (2 diffs)
-
src/c/bamg/Geometry.cpp (modified) (5 diffs)
-
src/c/bamg/Geometry.h (modified) (1 diff)
-
src/c/bamg/Mesh.cpp (modified) (16 diffs)
-
src/c/bamg/Mesh.h (modified) (6 diffs)
-
src/c/bamg/Triangle.cpp (modified) (1 diff)
-
src/c/classes/AmrBamg.cpp (modified) (2 diffs)
-
src/c/modules/Bamgx/Bamgx.cpp (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/m4/issm_options.m4
r21767 r21837 2060 2060 fi 2061 2061 2062 dnl check that fortran is provided if GiaIvins is on 2063 if test "$HAVE_GIAIVINS" = "yes" && test "$HAVE_FORTRAN" = "no" ; then 2064 AC_MSG_ERROR([need fortran compiler to compile GiaIvins (or use --without-GiaIvins )]); 2065 fi 2066 2062 2067 dnl check that if we have MPI, we have metis 2063 2068 if test "$HAVE_METIS" = "yes" && test "$HAVE_MPI" = "no" ; then -
issm/trunk-jpl/src/c/Makefile.am
r21802 r21837 454 454 endif 455 455 #}}} 456 #Gia sources {{{456 #Gia sources (only if have fortran){{{ 457 457 if GIAIVINS 458 if FORTRAN 458 459 issm_sources += ./cores/gia_core.cpp\ 459 460 ./analyses/GiaIvinsAnalysis.cpp\ … … 466 467 ./modules/GiaDeflectionCorex/stot.f\ 467 468 ./modules/GiaDeflectionCorex/what0.f 469 endif 468 470 endif 469 471 #}}} -
issm/trunk-jpl/src/c/bamg/BamgOpts.cpp
r21802 r21837 11 11 this->gradation = 0; 12 12 this->Hessiantype = 0; 13 this->MaxCornerAngle = 0;14 13 this->maxnbv = 0; 15 14 this->maxsubdiv = 0; … … 19 18 this->omega = 0; 20 19 this->power = 0; 21 this->random = 0;22 20 this->verbose = 0; 23 21 24 22 this->Crack = 0; 25 this->geometricalmetric = 0;26 23 this->KeepVertices = 0; 27 24 this->splitcorners = 0; … … 68 65 if (this->Crack!=0 && this->Crack!=1) _error_("'Crack' supported options are 0 and 1"); 69 66 if (this->KeepVertices!=0 && this->KeepVertices!=1) _error_("'KeepVertices' supported options are 0 and 1"); 70 if (this->geometricalmetric!=0 && this->geometricalmetric!=1) _error_("'geometricalmetric' supported options are 0 and 1");71 67 72 68 if (this->hmin<=0) _error_("'hmin' option should be >0"); -
issm/trunk-jpl/src/c/bamg/BamgOpts.h
r16967 r21837 17 17 double gradation; 18 18 int Hessiantype; 19 double MaxCornerAngle;20 19 int maxnbv; 21 20 double maxsubdiv; … … 25 24 double omega; 26 25 double power; 27 bool random;28 26 int verbose; 29 27 30 28 /*Flags*/ 31 29 int Crack; 32 int geometricalmetric;33 30 int KeepVertices; 34 31 int splitcorners; -
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r21629 r21837 187 187 } 188 188 189 //MaxCornerAngle190 if (bamgopts->MaxCornerAngle){191 if(verbose>5) _printf_(" processing MaxCornerAngle\n");192 MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;193 }194 195 189 //TangentAtEdges 196 190 if (bamggeom->TangentAtEdges){ … … 401 395 _printf_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")\n"); 402 396 _printf_(" coefIcoor: " << coefIcoor << "\n"); 403 _printf_(" MaxCornerAngle: " << MaxCornerAngle << "\n");404 397 405 398 return; … … 419 412 nbcurves=0; 420 413 subdomains=NULL; 421 MaxCornerAngle = 10*Pi/180; //default is 10 degres422 414 } 423 415 /*}}}*/ … … 614 606 } 615 607 616 // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) 617 if(ord==2) { 608 //all vertices provided in geometry are corners (ord = number of edges holding i) 609 vertices[i].SetCorner() ; 610 if(ord==2){ 618 611 long n1 = head_v[i]; 619 612 long n2 = next_p[n1]; 620 613 long i1 = n1/2, i2 = n2/2; // edge number 621 long j1 = n1%2, j2 = n2%2; // vertex in the edge 622 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; 623 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; 624 float da12 = Abs(angle2-angle1); 625 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) { 626 vertices[i].SetCorner() ; 627 } 628 // if the edge type/referencenumber a changing then is SetRequired(); 614 long j1 = n1%2, j2 = n2%2; // vertex in the edge 615 /* if the edge type/referencenumber a changing then is SetRequired();*/ 629 616 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 630 617 vertices[i].SetRequired(); … … 633 620 vertices[i].SetRequired(); 634 621 } 635 }636 if(ord != 2) {637 vertices[i].SetCorner();638 622 } 639 623 -
issm/trunk-jpl/src/c/bamg/Geometry.h
r21629 r21837 32 32 R2 pmin,pmax; // domain extrema coordinates 33 33 double coefIcoor; // coef to integer coordinates; 34 double MaxCornerAngle;35 34 36 35 //Constructor/Destructors -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r21821 r21837 980 980 981 981 /*Methods*/ 982 void Mesh::AddGeometryMetric(BamgOpts* bamgopts){/*{{{*/983 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/IntersectGeomMetric)*/984 985 /*Get options*/986 double anisomax = bamgopts->anisomax;987 double errg = bamgopts->errg;988 989 double ss[2]={0.00001,0.99999};990 double errC = 2*sqrt(2*errg);991 double hmax = Gh.MaximalHmax();992 double hmin = Gh.MinimalHmin();993 994 //check that hmax is positive995 _assert_(hmax>0);996 997 //errC cannot be higher than 1998 if(errC>1) errC=1;999 1000 //Set all vertices to "on"1001 SetVertexFieldOn();1002 1003 //loop over all the vertices on edges1004 for(int i=0;i<nbe;i++){1005 for (int j=0;j<2;j++){1006 1007 BamgVertex V;1008 VertexOnGeom GV;1009 Gh.ProjectOnCurve(edges[i],ss[j],V,GV);1010 1011 GeomEdge* eg = GV;1012 double s = GV;1013 R2 tg;1014 double R1= eg->R1tg(s,tg);1015 double ht=hmax;1016 // err relative to the length of the edge1017 if (R1>1.0e-20) {1018 ht = Min(Max(errC/R1,hmin),hmax);1019 }1020 double hn=Min(hmax,ht*anisomax);1021 1022 if (ht<=0 || hn<=0){1023 _error_("ht<=0 || hn<=0");1024 }1025 EigenMetric Vp(1/(ht*ht),1/(hn*hn),tg);1026 Metric MVp(Vp);1027 edges[i][j].m.IntersectWith(MVp);1028 }1029 }1030 // the problem is for the vertex on vertex1031 }1032 /*}}}*/1033 982 void Mesh::AddMetric(BamgOpts* bamgopts){/*{{{*/ 1034 983 // Hessiantype = 0 => H is computed using double L2 projection … … 1231 1180 int i,j,k,kk,it,jt; 1232 1181 int verbose=0; 1233 double cutoffradian=10*Pi/180;1234 1182 1235 1183 /*Recover options*/ 1236 if 1184 if(bamgopts){ 1237 1185 verbose=bamgopts->verbose; 1238 cutoffradian=bamgopts->MaxCornerAngle*Pi/180;1239 1186 } 1240 1187 … … 1243 1190 1244 1191 //check that the mesh is not empty 1245 if (nbt<=0 || nbv <=0 ) { 1246 _error_("nbt or nbv is negative (Mesh empty?)"); 1247 } 1248 1249 //Gh is the geometry of the mesh (this), initialize MaxCornerAngle 1250 if (cutoffradian>=0) Gh.MaxCornerAngle = cutoffradian; 1192 if(nbt<=0 || nbv <=0 ) _error_("nbt or nbv is negative (Mesh empty?)"); 1251 1193 1252 1194 /*Construction of the edges*/ … … 1447 1389 1448 1390 //check that nbsubdomains is empty 1449 if (nbsubdomains){ 1450 _error_("nbsubdomains should be 0"); 1451 } 1391 if(nbsubdomains) _error_("nbsubdomains should be 0"); 1452 1392 nbsubdomains=0; 1453 1393 … … 2649 2589 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/PreInit)*/ 2650 2590 2651 /* initialize random seed: */2652 srand(19999999);2653 2654 2591 /*Initialize fields*/ 2655 NbRef=0; 2656 quadtree=NULL; 2657 nbv=0; 2658 nbt=0; 2659 nbe=0; 2660 edges=NULL; 2661 nbsubdomains=0; 2662 subdomains=NULL; 2663 maxnbv=maxnbv_in; 2664 maxnbt=2 *maxnbv_in-2; 2665 NbVertexOnBThVertex=0; 2666 VertexOnBThVertex=NULL; 2667 NbVertexOnBThEdge=0; 2668 VertexOnBThEdge=NULL; 2669 NbCrackedVertices=0; 2670 CrackedVertices =NULL; 2671 NbCrackedEdges =0; 2672 CrackedEdges =NULL; 2673 NbVerticesOnGeomVertex=0; 2674 VerticesOnGeomVertex=NULL; 2675 NbVerticesOnGeomEdge=0; 2676 VerticesOnGeomEdge=NULL; 2592 this->NbRef = 0; 2593 this->quadtree = NULL; 2594 this->nbv = 0; 2595 this->nbt = 0; 2596 this->nbe = 0; 2597 this->edges = NULL; 2598 this->nbsubdomains = 0; 2599 this->subdomains = NULL; 2600 this->maxnbv = maxnbv_in; 2601 this->maxnbt = 2 *maxnbv_in-2; 2602 this->NbVertexOnBThVertex = 0; 2603 this->VertexOnBThVertex = NULL; 2604 this->NbVertexOnBThEdge = 0; 2605 this->VertexOnBThEdge = NULL; 2606 this->NbCrackedVertices = 0; 2607 this->CrackedVertices = NULL; 2608 this->NbCrackedEdges = 0; 2609 this->CrackedEdges = NULL; 2610 this->NbVerticesOnGeomVertex = 0; 2611 this->VerticesOnGeomVertex = NULL; 2612 this->NbVerticesOnGeomEdge = 0; 2613 this->VerticesOnGeomEdge = NULL; 2614 2615 /*Initialize random seed*/ 2616 this->randomseed = 1; 2677 2617 2678 2618 /*Allocate if maxnbv_in>0*/ 2679 2619 if(maxnbv_in){ 2680 vertices=new BamgVertex[maxnbv];2681 _assert_(vertices);2682 orderedvertices=new BamgVertex* [maxnbv];2683 _assert_(orderedvertices);2684 triangles=new Triangle[maxnbt];2685 _assert_(t riangles);2620 this->vertices=new BamgVertex[this->maxnbv]; 2621 this->orderedvertices=new BamgVertex* [this->maxnbv]; 2622 this->triangles=new Triangle[this->maxnbt]; 2623 _assert_(this->vertices); 2624 _assert_(this->orderedvertices); 2625 _assert_(this->triangles); 2686 2626 } 2687 2627 else{ 2688 vertices=NULL;2689 orderedvertices=NULL;2690 t riangles=NULL;2691 maxnbt=0;2628 this->vertices = NULL; 2629 this->orderedvertices = NULL; 2630 this->triangles = NULL; 2631 this->maxnbt = 0; 2692 2632 } 2693 2633 } 2694 2634 /*}}}*/ 2695 void Mesh::Insert( bool random){/*{{{*/2635 void Mesh::Insert(){/*{{{*/ 2696 2636 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Insert)*/ 2697 2637 … … 2740 2680 //Get Prime number 2741 2681 const long PrimeNumber= BigPrimeNumber(nbv); 2742 int temp = rand(); if(!random) temp = 756804110; 2743 int k0=temp%nbv; 2682 int k0=this->RandomNumber(nbv); 2744 2683 2745 2684 //Build orderedvertices … … 2825 2764 } 2826 2765 /*}}}*/ 2827 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap ,bool random) {/*{{{*/2766 long Mesh::InsertNewPoints(long nbvold,long & NbTSwap) {/*{{{*/ 2828 2767 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/InsertNewPoints)*/ 2829 2768 … … 2845 2784 /*construction of a random order*/ 2846 2785 const long PrimeNumber= BigPrimeNumber(nbv) ; 2847 //remainder of the division of rand() by nbvnew 2848 int temp = rand(); if(!random) temp = 1098566905; 2849 long k3 = temp%nbvnew; 2786 long k3 = long(this->RandomNumber(nbvnew)); 2850 2787 //loop over the new points 2851 2788 for (int is3=0; is3<nbvnew; is3++){ … … 3088 3025 //if(pointsoutside) _printf_("WARNING: One or more points of the initial mesh fall outside of the geometric boundary\n"); 3089 3026 Bh.CreateSingleVertexToTriangleConnectivity(); 3090 InsertNewPoints(nbvold,NbTSwap ,bamgopts->random);3027 InsertNewPoints(nbvold,NbTSwap); 3091 3028 } 3092 3029 else Bh.CreateSingleVertexToTriangleConnectivity(); … … 3146 3083 }// for triangle 3147 3084 3148 if (!InsertNewPoints(nbvold,NbTSwap ,bamgopts->random)) break;3085 if (!InsertNewPoints(nbvold,NbTSwap)) break; 3149 3086 for (i=nbtold;i<nbt;i++) first_np_or_next_t[i]=iter; 3150 3087 Headt = nbt; // empty list … … 4095 4032 4096 4033 /*Insert Vertices*/ 4097 Insert( true);4034 Insert(); 4098 4035 } 4099 4036 /*}}}*/ … … 4396 4333 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4397 4334 if (verbose>3) _printf_(" Inserting boundary points\n"); 4398 Insert( bamgopts->random);4335 Insert(); 4399 4336 4400 4337 //Force the boundary … … 4722 4659 if (verbose>3) _printf_(" Creating initial Constrained Delaunay Triangulation...\n"); 4723 4660 if (verbose>3) _printf_(" Inserting boundary points\n"); 4724 Insert( bamgopts->random);4661 Insert(); 4725 4662 4726 4663 //Force the boundary … … 4735 4672 NewPoints(BTh,bamgopts,KeepVertices) ; 4736 4673 if (verbose>4) _printf_(" -- current number of vertices = " << nbv << "\n"); 4674 } 4675 /*}}}*/ 4676 int Mesh::RandomNumber(int max){/*{{{*/ 4677 /* Generate a random number from 0 to max-1 using linear congruential 4678 * random number generator*/ 4679 4680 this->randomseed = (this->randomseed * 1366l + 150889l) % 714025l; 4681 return this->randomseed / (714025l / max + 1); 4682 } /*}}}*/ 4683 int Mesh::ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/ 4684 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/ 4685 4686 int NbSwap =0; 4687 if (!a.t || !b.t){ // the 2 vertex is in a mesh 4688 _error_("!a.t || !b.t"); 4689 } 4690 int k=0; 4691 taret=AdjacentTriangle(0,0); // erreur 4692 4693 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]); 4694 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; 4695 // we turn around a in the direct direction 4696 4697 long long det2 = v2 ? det(*v2,a,b): -1 , det1; 4698 if(v2) // normal case 4699 det2 = det(*v2,a,b); 4700 else { // no chance infini vertex try the next 4701 tta= Previous(Adj(tta)); 4702 v2 = tta.EdgeVertex(0); 4703 vbegin =v2; 4704 if (!v2){ 4705 _error_("!v2"); 4706 } 4707 det2 = det(*v2,a,b); 4708 } 4709 4710 while (v2 != &b) { 4711 AdjacentTriangle tc = Previous(Adj(tta)); 4712 v1 = v2; 4713 v2 = tc.EdgeVertex(0); 4714 det1 = det2; 4715 det2 = v2 ? det(*v2,a,b): det2; 4716 4717 if((det1 < 0) && (det2 >0)) { 4718 // try to force the edge 4719 BamgVertex * va = &a, *vb = &b; 4720 tc = Previous(tc); 4721 if (!v1 || !v2){ 4722 _error_("!v1 || !v2"); 4723 } 4724 long long detss = 0,l=0; 4725 while ((this->SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) 4726 if(l++ > 10000000) { 4727 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big"); 4728 } 4729 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); 4730 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){ 4731 tc.SetLock(); 4732 a.Optim(1,0); 4733 b.Optim(1,0); 4734 taret = tc; 4735 return NbSwap; 4736 } 4737 else 4738 { 4739 taret = tc; 4740 return -2; // error boundary is crossing 4741 } 4742 } 4743 tta = tc; 4744 k++; 4745 if (k>=2000){ 4746 _error_("k>=2000"); 4747 } 4748 if ( vbegin == v2 ) return -1;// error 4749 } 4750 4751 tta.SetLock(); 4752 taret=tta; 4753 a.Optim(1,0); 4754 b.Optim(1,0); 4755 return NbSwap; 4756 } 4757 /*}}}*/ 4758 int Mesh::SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/ 4759 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/ 4760 // l'arete ta coupe l'arete pva pvb 4761 // de cas apres le swap sa coupe toujours 4762 // on cherche l'arete suivante 4763 // on suppose que detsa >0 et detsb <0 4764 // attention la routine echange pva et pvb 4765 4766 if(tt1.Locked()) return 0; // frontiere croise 4767 4768 AdjacentTriangle tt2 = Adj(tt1); 4769 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent 4770 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles 4771 if ( a1<0 || a1>=3 ){ 4772 _error_("a1<0 || a1>=3"); 4773 } 4774 4775 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; 4776 BamgVertex & s1= (*t1)[OppositeVertex[a1]]; 4777 BamgVertex & s2= (*t2)[OppositeVertex[a2]]; 4778 4779 long long dets2 = det(*pva,*pvb,s2); 4780 long long det1=t1->det , det2=t2->det ; 4781 long long detT = det1+det2; 4782 if ((det1<=0 ) || (det2<=0)){ 4783 _error_("(det1<=0 ) || (det2<=0)"); 4784 } 4785 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb 4786 _error_("(detsa>=0) || (detsb<=0)"); 4787 } 4788 long long ndet1 = bamg::det(s1,sa,s2); 4789 long long ndet2 = detT - ndet1; 4790 4791 int ToSwap =0; //pas de swap 4792 if ((ndet1 >0) && (ndet2 >0)) 4793 { // on peut swaper 4794 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) 4795 ToSwap =1; 4796 else // swap alleatoire 4797 if (BinaryRand()) 4798 ToSwap =2; 4799 } 4800 if (ToSwap) NbSwap++, 4801 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); 4802 4803 int ret=1; 4804 4805 if (dets2 < 0) {// haut 4806 dets1 = ToSwap ? dets1 : detsa ; 4807 detsa = dets2; 4808 tt1 = Previous(tt2) ;} 4809 else if (dets2 > 0){// bas 4810 dets1 = ToSwap ? dets1 : detsb ; 4811 detsb = dets2; 4812 //xxxx tt1 = ToSwap ? tt1 : Next(tt2); 4813 if(!ToSwap) tt1 = Next(tt2); 4814 } 4815 else { // changement de direction 4816 ret = -1; 4817 Exchange(pva,pvb); 4818 Exchange(detsa,detsb); 4819 Exchange(dets1,dets2); 4820 Exchange(tt1,tt2); 4821 dets1=-dets1; 4822 dets2=-dets2; 4823 detsa=-detsa; 4824 detsb=-detsb; 4825 4826 if(ToSwap){ 4827 if (dets2 < 0) {// haut 4828 dets1 = (ToSwap ? dets1 : detsa) ; 4829 detsa = dets2; 4830 tt1 = Previous(tt2) ;} 4831 else if(dets2 > 0){// bas 4832 dets1 = (ToSwap ? dets1 : detsb) ; 4833 detsb = dets2; 4834 if(!ToSwap) tt1 = Next(tt2); 4835 } 4836 else {// on a fin ??? 4837 tt1 = Next(tt2); 4838 ret =0;} 4839 } 4840 4841 } 4842 return ret; 4737 4843 } 4738 4844 /*}}}*/ … … 4776 4882 return edge; 4777 4883 } 4778 }4779 /*}}}*/4780 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) { /*{{{*/4781 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/ForceEdge)*/4782 4783 int NbSwap =0;4784 if (!a.t || !b.t){ // the 2 vertex is in a mesh4785 _error_("!a.t || !b.t");4786 }4787 int k=0;4788 taret=AdjacentTriangle(0,0); // erreur4789 4790 AdjacentTriangle tta(a.t,EdgesVertexTriangle[a.IndexInTriangle][0]);4791 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;4792 // we turn around a in the direct direction4793 4794 long long det2 = v2 ? det(*v2,a,b): -1 , det1;4795 if(v2) // normal case4796 det2 = det(*v2,a,b);4797 else { // no chance infini vertex try the next4798 tta= Previous(Adj(tta));4799 v2 = tta.EdgeVertex(0);4800 vbegin =v2;4801 if (!v2){4802 _error_("!v2");4803 }4804 det2 = det(*v2,a,b);4805 }4806 4807 while (v2 != &b) {4808 AdjacentTriangle tc = Previous(Adj(tta));4809 v1 = v2;4810 v2 = tc.EdgeVertex(0);4811 det1 = det2;4812 det2 = v2 ? det(*v2,a,b): det2;4813 4814 if((det1 < 0) && (det2 >0)) {4815 // try to force the edge4816 BamgVertex * va = &a, *vb = &b;4817 tc = Previous(tc);4818 if (!v1 || !v2){4819 _error_("!v1 || !v2");4820 }4821 long long detss = 0,l=0;4822 while ((SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap)))4823 if(l++ > 10000000) {4824 _error_("Loop in forcing Egde, nb de swap=" << NbSwap << ", nb of try swap (" << l << ") too big");4825 }4826 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);4827 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){4828 tc.SetLock();4829 a.Optim(1,0);4830 b.Optim(1,0);4831 taret = tc;4832 return NbSwap;4833 }4834 else4835 {4836 taret = tc;4837 return -2; // error boundary is crossing4838 }4839 }4840 tta = tc;4841 k++;4842 if (k>=2000){4843 _error_("k>=2000");4844 }4845 if ( vbegin == v2 ) return -1;// error4846 }4847 4848 tta.SetLock();4849 taret=tta;4850 a.Optim(1,0);4851 b.Optim(1,0);4852 return NbSwap;4853 4884 } 4854 4885 /*}}}*/ … … 4897 4928 } // end swap 4898 4929 /*}}}*/ 4899 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb ,AdjacentTriangle & tt1,long long & dets1, long long & detsa,long long & detsb, int & NbSwap) {/*{{{*/4900 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SwapForForcingEdge)*/4901 // l'arete ta coupe l'arete pva pvb4902 // de cas apres le swap sa coupe toujours4903 // on cherche l'arete suivante4904 // on suppose que detsa >0 et detsb <04905 // attention la routine echange pva et pvb4906 4907 if(tt1.Locked()) return 0; // frontiere croise4908 4909 AdjacentTriangle tt2 = Adj(tt1);4910 Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent4911 short a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles4912 if ( a1<0 || a1>=3 ){4913 _error_("a1<0 || a1>=3");4914 }4915 4916 BamgVertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]];4917 BamgVertex & s1= (*t1)[OppositeVertex[a1]];4918 BamgVertex & s2= (*t2)[OppositeVertex[a2]];4919 4920 long long dets2 = det(*pva,*pvb,s2);4921 long long det1=t1->det , det2=t2->det ;4922 long long detT = det1+det2;4923 if ((det1<=0 ) || (det2<=0)){4924 _error_("(det1<=0 ) || (det2<=0)");4925 }4926 if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb4927 _error_("(detsa>=0) || (detsb<=0)");4928 }4929 long long ndet1 = bamg::det(s1,sa,s2);4930 long long ndet2 = detT - ndet1;4931 4932 int ToSwap =0; //pas de swap4933 if ((ndet1 >0) && (ndet2 >0))4934 { // on peut swaper4935 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))4936 ToSwap =1;4937 else // swap alleatoire4938 if (BinaryRand())4939 ToSwap =2;4940 }4941 if (ToSwap) NbSwap++,4942 bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2);4943 4944 int ret=1;4945 4946 if (dets2 < 0) {// haut4947 dets1 = ToSwap ? dets1 : detsa ;4948 detsa = dets2;4949 tt1 = Previous(tt2) ;}4950 else if (dets2 > 0){// bas4951 dets1 = ToSwap ? dets1 : detsb ;4952 detsb = dets2;4953 //xxxx tt1 = ToSwap ? tt1 : Next(tt2);4954 if(!ToSwap) tt1 = Next(tt2);4955 }4956 else { // changement de direction4957 ret = -1;4958 Exchange(pva,pvb);4959 Exchange(detsa,detsb);4960 Exchange(dets1,dets2);4961 Exchange(tt1,tt2);4962 dets1=-dets1;4963 dets2=-dets2;4964 detsa=-detsa;4965 detsb=-detsb;4966 4967 if(ToSwap){4968 if (dets2 < 0) {// haut4969 dets1 = (ToSwap ? dets1 : detsa) ;4970 detsa = dets2;4971 tt1 = Previous(tt2) ;}4972 else if(dets2 > 0){// bas4973 dets1 = (ToSwap ? dets1 : detsb) ;4974 detsb = dets2;4975 if(!ToSwap) tt1 = Next(tt2);4976 }4977 else {// on a fin ???4978 tt1 = Next(tt2);4979 ret =0;}4980 }4981 4982 }4983 return ret;4984 }4985 /*}}}*/4986 4930 } -
issm/trunk-jpl/src/c/bamg/Mesh.h
r21791 r21837 40 40 double coefIcoor; // coef to integer 41 41 ListofIntersectionTriangles lIntTria; 42 int randomseed; //used for random number generation 42 43 43 44 long NbVerticesOnGeomVertex; … … 77 78 R2 I2ToR2(const I2 & P) const; 78 79 void AddVertex(BamgVertex & s,Triangle * t,long long * =0) ; 79 void Insert( bool random);80 void Insert(); 80 81 void Echo(void); 81 82 void ForceBoundary(); … … 91 92 void MaxSubDivision(double maxsubdiv); 92 93 void NewPoints(Mesh &,BamgOpts* bamgopts,int KeepVertices=1); 93 long InsertNewPoints(long nbvold,long & NbTSwap ,bool random);94 long InsertNewPoints(long nbvold,long & NbTSwap); 94 95 void TrianglesRenumberBySubDomain(bool justcompress=false); 95 96 void SmoothingVertex(int =3,double=0.3); … … 113 114 void BuildMetric0(BamgOpts* bamgopts); 114 115 void BuildMetric1(BamgOpts* bamgopts); 115 void AddGeometryMetric(BamgOpts* bamgopts);116 116 void BuildGeometryFromMesh(BamgOpts* bamgopts=NULL); 117 int RandomNumber(int max); 117 118 void ReconstructExistingMesh(); 118 119 … … 143 144 void Triangulate(double* x,double* y,int nods); 144 145 void Init(long); 146 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; 147 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , 148 AdjacentTriangle & tt1,long long & dets1, 149 long long & detsa,long long & detsb, int & nbswap); 145 150 }; 146 151 … … 150 155 Triangle *t2,short a2, 151 156 BamgVertex *s1,BamgVertex *s2,long long det1,long long det2); 152 int SwapForForcingEdge(BamgVertex * & pva ,BamgVertex * & pvb , 153 AdjacentTriangle & tt1,long long & dets1, 154 long long & detsa,long long & detsb, int & nbswap); 155 int ForceEdge(BamgVertex &a, BamgVertex & b,AdjacentTriangle & taret) ; 157 156 158 inline AdjacentTriangle Previous(const AdjacentTriangle & ta){ 157 159 return AdjacentTriangle(ta.t,PreviousEdge[ta.a]); -
issm/trunk-jpl/src/c/bamg/Triangle.cpp
r21791 r21837 257 257 258 258 long long detMinNew=Min(det1,det2); 259 // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test260 259 if (! OnSwap &&(detMinNew>0)) { 261 260 OnSwap = detMin ==0; -
issm/trunk-jpl/src/c/classes/AmrBamg.cpp
r21822 r21837 32 32 this->options->gradation = 1.5; 33 33 this->options->Hessiantype = 0; 34 this->options->MaxCornerAngle = 1.e-12;35 34 this->options->maxnbv = 1e6; 36 35 this->options->maxsubdiv = 10; … … 40 39 this->options->omega = 1.8; 41 40 this->options->power = 1; 42 this->options->random = 0;43 41 this->options->verbose = 0; 44 42 this->options->Crack = 0; 45 this->options->geometricalmetric = 0;46 43 this->options->KeepVertices = 1; /*!!!!! VERY IMPORTANT !!!!!*/ 47 44 this->options->splitcorners = 1; -
issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
r21623 r21837 149 149 } 150 150 151 //Add geometry metric if provided152 if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts);153 154 151 //Smoothe metric 155 152 BTh.SmoothMetric(bamgopts->gradation);
Note:
See TracChangeset
for help on using the changeset viewer.
