![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 24113
- Timestamp:
- 07/30/19 07:18:02 (6 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 4 edited
-
bamg/BamgOpts.cpp (modified) (2 diffs)
-
bamg/BamgOpts.h (modified) (2 diffs)
-
bamg/Geometry.cpp (modified) (29 diffs)
-
modules/Bamgx/Bamgx.cpp (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/bamg/BamgOpts.cpp
r21864 r24113 28 28 this->hminVertices=NULL; this->hminVerticesSize[0]=this->hminVerticesSize[1]=0; 29 29 this->hmaxVertices=NULL; this->hmaxVerticesSize[0]=this->hmaxVerticesSize[1]=0; 30 this->hVertices=NULL; this->hVertices Size[0]=this->hVerticesSize[1]=0;30 this->hVertices=NULL; this->hVerticesLength=0; 31 31 this->metric=NULL; this->metricSize[0]=this->metricSize[1]=0; 32 32 this->field=NULL; this->fieldSize[0]=this->fieldSize[1]=0; … … 70 70 if (this->hminVertices && this->hminVerticesSize[1]!=1) _error_("'hminVertices' should be a column"); 71 71 if (this->hmaxVertices && this->hmaxVerticesSize[1]!=1) _error_("'hmaxVertices' should be a column"); 72 if (this->hVertices && this->hVertices Size[1]!=1) _error_("'hVertices' should be a column");72 if (this->hVertices && this->hVerticesLength<2) _error_("'hVertices' should be a vector"); 73 73 if (this->metric && (this->metricSize[1]!=1 && this->metricSize[1]!=3)) _error_("'metric' should have either 1 (iso) or 3 (aniso) columns."); 74 74 if (this->field){ -
issm/trunk-jpl/src/c/bamg/BamgOpts.h
r21864 r24113 1 1 /*!\file: BamgOpts.h 2 2 * \brief place holder for optimization function arguments 3 */ 3 */ 4 4 5 5 #ifndef _BAMGOPTS_H_ … … 38 38 int hmaxVerticesSize[2]; 39 39 double* hmaxVertices; 40 int hVertices Size[2];40 int hVerticesLength; 41 41 double* hVertices; 42 42 int metricSize[2]; -
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r22897 r24113 57 57 nbcurves=0; 58 58 59 double Hmin = HUGE_VAL;// the infinie value 59 double Hmin = HUGE_VAL;// the infinie value 60 60 int i,j,n,i0,i1,i2,i3; 61 61 … … 96 96 /*coefIcoor is the coefficient used for integer coordinates: 97 97 * (x-pmin.x) 98 * Icoor x = (2^30 -1) ------------ 98 * Icoor x = (2^30 -1) ------------ 99 99 * D 100 100 * where D is the longest side of the domain (direction x or y) … … 155 155 } 156 156 157 // definition the default of the given mesh size 157 // definition the default of the given mesh size 158 158 for (i=0;i<nbv;i++) { 159 if (vertices[i].color > 0) 159 if (vertices[i].color > 0) 160 160 vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color); 161 else 161 else 162 162 vertices[i].m=Metric(Hmin); 163 163 } … … 170 170 171 171 //hVertices 172 if(bamgopts->hVertices && bamgopts->hVertices Size[0]==nbv){172 if(bamgopts->hVertices && bamgopts->hVerticesLength==nbv){ 173 173 if(verbose>5) _printf_(" processing hVertices\n"); 174 174 for (i=0;i< nbv;i++){ … … 211 211 if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column"); 212 212 n=bamggeom->CornersSize[0]; 213 for (i=0;i<n;i++) { 213 for (i=0;i<n;i++) { 214 214 j=(int)bamggeom->Corners[i]-1; //for C indexing 215 215 if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]"); … … 225 225 if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column"); 226 226 n=bamggeom->RequiredVerticesSize[0]; 227 for (i=0;i<n;i++) { 227 for (i=0;i<n;i++) { 228 228 j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing 229 229 if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]"); … … 237 237 if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column"); 238 238 n=bamggeom->RequiredEdgesSize[0]; 239 for (i=0;i<n;i++) { 239 for (i=0;i<n;i++) { 240 240 j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing 241 241 if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]"); 242 edges[j].SetRequired(); 242 edges[j].SetRequired(); 243 243 } 244 244 } … … 457 457 the vertices by groups of 5: 458 458 All the coordinates are transformed to ]0,1[^2 459 then, the integer coordinates are computed using 459 then, the integer coordinates are computed using 460 460 the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/ 461 vertices[i].i=R2ToI2(vertices[i].r); 461 vertices[i].i=R2ToI2(vertices[i].r); 462 462 463 463 /*find nearest vertex already present in the quadtree (NULL if empty)*/ 464 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 464 BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y); 465 465 466 466 /*if there is a vertex found that is to close to vertices[i] -> error*/ … … 506 506 * head_F[j]=i; 507 507 * } 508 * 508 * 509 509 * Then, we can go through all the elements that have for image j: 510 510 * for(i=head_F[j]; i!=-1; i=next_F[i]) … … 512 512 * stop the loop when i=-1 (end of the chain) 513 513 * iterate using i=next_F[i] (next element that have for image j) 514 * 514 * 515 515 * 2. How to use this algorithm here? 516 * 516 * 517 517 * Here F is a function that associates two vertices v0 and v1 for a given edge E 518 518 * We want to build the reciprocal function: what are the edges that contains 519 519 * a vertex v? 520 520 * To do so, we use the same chaining algorithm but there is a difficulty 521 * coming from the fact that for F we have a couple of vertices and not one 521 * coming from the fact that for F we have a couple of vertices and not one 522 522 * vertex. 523 * To overcome this difficulty, we use a global indexing exactly like in 523 * To overcome this difficulty, we use a global indexing exactly like in 524 524 * C/C++ so that 525 525 * a member of a 2-column-table can be described by one index p=i*2+j … … 566 566 //sort head_v by order of increasing edges angle 567 567 for (i=0;i<nbv;i++) { 568 int exch=1,ord=0; 568 int exch=1,ord=0; 569 569 570 570 //exchange vertices position in head_v and next_p till tey are sorted 571 571 while (exch){ 572 long *p=head_v+i; 573 long *po=p; 574 long n=*p; 572 long *p=head_v+i; 573 long *po=p; 574 long n=*p; 575 575 float angleold=-1000 ; // angle = - infinity 576 576 ord=0; exch=0; … … 584 584 585 585 //Next vertex index 586 n=*pn; 586 n=*pn; 587 587 588 588 //compute angle between horizontal axis and v0->v1 589 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; 589 float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1]; 590 590 591 591 //exchange if the current edge angle is smaller than the previous one … … 593 593 exch=1; 594 594 *pn=*po; // next_p[n] = n + 1 595 *po=*p; // 595 *po=*p; // 596 596 *p=n; // next_p[n+1] = n 597 597 po=pn; // po now point toward pn (invert next and current) … … 634 634 else{ 635 635 /*all vertices provided in geometry are corners (ord = number of edges holding i)*/ 636 vertices[i].SetCorner() ; 637 if(ord==2){ 636 vertices[i].SetCorner() ; 637 if(ord==2){ 638 638 long n1 = head_v[i]; 639 639 long n2 = next_p[n1]; … … 652 652 /*close the list around the vertex to have a circular loop*/ 653 653 long no=-1, ne = head_v[i]; 654 while (ne >=0) ne = next_p[no=ne]; 654 while (ne >=0) ne = next_p[no=ne]; 655 655 if(no>=0) next_p[no] = head_v[i]; 656 656 } … … 662 662 for (j=0;j<2;j++){ 663 663 664 long n1 = next_p[k++]; 664 long n1 = next_p[k++]; 665 665 long i1 = n1/2 ,j1=n1%2; 666 666 … … 695 695 ltg= lAB; 696 696 } 697 //else: a Corner no tangent => nothing to do 697 //else: a Corner no tangent => nothing to do 698 698 } 699 699 else{ … … 710 710 if (ltg2[0]!=0) edges[i].SetTgA(); 711 711 if (ltg2[1]!=0) edges[i].SetTgB(); 712 } 712 } 713 713 714 714 /* generation of all curves (from corner to corner)*/ … … 726 726 for (i=0;i<nbe;i++){ 727 727 728 GeomEdge & ei=edges[i]; 728 GeomEdge & ei=edges[i]; 729 729 for(j=0;j<2;j++){ 730 730 /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)): 731 731 * we do have the first edge of a new curve*/ 732 if (!ei.Mark() && (level || ei[j].Required())) { 732 if (!ei.Mark() && (level || ei[j].Required())) { 733 733 int k0=j,k1; 734 734 GeomEdge *e=&ei; 735 GeomVertex *a=(*e)(k0); // begin 735 GeomVertex *a=(*e)(k0); // begin 736 736 if(curves){ 737 737 curves[nbcurves].FirstEdge=e; … … 739 739 } 740 740 int nee=0; 741 for(;;){ 741 for(;;){ 742 742 nee++; 743 k1 = 1-k0; // next vertex of the edge 743 k1 = 1-k0; // next vertex of the edge 744 744 e->SetMark(); 745 745 nb_marked_edges++; … … 765 765 } 766 766 } 767 } 767 } 768 768 } 769 769 _assert_(nb_marked_edges && nbe); 770 770 //allocate if first step 771 771 if(step==0) curves=new Curve[nbcurves]; 772 } 772 } 773 773 774 774 /*clean up*/ … … 812 812 813 813 //Get edge direction and swap v0 and v1 if necessary 814 R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; 814 R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; 815 815 int OppositeSens = (V01,AB)<0; 816 816 int direction0=0,direction1=1; … … 828 828 directionge[bge]=1; 829 829 830 while(eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ 830 while(eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){ 831 831 if (bge<=0) { 832 832 if(NbTry) { … … 846 846 direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]); 847 847 } 848 while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) { 849 if(tge>=mxe ) { 848 while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) { 849 if(tge>=mxe ) { 850 850 _printf_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); 851 851 NbTry++; … … 871 871 872 872 double sg; 873 if (eg0 == eg1) { 873 if (eg0 == eg1) { 874 874 double s0=vg0,s1=vg1; 875 875 sg = s0*(1.0-s) + s*s1; … … 886 886 lge[i]=ll += Norme2(AA-BB); 887 887 AA=BB ;} 888 lge[tge]=ll+=Norme2(AA-V1); 888 lge[tge]=ll+=Norme2(AA-V1); 889 889 // search the geometrical edge 890 890 _assert_(s<=1.0); … … 900 900 } 901 901 on=ge[i]; 902 if (i==tge) 902 if (i==tge) 903 903 s1=vg1; 904 904 905 905 s =(ls-l0)/(l1-l0); 906 sg =s0*(1.0-s)+s*s1; 907 } 906 sg =s0*(1.0-s)+s*s1; 907 } 908 908 _assert_(on); 909 909 V.r= on->F(sg); … … 915 915 /*coefIcoor is the coefficient used for integer coordinates: 916 916 * (x-pmin.x) 917 * Icoor x = (2^30 -1) ------------ 917 * Icoor x = (2^30 -1) ------------ 918 918 * D 919 919 * where D is the longest side of the domain (direction x or y) … … 927 927 for (int i=0;i<nbe;i++) edges[i].SetUnMark(); 928 928 }/*}}}*/ 929 } 929 } -
issm/trunk-jpl/src/c/modules/Bamgx/Bamgx.cpp
r22252 r24113 21 21 int i; 22 22 int noerr=1; 23 double hminaniso=1e-100; 23 double hminaniso=1e-100; 24 24 Mesh* Thr=NULL; 25 25 Mesh* Thb=NULL; … … 31 31 verbosity=bamgopts->verbose; 32 32 33 // no metric -> no smoothing 34 if (bamgopts->metric==NULL) nbsmooth=0; 33 // no metric -> no smoothing 34 if (bamgopts->metric==NULL) nbsmooth=0; 35 35 36 36 /*If no mesh in input, generate one*/ … … 84 84 /*Anisotropic mesh adaptation {{{*/ 85 85 86 // read background mesh 86 // read background mesh 87 87 if (verbosity>0) _printf_("Anisotropic mesh adaptation\n"); 88 88 if (verbosity>1) _printf_(" Processing initial mesh and geometry...\n"); 89 Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts); 89 Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts); 90 90 91 91 //Make Quadtree from background mesh … … 101 101 BTh.ReadMetric(bamgopts); 102 102 } 103 else { 103 else { 104 104 if (verbosity>1) _printf_(" Generating initial Metric...\n"); 105 105 Metric Mhmax(bamgopts->hmax); … … 114 114 115 115 // change using hVertices if provided 116 if(bamgopts->hVertices && bamgopts->hVertices Size[0]==BTh.nbv){116 if(bamgopts->hVertices && bamgopts->hVerticesLength==BTh.nbv){ 117 117 if (verbosity>1) _printf_(" Merging Metric with hVertices...\n"); 118 118 for (i=0;i<BTh.nbv;i++){
Note:
See TracChangeset
for help on using the changeset viewer.
