Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpandMeshByRotation.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpandMeshByRotation.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: Take a mesh and expand it by rotating it 180 degrees around x-width/2, 0
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 
41 #include <cstdio>
42 #include <cstdlib>
43 #include <sstream>
44 #include <map>
45 #include <iomanip>
46 #include <tinyxml.h>
47 
48 using namespace std;
49 using namespace Nektar;
50 
51 void ExpandVertices(TiXmlElement* mesh, map<int,int> jointVerts, map<int,int> &fullVerts);
52 
53 void ExpandEdges(TiXmlElement* mesh, map<int,int> &newVerts,
54  map<int,int> jointEdges, map<int,int> &newEdges);
55 
56 void ExpandElmts(TiXmlElement* mesh, map<int,int> &newEdges, int &nOrigElmts);
57 
58 void ExpandComposites(TiXmlElement* mesh, map<int,int> fullEdges, int nOrigElmts);
59 
60 int main(int argc, char *argv[])
61 {
62  int i;
63 
64  if(argc !=3)
65  {
66  fprintf(stderr,"Usage: ./ExpandMeshByRotation meshfile outfile\n");
67  exit(1);
68  }
69 
70  //------------------------------------------------------------
71  // Create Session file which also reads meshfile.
73  = LibUtilities::SessionReader::CreateInstance(argc-1, argv);
74  //-----------------------------------------------------------
75 
76  //-------------------------------------------------------------
77  // Read in mesh from input file
78  string meshfile(argv[argc-2]);
79  SpatialDomains::MeshGraphSharedPtr graphShPt=SpatialDomains::MeshGraph::Read(vSession);
80  SpatialDomains::Composite composite;
81  composite = graphShPt->GetComposite(300);
82  std::map<int,int> jointEdges, jointVerts, newVerts, newEdges;
83  int compsize = composite->size();
84  for(i = 0; i < compsize; ++i)
85  {
86  SpatialDomains::Geometry1DSharedPtr tmp1 = boost::dynamic_pointer_cast<SpatialDomains::Geometry1D>((*composite)[i]);
87  SpatialDomains::Geometry1DSharedPtr tmp2 = boost::dynamic_pointer_cast<SpatialDomains::Geometry1D>((*composite)[compsize-1-i]);
88  jointEdges[tmp1->GetEid() ] = tmp2->GetEid();
89  jointVerts[tmp1->GetVid(0)] = tmp2->GetVid(1);
90  jointVerts[tmp1->GetVid(1)] = tmp2->GetVid(0);
91  }
92 
93 
94  //------------------------------------------------------------
95  TiXmlDocument& meshdoc = vSession->GetDocument();
96 
97  TiXmlHandle docHandle(&meshdoc);
98  TiXmlElement* mesh = docHandle.FirstChildElement("NEKTAR").FirstChildElement("GEOMETRY").Element();
99 
100  int nOrigElmts;
101  //------------------------------------------------------------
102  // Expand Mesh
103  ExpandVertices(mesh,jointVerts,newVerts);
104 
105  ExpandEdges(mesh,newVerts,jointEdges,newEdges);
106 
107  ExpandElmts(mesh, newEdges,nOrigElmts);
108 
109  ExpandComposites( mesh, newEdges, nOrigElmts);
110 
111  meshdoc.SaveFile(argv[argc-1]);
112 
113 }
114 
115 void ExpandVertices(TiXmlElement* mesh, map<int,int> jointVerts, map<int,int> &newVerts)
116 {
117 
118  TiXmlElement* element = mesh->FirstChildElement("VERTEX");
119  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
120 
121  TiXmlElement *vertex = element->FirstChildElement("V");
122 
123  int indx;
124  int nextVertexNumber = -1;
125  int err; /// Error value returned by TinyXML.
126 
127  vector<NekDouble> xpts,ypts,zpts;
128  NekDouble xval, yval, zval;
129 
130  NekDouble yoffset = 0.0;
131  while (vertex)
132  {
133  nextVertexNumber++;
134 
135  TiXmlAttribute *vertexAttr = vertex->FirstAttribute();
136  std::string attrName(vertexAttr->Name());
137 
138  ASSERTL0(attrName == "ID", (std::string("Unknown attribute name: ") + attrName).c_str());
139 
140  err = vertexAttr->QueryIntValue(&indx);
141  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
142 
143  // Now read body of vertex
144  std::string vertexBodyStr;
145 
146  TiXmlNode *vertexBody = vertex->FirstChild();
147 
148  while (vertexBody)
149  {
150  // Accumulate all non-comment body data.
151  if (vertexBody->Type() == TiXmlNode::TEXT)
152  {
153  vertexBodyStr += vertexBody->ToText()->Value();
154  vertexBodyStr += " ";
155  }
156 
157  vertexBody = vertexBody->NextSibling();
158  }
159 
160  ASSERTL0(!vertexBodyStr.empty(), "Vertex definitions must contain vertex data.");
161 
162  // Get vertex data from the data string.
163  std::istringstream vertexDataStrm(vertexBodyStr.c_str());
164 
165  try
166  {
167  while(!vertexDataStrm.fail())
168  {
169  vertexDataStrm >> xval >> yval >> zval;
170  }
171  xpts.push_back(xval);
172  ypts.push_back(yval +yoffset);
173  zpts.push_back(zval);
174  }
175  catch(...)
176  {
177  ASSERTL0(false, "Unable to read VERTEX data.");
178  }
179  vertex = vertex->NextSiblingElement("V");
180  }
181 
182  // Add in newvertices
183  int npts = xpts.size();
184  NekDouble xmax = Vmath::Vmax(npts,&xpts[0],1);
185  int cnt = npts;
186  int cnt1 = npts;
187 
188  for(int i = 0; i < npts; ++i)
189  {
190  if(jointVerts.count(i) == 0)
191  {
192  stringstream s;
193  xval = xmax - xpts[i];
194  yval = -ypts[i];
195  s << scientific << setprecision(8)
196  << xval << " " << yval << " " << zpts[i];
197  TiXmlElement * v = new TiXmlElement( "V" );
198  v->SetAttribute("ID",cnt1);
199  v->LinkEndChild(new TiXmlText(s.str()));
200  element->LinkEndChild(v);
201  newVerts[cnt++] = cnt1++;
202  }
203  else
204  {
205  newVerts[cnt++] = jointVerts[i];
206  }
207  }
208 }
209 
210 void ExpandEdges(TiXmlElement* mesh, map<int,int> &newVerts, map<int,int> jointEdges, map<int,int> &newEdges)
211 {
212  /// Look for elements in ELEMENT block.
213  TiXmlElement* field = mesh->FirstChildElement("EDGE");
214 
215  ASSERTL0(field, "Unable to find EDGE tag in file.");
216 
217  /// All elements are of the form: "<E ID="#"> ... </E>", with
218  /// ? being the element type.
219  /// Read the ID field first.
220  TiXmlElement *edge = field->FirstChildElement("E");
221 
222  /// Since all edge data is one big text block, we need to
223  /// accumulate all TEXT data and then parse it. This
224  /// approach effectively skips all comments or other node
225  /// types since we only care about the edge list. We
226  /// cannot handle missing edge numbers as we could with
227  /// missing element numbers due to the text block format.
228  std::string edgeStr;
229  int indx;
230  int nextEdgeNumber = -1;
231 
232  map<int,int> edgeVert0,edgeVert1;
233  while(edge)
234  {
235  nextEdgeNumber++;
236 
237  int err = edge->QueryIntAttribute("ID",&indx);
238  ASSERTL0(err == TIXML_SUCCESS, "Unable to read edge attribute ID.");
239 
240  TiXmlNode *child = edge->FirstChild();
241  edgeStr.clear();
242  if (child->Type() == TiXmlNode::TEXT)
243  {
244  edgeStr += child->ToText()->ValueStr();
245  }
246 
247  /// Now parse out the edges, three fields at a time.
248  int vertex1, vertex2;
249  std::istringstream edgeDataStrm(edgeStr.c_str());
250 
251  try
252  {
253  while (!edgeDataStrm.fail())
254  {
255  edgeDataStrm >> vertex1 >> vertex2;
256  }
257  edgeVert0[indx] = vertex1;
258  edgeVert1[indx] = vertex2;
259  }
260  catch(...)
261  {
262  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read edge data: ") + edgeStr).c_str());
263  }
264  edge = edge->NextSiblingElement("E");
265  }
266 
267  int nedges = edgeVert0.size();
268  int cnt = nedges;
269  int cnt1 = nedges;
270  int norigverts = newVerts.size();
271 
272  for(int i = 0; i < nedges; ++i)
273  {
274  if(jointEdges.count(i) == 0)
275  {
276  stringstream s;
277 
278  s << setw(5) << newVerts[edgeVert0[i]+norigverts] << " "
279  << newVerts[edgeVert1[i]+norigverts] << " ";
280  TiXmlElement * e = new TiXmlElement( "E" );
281  e->SetAttribute("ID",cnt1);
282  e->LinkEndChild( new TiXmlText(s.str()) );
283  field->LinkEndChild(e);
284  newEdges[cnt++] = cnt1++;
285  }
286  else
287  {
288  newEdges[cnt++] = jointEdges[i];
289  }
290  }
291 }
292 
293 void ExpandElmts(TiXmlElement* mesh, map<int,int> &newEdges, int &nelmts)
294 {
295 
296  TiXmlElement* field = mesh->FirstChildElement("ELEMENT");
297 
298  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
299 
300  int nextElementNumber = -1;
301 
302  /// All elements are of the form: "<? ID="#"> ... </?>", with
303  /// ? being the element type.
304 
305  TiXmlElement *element = field->FirstChildElement();
306 
307  map<int,vector<int> > ElmtEdges;
308 
309  while (element)
310  {
311  std::string elementType(element->ValueStr());
312 
313  ASSERTL0(elementType == "Q" || elementType == "T",
314  (std::string("Unknown 2D element type: ") + elementType).c_str());
315 
316  /// These should be ordered.
317  nextElementNumber++;
318 
319  /// Read id attribute.
320  int indx;
321  int err = element->QueryIntAttribute("ID", &indx);
322  ASSERTL0(err == TIXML_SUCCESS, "Unable to read element attribute ID.");
323 
324  /// Read text element description.
325  TiXmlNode* elementChild = element->FirstChild();
326  std::string elementStr;
327  while(elementChild)
328  {
329  if (elementChild->Type() == TiXmlNode::TEXT)
330  {
331  elementStr += elementChild->ToText()->ValueStr();
332  }
333  elementChild = elementChild->NextSibling();
334  }
335 
336  ASSERTL0(!elementStr.empty(), "Unable to read element description body.");
337 
338  /// Parse out the element components corresponding to type of element.
339  if (elementType == "T")
340  {
341  // Read three edge numbers
342  int edge1, edge2, edge3;
343  std::istringstream elementDataStrm(elementStr.c_str());
344 
345  try
346  {
347  elementDataStrm >> edge1;
348  elementDataStrm >> edge2;
349  elementDataStrm >> edge3;
350 
351  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
352 
353  vector<int> edges;
354  edges.push_back(edge1);
355  edges.push_back(edge2);
356  edges.push_back(edge3);
357 
358  ElmtEdges[indx] = edges;
359  }
360  catch(...)
361  {
362  NEKERROR(ErrorUtil::efatal, (std::string("Unable to read element data for TRIANGLE: ") + elementStr).c_str());
363  }
364  }
365  else if (elementType == "Q")
366  {
367  // Read four edge numbers
368  int edge1, edge2, edge3, edge4;
369  std::istringstream elementDataStrm(elementStr.c_str());
370 
371  try
372  {
373  elementDataStrm >> edge1;
374  elementDataStrm >> edge2;
375  elementDataStrm >> edge3;
376  elementDataStrm >> edge4;
377 
378  ASSERTL0(!elementDataStrm.fail(), (std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
379 
380  vector<int> edges;
381  edges.push_back(edge1);
382  edges.push_back(edge2);
383  edges.push_back(edge3);
384  edges.push_back(edge4);
385 
386  ElmtEdges[indx] = edges;
387 
388  }
389  catch(...)
390  {
391  NEKERROR(ErrorUtil::efatal,(std::string("Unable to read element data for QUAD: ") + elementStr).c_str());
392  }
393  }
394 
395  /// Keep looking
396  element = element->NextSiblingElement();
397  }
398 
399  nelmts = ElmtEdges.size();
400  int nedges = newEdges.size();
401 
402  for(int i = 0; i < ElmtEdges.size(); ++i)
403  {
404  stringstream s;
405 
406  for(int j = 0; j < ElmtEdges[i].size() ; ++j)
407  {
408  s << setw(10) << newEdges[ElmtEdges[i][j]+nedges];
409  }
410 
411  TiXmlElement * f;
412  switch(ElmtEdges[i].size())
413  {
414  case 3:
415  f = new TiXmlElement("T");
416  break;
417  case 4:
418  f = new TiXmlElement("Q");
419  break;
420  default:
421  abort();
422  }
423  f->SetAttribute("ID", i+nelmts);
424  f->LinkEndChild( new TiXmlText(s.str()));
425  field->LinkEndChild(f);
426  }
427 }
428 
429 // Generate a Nektar++ string describing the vector of ints
430 string GetXmlString(char tag, vector<unsigned int> &ids)
431 {
432  stringstream st;
434  bool range = false;
435  int vId = ids[0];
436  int prevId = vId;
437 
438  st << " " << tag << "[" << vId;
439 
440  for (it = ids.begin()+1; it != ids.end(); ++it){
441  // store previous element ID and get current one
442  prevId = vId;
443  vId = (*it);
444 
445  // continue an already started range
446  if (prevId > -1 && vId == prevId + 1)
447  {
448  range = true;
449  // if this is the last element, it's the end of a range, so write
450  if (*it == ids.back())
451  {
452  st << "-" << vId;
453  }
454  continue;
455  }
456 
457  // terminate a range, if present
458  if (range)
459  {
460  st << "-" << prevId;
461  range = false;
462  }
463 
464  // write what will be either a single entry or start of new range
465  st << "," << vId;
466  }
467  // terminate
468  st << "] ";
469  return st.str();
470 }
471 
472 
473 void ExpandComposites(TiXmlElement * mesh, map<int,int> newEdges, int nOrigElmts)
474 {
475  TiXmlElement* field = mesh->FirstChildElement("COMPOSITE");
476  ASSERTL0(field, "Unable to find COMPOSITE tag in file.");
477 
478  int nextCompositeNumber = -1;
479 
480  /// All elements are of the form: "<C ID = "N"> ... </C>".
481  /// Read the ID field first.
482  TiXmlElement *composite = field->FirstChildElement("C");
483 
484  while (composite)
485  {
486  nextCompositeNumber++;
487 
488  int indx;
489  int err = composite->QueryIntAttribute("ID", &indx);
490  ASSERTL0(err == TIXML_SUCCESS, "Unable to read attribute ID.");
491 
492  TiXmlNode* compositeChild = composite->FirstChild();
493  // This is primarily to skip comments that may be present.
494  // Comments appear as nodes just like elements.
495  // We are specifically looking for text in the body
496  // of the definition.
497  while(compositeChild && compositeChild->Type() != TiXmlNode::TEXT)
498  {
499  compositeChild = compositeChild->NextSibling();
500  }
501 
502  ASSERTL0(compositeChild, "Unable to read composite definition body.");
503  std::string compositeStr = compositeChild->ToText()->ValueStr();
504 
505  /// Parse out the element components corresponding to type of element.
506  std::istringstream compositeDataStrm(compositeStr.c_str());
507 
508 
509  try
510  {
511 
512  std::string compositeElementStr;
513  compositeDataStrm >> compositeElementStr;
514 
515  std::istringstream tokenStream(compositeElementStr);
516  char type;
517 
518  tokenStream >> type;
519 
520  // in what follows we are assuming there is only one block of data
521  std::string::size_type indxBeg = compositeElementStr.find_first_of('[') + 1;
522  std::string::size_type indxEnd = compositeElementStr.find_last_of(']') - 1;
523 
524  ASSERTL0(indxBeg <= indxEnd, (std::string("Error reading index definition:") + compositeElementStr).c_str());
525 
526  std::string indxStr = compositeElementStr.substr(indxBeg, indxEnd - indxBeg + 1);
527  std::vector<unsigned int> seqVector;
529 
530  bool err = ParseUtils::GenerateSeqVector(indxStr.c_str(), seqVector);
531 
532  ASSERTL0(err, (std::string("Error reading composite elements: ") + indxStr).c_str());
533 
534  switch(type)
535  {
536  case 'E': // Expand edges using newEdges map and known values
537  {
538  int seqlen = seqVector.size();
539  int nedges = newEdges.size();
540 
541  map<unsigned int, unsigned int> seqMap;
542 
543  for(int i =0; i < seqlen; ++i) // set up a map of defined edges
544  {
545  seqMap[seqVector[i]] = 1;
546  }
547 
548  // if edge does not exist in composite add it to the list
549  for(int i =0; i < seqlen; ++i)
550  {
551  if(seqMap.count(newEdges[seqVector[i]+nedges]) == 0)
552  {
553  seqVector.push_back(newEdges[seqVector[i]+nedges]);
554  }
555  }
556  }
557  break;
558 
559  case 'T': case 'Q': // Expand Triangles & Quads with new elements
560  {
561  int seqlen = seqVector.size();
562 
563  for(int i = 0; i < seqlen; ++i)
564  {
565  seqVector.push_back(seqVector[i]+nOrigElmts);
566  }
567 
568  break;
569  }
570  default:
571  NEKERROR(ErrorUtil::efatal, (std::string("Unrecognized composite token: ") + compositeElementStr).c_str());
572  }
573 
574  // now redefine string in composite
575 
576 
577  composite->ReplaceChild(compositeChild, TiXmlText(GetXmlString(type,seqVector)));
578 
579 
580  }
581  catch(...)
582  {
584  (std::string("Unable to read COMPOSITE data for composite: ") + compositeStr).c_str());
585  }
586 
587  /// Keep looking
588  composite = composite->NextSiblingElement("C");
589  }
590 }
591 
592