36 #include <boost/geometry.hpp>
37 #include <boost/geometry/algorithms/assign.hpp>
38 #include <boost/geometry/geometries/point_xy.hpp>
39 #include <boost/geometry/geometries/polygon.hpp>
44 namespace bg = boost::geometry;
45 typedef bg::model::d2::point_xy<double>
point_xy;
51 namespace NekMeshUtils
55 vector<CADSystem::EdgeLoopSharedPtr> &ein)
60 vector<vector<Array<OneD, NekDouble> > > loopt;
61 for (
int i = 0; i < ein.size(); i++)
63 vector<Array<OneD, NekDouble> > loop;
64 for (
int j = 0; j < ein[i]->edges.size(); j++)
67 NekDouble dt = (bnds[1] - bnds[0]) / (np - 1);
70 for (
int k = 0; k < np - 1; k++)
80 for (
int k = np - 1; k > 0; k--)
89 loopt.push_back(loop);
92 vector<bg::model::polygon<point_xy, false, true> > polygons;
94 for (
int i = 0; i < loopt.size(); i++)
96 bg::model::polygon<point_xy, false, true> polygon;
97 vector<point_xy> points;
98 for (
int j = 0; j < loopt[i].size(); j++)
100 points.push_back(
point_xy(loopt[i][j][0], loopt[i][j][1]));
103 points.push_back(
point_xy(loopt[i][0][0], loopt[i][0][1]));
105 bg::assign_points(polygon, points);
112 bg::centroid(polygon, cen);
115 ein[i]->center[0] = cen.x();
116 ein[i]->center[1] = cen.y();
118 polygons.push_back(polygon);
126 for (
int i = 0; i < ein.size() - 1; i++)
128 if (fabs(ein[i]->area) < fabs(ein[i + 1]->area))
131 swap(ein[i], ein[i + 1]);
132 swap(loopt[i], loopt[i + 1]);
133 swap(polygons[i], polygons[i + 1]);
141 for (
int i = 1; i < ein.size(); i++)
143 point_xy p(ein[i]->center[0], ein[i]->center[1]);
145 if (!bg::within(p, polygons[i]))
151 NekDouble mag = sqrt((n1[0] - n2[0]) * (n1[0] - n2[0]) +
152 (n1[1] - n2[1]) * (n1[1] - n2[1]));
154 N[0] = (n2[1] - n1[1]) / mag;
155 N[1] = -1.0 * (n2[0] - n1[0]) / mag;
165 ASSERTL0(boost::geometry::within(p, polygons[i]),
"point is not side loop");
#define ASSERTL0(condition, msg)
bg::model::d2::point_xy< double > point_xy
boost::shared_ptr< CADSurf > CADSurfSharedPtr