![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22380
- Timestamp:
- 01/29/18 14:39:39 (7 years ago)
- Location:
- issm/trunk-jpl/src/c/bamg
- Files:
-
- 3 edited
-
Geometry.cpp (modified) (2 diffs)
-
Geometry.h (modified) (1 diff)
-
Mesh.cpp (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/bamg/Geometry.cpp
r22269 r22380 438 438 return c - curves; 439 439 }/*}}}*/ 440 void Geometry::PostRead( ){/*{{{*/440 void Geometry::PostRead(bool checkcurve){/*{{{*/ 441 441 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/ 442 442 … … 606 606 } 607 607 608 //all vertices provided in geometry are corners (ord = number of edges holding i) 609 vertices[i].SetCorner() ; 610 if(ord==2){ 611 long n1 = head_v[i]; 612 long n2 = next_p[n1]; 613 long i1 = n1/2, i2 = n2/2; // edge number 614 long j1 = n1%2, j2 = n2%2; // vertex in the edge 615 /* if the edge type/referencenumber a changing then is SetRequired();*/ 616 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 617 vertices[i].SetRequired(); 618 } 619 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) { 620 vertices[i].SetRequired(); 608 /*Do we want to check for curve? Default is no, but if we are reconstructing, it's better to turn it on with a small threshold*/ 609 if(checkcurve){ 610 /*angular test on current vertex to guess whether it is a corner (ord = number of edges holding i) */ 611 if(ord==2){ 612 IssmDouble MaxCornerAngle = 1*Pi/180; /*default is 1 degree*/ 613 long n1 = head_v[i]; 614 long n2 = next_p[n1]; 615 long i1 = n1/2, i2 = n2/2; // edge number 616 long j1 = n1%2, j2 = n2%2; // vertex in the edge 617 float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1]; 618 float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2]; 619 float da12 = Abs(angle2-angle1); 620 if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) { 621 vertices[i].SetCorner() ; 622 } 623 /* if the edge type/referencenumber a changing then is SetRequired();*/ 624 if(edges[i1].type != edges[i2].type || edges[i1].Required()){ 625 vertices[i].SetRequired(); 626 } 627 if(edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) { 628 vertices[i].SetRequired(); 629 } 630 } 631 if(ord!=2) vertices[i].SetCorner(); 632 } 633 else{ 634 /*all vertices provided in geometry are corners (ord = number of edges holding i)*/ 635 vertices[i].SetCorner() ; 636 if(ord==2){ 637 long n1 = head_v[i]; 638 long n2 = next_p[n1]; 639 long i1 = n1/2, i2 = n2/2; // edge number 640 long j1 = n1%2, j2 = n2%2; // vertex in the edge 641 /* if the edge type/referencenumber a changing then is SetRequired();*/ 642 if (edges[i1].type != edges[i2].type || edges[i1].Required()){ 643 vertices[i].SetRequired(); 644 } 645 if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) { 646 vertices[i].SetRequired(); 647 } 621 648 } 622 649 } -
issm/trunk-jpl/src/c/bamg/Geometry.h
r21864 r22380 52 52 void ReadGeometry(BamgGeom *bamggeom, BamgOpts*bamgopts); 53 53 void Init(void); 54 void PostRead( );54 void PostRead(bool checkcurve=false); 55 55 long GetId(const GeomVertex &t) const; 56 56 long GetId(const GeomVertex *t) const; -
issm/trunk-jpl/src/c/bamg/Mesh.cpp
r22253 r22380 30 30 _printf_("WARNING: mesh present but no geometry found. Reconstructing...\n"); 31 31 BuildGeometryFromMesh(bamgopts); 32 Gh.PostRead( );32 Gh.PostRead(true); 33 33 } 34 34 … … 1215 1215 /*Construction of the edges*/ 1216 1216 1217 / /initialize st and edge41217 /*initialize st and edge4*/ 1218 1218 SetOfEdges4* edge4= new SetOfEdges4(nbt*3,nbv); 1219 1219 long* st = new long[nbt*3]; 1220 1220 1221 / /initialize st as -1 (chaining algorithm)1221 /*initialize st as -1 (chaining algorithm)*/ 1222 1222 for (i=0;i<nbt*3;i++) st[i]=-1; 1223 1223 1224 / /build edge4 (chain)1224 /*build edge4 (chain)*/ 1225 1225 for (i=0;i<nbe;i++){ 1226 1226 edge4->SortAndAdd(GetId(edges[i][0]),GetId(edges[i][1])); 1227 1227 } 1228 / /check that there is no double edge1228 /*check that there is no double edge*/ 1229 1229 if (nbe != edge4->nb()){ 1230 1230 delete [] st; … … 1234 1234 long nbeold = nbe; 1235 1235 1236 //Go through the triangles and a ssthe edges in edge4 if they are not there yet1236 //Go through the triangles and add the edges in edge4 if they are not there yet 1237 1237 for (i=0;i<nbt;i++){ 1238 1238 //3 edges per triangle … … 1251 1251 _error_("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)"); 1252 1252 } 1253 / /OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list1253 /*OK, the element is not on boundary, is belongs to 2 triangles -> build Adjacent triangles list*/ 1254 1254 triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); 1255 1255 if (invisible) triangles[i].SetHidden(j); 1256 / / if k < nbe mark the edge as on Boundary (Locked)1256 /* if k < nbe mark the edge as on Boundary (Locked)*/ 1257 1257 if (k<nbe) { 1258 1258 triangles[i].SetLocked(j); 1259 1259 } 1260 / /set st[k] as negative so that the edge should not be called again1260 /*set st[k] as negative so that the edge should not be called again*/ 1261 1261 st[k]=-2-st[k]; 1262 1262 } 1263 //else (see 3 lines above), the edge has been called more than twice: return error1264 1263 else { 1264 /*else (see 3 lines above), the edge has been called more than twice: return error*/ 1265 1265 _printf_("The edge (" << GetId(triangles[i][VerticesOfTriangularEdge[j][0]]) << "," << GetId(triangles[i][VerticesOfTriangularEdge[j][1]]) << ") belongs to more than 2 triangles (" << k << ")\n"); 1266 1266 _printf_("Edge " << j << " of triangle " << i << "\n"); … … 1286 1286 } 1287 1287 1288 / / check consistency of edge[].adj and geometrical required vertices1288 /*check consistency of edge[].adj and geometrical required vertices*/ 1289 1289 k=0; kk=0; 1290 1290 for (i=0;i<nbedges;i++){ … … 1304 1304 1305 1305 /*Constructions of edges*/ 1306 1307 1306 k += kk; 1308 1307 kk=0;
Note:
See TracChangeset
for help on using the changeset viewer.
