Mesh extraction with Polygon
During the data acquitition, it is possible that some areas aren’t completed or present some errors in it. These errors can trouble the solve of our problem, that’s why we want to exclude them.
The first method we went into is the extraction by polygons. As the name suppose, we define polygons that cover the problematic areas and we delete them in our mesh.
1. Method & Implementation
We define the polygons in a text file like that
1 // number of polygons
4 // number of nodes
132.812 235.938 // coordinates
349.609 93.3594
496.094 261.328
328.125 360.938We will then read and recover these data during the construction of our mesh, then build with the help of boost::geometry a polygon object, if we activate the option withPoly.
typedef boost::geometry::model::d2::point_xy<double> poly_point_type;
typedef boost::geometry::model::segment<poly_point_type> segment_type;
typedef boost::geometry::model::polygon<poly_point_type> polygon_type;
....
// polygon build
std::vector<polygon_type> list_poly ;
//std::cout << list_poly.size() << std::endl;
polygon_type poly;
bool inPoly = false;
bool testInPoly = false;
if (withPoly)
{
polygon_type p ;
std::ifstream flux(pathPoly) ;
int nbPolygons ;
int nbVertices ;
double coordX, coordY ;
flux >> nbPolygons ;
//std::cout<< "nP " << nbPolygons << std::endl;
for (int j = 0 ; j < nbPolygons ; j ++) {
flux >> nbVertices ;
//std::cout<< "nV: " << nbVertices << std::endl;
for (int i = 0 ; i < nbVertices ; i ++) {
flux >> coordX ;
flux >> coordY ;
p.outer().push_back(poly_point_type(coordX*pixelsize, coordY*pixelsize)) ;
}
p.outer().push_back(poly_point_type(p.outer().at(0).x(), p.outer().at(0).y())) ;
list_poly.push_back(p) ;
p.clear() ;
}
}Then for each mesh node, we test if it is within our polygon object. If it is, we didn’t add it to the mesh. If not, we build it normally.
if (withPoly)
{
poly_point_type p(coords[0],coords[1]);
for (int k=0;k<list_poly.size();k++)
{
poly= list_poly[k];
testInPoly=boost::geometry::within(p,poly);
inPoly=(inPoly)||(testInPoly);
}
}
if ( (!withPoly)||(!inPoly))
{
point_type pt( ptid,coords );
pt.setProcessId( partId );
pt.setProcessIdInPartition( partId );
this->addPoint( pt );
}
inPoly=false;We apply the same idea on elements : we test each node of the current element. If at least one of them is include in our polygons, we don’t create the element and we pass to the next one. If they are all out, we add the element normally.
for (int k=0;k<list_poly.size();k++)
{
poly= list_poly[k];
bool testInPoly1=boost::geometry::within(p1,poly);
bool testInPoly2=boost::geometry::within(p2,poly);
bool testInPoly3=boost::geometry::within(p3,poly);
bool testInPoly4=boost::geometry::within(p4,poly);
inPoly=(inPoly)||(testInPoly1)||(testInPoly2)||(testInPoly3)||(testInPoly4);
}
}
if ( (!withPoly)||(!inPoly))
{
size_type eid = (M_ny-1)*i+j; // StructuredMesh Id
element_type e;
e.setProcessIdInPartition( partId );
e.setProcessId( partId );
size_type ptid[4] = { (M_ny)*(i+1)+j, // 0
(M_ny)*(i+1)+j+1, // 1
(M_ny)*i+j+1, // 2
(M_ny)*i+j }; // 3
for( uint16_type k = 0; k < 4; ++k )
e.setPoint( k, this->point( ptid[k] ) );
if ( neighborProcessId != invalid_rank_type_value )
e.addNeighborPartitionId( neighborProcessId );
auto const& eltInserted = this->addElement( e, true ); // e.id() is defined by Feel++
idStructuredMeshToFeelMesh.insert( std::make_pair(eid,eltInserted.id()));
}
inPoly=false;