Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessInsertSurface.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessInsertSurface.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Calculate Jacobians of elements.
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 #include "ProcessInsertSurface.h"
38 
39 #include <boost/geometry.hpp>
40 #include <boost/geometry/geometries/point.hpp>
41 #include <boost/geometry/geometries/box.hpp>
42 #include <boost/geometry/index/rtree.hpp>
43 
44 namespace bg = boost::geometry;
45 namespace bgi = boost::geometry::index;
46 
47 using namespace std;
48 using namespace Nektar::NekMeshUtils;
49 
50 namespace Nektar
51 {
52 namespace Utilities
53 {
54 
55 ModuleKey ProcessInsertSurface::className = GetModuleFactory().RegisterCreatorFunction(
56  ModuleKey(eProcessModule, "insertsurface"),
57  ProcessInsertSurface::create,
58  "Insert high-order surface mesh into current working mesh.");
59 
60 ProcessInsertSurface::ProcessInsertSurface(MeshSharedPtr m) : ProcessModule(m)
61 {
62  m_config["mesh"] =
63  ConfigOption(false, "", "Mesh to be inserted.");
64  m_config["nonconforming"] =
65  ConfigOption(false,"", "Relax tests for nonconforming boundries");
66 }
67 
69 {
70 }
71 
73 {
74  typedef bg::model::point<NekDouble, 3, bg::cs::cartesian> Point;
75  typedef pair<Point, unsigned int> PointI;
76 
77  if (m_mesh->m_verbose)
78  {
79  cout << "ProcessInsertSurface: Inserting mesh... " << endl;
80  }
81 
82  string file = m_config["mesh"].as<string>();
83  bool nonconform = m_config["nonconforming"].beenSet;
84 
85  if (m_mesh->m_verbose)
86  {
87  cout << "inserting surface from " << file << endl;
88  }
89  MeshSharedPtr inMsh = boost::shared_ptr<Mesh>(new Mesh());
90  inMsh->m_verbose = m_mesh->m_verbose;
92  ModuleKey(eInputModule, "xml"), inMsh);
93  mod->RegisterConfig("infile", file);
94  mod->Process();
95 
96  //build ann tree of surface verticies from inMsh
97  //match surface vertices in ccm mesh to inMsh and copy information
98 
99  //tolerance of matching vertices
100  NekDouble tol = 1e-5;
101 
102  NodeSet surfaceNodes;
103  for(int i = 0; i < inMsh->m_element[2].size(); i++)
104  {
105  vector<NodeSharedPtr> ns = inMsh->m_element[2][i]->GetVertexList();
106  for(int j = 0; j < ns.size(); j++)
107  {
108  surfaceNodes.insert(ns[j]);
109  }
110  }
111 
112  vector<NodeSharedPtr> inMshnodeList(surfaceNodes.begin(), surfaceNodes.end());
113 
114  vector<PointI> dataPts;
115  for(int i = 0; i < inMshnodeList.size(); i++)
116  {
117  dataPts.push_back(make_pair(Point( inMshnodeList[i]->m_x,
118  inMshnodeList[i]->m_y,
119  inMshnodeList[i]->m_z), i));
120  }
121 
122  //Build tree
123  bgi::rtree<PointI, bgi::rstar<16> > rtree;
124  rtree.insert(dataPts.begin(), dataPts.end());
125 
126  surfaceNodes.clear();
127  for(int i = 0; i < m_mesh->m_element[2].size(); i++)
128  {
129  vector<NodeSharedPtr> ns = m_mesh->m_element[2][i]->GetVertexList();
130  for(int j = 0; j < ns.size(); j++)
131  {
132  surfaceNodes.insert(ns[j]);
133  }
134  }
135 
136  if(!nonconform)
137  {
138  ASSERTL0(surfaceNodes.size() == inMshnodeList.size(),
139  "surface mesh node count mismatch, will not work");
140  }
141 
142  EdgeSet surfEdges;
143  for(int i = 0; i < m_mesh->m_element[2].size(); i++)
144  {
145  FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
146  vector<EdgeSharedPtr> es = f->m_edgeList;
147  for(int j = 0; j < es.size(); j++)
148  {
149  surfEdges.insert(es[j]);
150  }
151  }
152 
154  for(it = surfEdges.begin(); it != surfEdges.end(); it++)
155  {
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));
159 
160  NekDouble dist1 = bg::distance(result[0].first, queryPt1);
161  if(nonconform)
162  {
163  if(dist1 > tol)
164  {
165  continue;
166  }
167  }
168  else
169  {
170  ASSERTL0(dist1 < tol, "cannot locate point accurately enough");
171  }
172 
173  NodeSharedPtr inN1 = inMshnodeList[result[0].second];
174 
175  Point queryPt2((*it)->m_n2->m_x, (*it)->m_n2->m_y, (*it)->m_n2->m_z);
176  result.clear();
177  rtree.query(bgi::nearest(queryPt2, 1), std::back_inserter(result));
178 
179  NekDouble dist2 = bg::distance(result[0].first, queryPt2);
180  if(nonconform)
181  {
182  if(dist2 > tol)
183  {
184  continue;
185  }
186  }
187  else
188  {
189  ASSERTL0(dist2 < tol, "cannot locate point accurately enough");
190  }
191  NodeSharedPtr inN2 = inMshnodeList[result[0].second];
192 
193  EdgeSharedPtr tst = boost::shared_ptr<Edge>(new Edge(inN1,inN2));
194 
195  EdgeSet::iterator f = inMsh->m_edgeSet.find(tst);
196 
197  ASSERTL0(f != inMsh->m_edgeSet.end(),"could not find edge in input");
198 
199  (*it)->m_edgeNodes = (*f)->m_edgeNodes;
200  (*it)->m_curveType = (*f)->m_curveType;
201 
202  if((*f)->m_n1->Distance((*it)->m_n1) > tol)
203  {
204  reverse((*it)->m_edgeNodes.begin(),(*it)->m_edgeNodes.end());
205  }
206  }
207 }
208 }
209 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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.
Definition: NekFactory.hpp:162
Represents an edge which joins two points.
Definition: Edge.h:58
STL namespace.
pair< ModuleType, string > ModuleKey
virtual void Process()
Write mesh to output file.
boost::unordered_set< NodeSharedPtr, NodeHash > NodeSet
Definition: Node.h:441
A 0-dimensional vertex.
Definition: Point.h:49
boost::shared_ptr< Module > ModuleSharedPtr
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:135
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
Abstract base class for processing modules.
boost::unordered_set< EdgeSharedPtr, EdgeHash > EdgeSet
Definition: Edge.h:161
boost::shared_ptr< Face > FaceSharedPtr
Definition: Face.h:148
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215