36 #include <boost/geometry.hpp>
37 #include <boost/geometry/geometries/point.hpp>
38 #include <boost/geometry/geometries/box.hpp>
39 #include <boost/geometry/index/rtree.hpp>
44 namespace bg = boost::geometry;
45 namespace bgi = boost::geometry::index;
57 ProcessInsertSurface::create,
58 "Insert high-order surface mesh into current working mesh.");
65 ConfigOption(
false,
"",
"Relax tests for nonconforming boundries");
74 typedef bg::model::point<NekDouble, 3, bg::cs::cartesian>
Point;
75 typedef pair<Point, unsigned int> PointI;
79 cout <<
"ProcessInsertSurface: Inserting mesh... " << endl;
82 string file =
m_config[
"mesh"].as<
string>();
83 bool nonconform =
m_config[
"nonconforming"].beenSet;
87 cout <<
"inserting surface from " << file << endl;
90 inMsh->m_verbose =
m_mesh->m_verbose;
93 mod->RegisterConfig(
"infile", file);
103 for(
int i = 0; i < inMsh->m_element[2].size(); i++)
105 vector<NodeSharedPtr> ns = inMsh->m_element[2][i]->GetVertexList();
106 for(
int j = 0; j < ns.size(); j++)
108 surfaceNodes.insert(ns[j]);
112 vector<NodeSharedPtr> inMshnodeList(surfaceNodes.begin(), surfaceNodes.end());
114 vector<PointI> dataPts;
115 for(
int i = 0; i < inMshnodeList.size(); i++)
117 dataPts.push_back(make_pair(Point( inMshnodeList[i]->m_x,
118 inMshnodeList[i]->m_y,
119 inMshnodeList[i]->m_z), i));
123 bgi::rtree<PointI, bgi::rstar<16> > rtree;
124 rtree.insert(dataPts.begin(), dataPts.end());
126 surfaceNodes.clear();
127 for(
int i = 0; i <
m_mesh->m_element[2].size(); i++)
129 vector<NodeSharedPtr> ns =
m_mesh->m_element[2][i]->GetVertexList();
130 for(
int j = 0; j < ns.size(); j++)
132 surfaceNodes.insert(ns[j]);
138 ASSERTL0(surfaceNodes.size() == inMshnodeList.size(),
139 "surface mesh node count mismatch, will not work");
143 for(
int i = 0; i <
m_mesh->m_element[2].size(); i++)
146 vector<EdgeSharedPtr> es = f->m_edgeList;
147 for(
int j = 0; j < es.size(); j++)
149 surfEdges.insert(es[j]);
154 for(it = surfEdges.begin(); it != surfEdges.end(); it++)
156 Point queryPt1((*it)->m_n1->m_x, (*it)->m_n1->m_y, (*it)->m_n1->m_z);
157 vector<PointI> result;
158 rtree.query(bgi::nearest(queryPt1, 1), std::back_inserter(result));
160 NekDouble dist1 = bg::distance(result[0].first, queryPt1);
170 ASSERTL0(dist1 < tol,
"cannot locate point accurately enough");
175 Point queryPt2((*it)->m_n2->m_x, (*it)->m_n2->m_y, (*it)->m_n2->m_z);
177 rtree.query(bgi::nearest(queryPt2, 1), std::back_inserter(result));
179 NekDouble dist2 = bg::distance(result[0].first, queryPt2);
189 ASSERTL0(dist2 < tol,
"cannot locate point accurately enough");
197 ASSERTL0(f != inMsh->m_edgeSet.end(),
"could not find edge in input");
199 (*it)->m_edgeNodes = (*f)->m_edgeNodes;
200 (*it)->m_curveType = (*f)->m_curveType;
202 if((*f)->m_n1->Distance((*it)->m_n1) > tol)
204 reverse((*it)->m_edgeNodes.begin(),(*it)->m_edgeNodes.end());
#define ASSERTL0(condition, msg)
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Represents an edge which joins two points.
pair< ModuleType, string > ModuleKey
MeshSharedPtr m_mesh
Mesh object.
virtual void Process()
Write mesh to output file.
boost::unordered_set< NodeSharedPtr, NodeHash > NodeSet
boost::shared_ptr< Module > ModuleSharedPtr
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
std::map< std::string, ConfigOption > m_config
List of configuration values.
virtual ~ProcessInsertSurface()
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Abstract base class for processing modules.
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
boost::shared_ptr< Face > FaceSharedPtr
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.