Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
OutputNekpp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: OutputNekpp.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: Nektar++ file format output.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <set>
37 #include <string>
38 using namespace std;
39 
40 #include <boost/iostreams/filtering_streambuf.hpp>
41 #include <boost/iostreams/copy.hpp>
42 #include <boost/iostreams/filter/gzip.hpp>
43 namespace io = boost::iostreams;
44 
45 #include <tinyxml.h>
49 
51 #include "OutputNekpp.h"
52 
53 using namespace Nektar::NekMeshUtils;
54 
55 namespace Nektar
56 {
57 namespace Utilities
58 {
59 ModuleKey OutputNekpp::className =
61  OutputNekpp::create,
62  "Writes a Nektar++ xml file.");
63 
64 OutputNekpp::OutputNekpp(MeshSharedPtr m) : OutputModule(m)
65 {
66  m_config["z"] = ConfigOption(
67  true, "0", "Compress output file and append a .gz extension.");
68  m_config["test"] = ConfigOption(
69  true, "0", "Attempt to load resulting mesh and create meshgraph.");
70  m_config["uncompress"] = ConfigOption(true, "0", "Uncompress xml sections");
71 }
72 
74 {
75 }
76 
77 template <typename T> void TestElmts(SpatialDomains::MeshGraphSharedPtr &graph)
78 {
79  const std::map<int, boost::shared_ptr<T> > &tmp =
80  graph->GetAllElementsOfType<T>();
81  typename std::map<int, boost::shared_ptr<T> >::const_iterator it1, it2;
82 
83  SpatialDomains::CurveMap &curvedEdges = graph->GetCurvedEdges();
84  SpatialDomains::CurveMap &curvedFaces = graph->GetCurvedFaces();
85 
86  for (it1 = tmp.begin(), it2 = tmp.end(); it1 != it2; ++it1)
87  {
88  SpatialDomains::GeometrySharedPtr geom = it1->second;
89  geom->FillGeom();
90  geom->Reset(curvedEdges, curvedFaces);
91  }
92 }
93 
95 {
96  if (m_mesh->m_verbose)
97  {
98  cout << "OutputNekpp: Writing file..." << endl;
99  }
100 
101  TiXmlDocument doc;
102  TiXmlDeclaration *decl = new TiXmlDeclaration("1.0", "utf-8", "");
103  doc.LinkEndChild(decl);
104 
105  TiXmlElement *root = new TiXmlElement("NEKTAR");
106  doc.LinkEndChild(root);
107 
108  // Begin <GEOMETRY> section
109  TiXmlElement *geomTag = new TiXmlElement("GEOMETRY");
110  geomTag->SetAttribute("DIM", m_mesh->m_expDim);
111  geomTag->SetAttribute("SPACE", m_mesh->m_spaceDim);
112  root->LinkEndChild(geomTag);
113 
114  WriteXmlNodes(geomTag);
115  WriteXmlEdges(geomTag);
116  WriteXmlFaces(geomTag);
117  WriteXmlElements(geomTag);
118  WriteXmlCurves(geomTag);
119  WriteXmlComposites(geomTag);
120  WriteXmlDomain(geomTag);
121  WriteXmlExpansions(root);
122  WriteXmlConditions(root);
123 
124  // Extract the output filename and extension
125  string filename = m_config["outfile"].as<string>();
126 
127  // Compress output and append .gz extension
128  if (m_config["z"].as<bool>())
129  {
130  filename += ".gz";
131  ofstream fout(filename.c_str(),
132  std::ios_base::out | std::ios_base::binary);
133 
134  std::stringstream decompressed;
135  decompressed << doc;
136  io::filtering_streambuf<io::output> out;
137  out.push(io::gzip_compressor());
138  out.push(fout);
139  io::copy(decompressed, out);
140 
141  fout.close();
142  }
143  else
144  {
145  doc.SaveFile(filename);
146  }
147 
148  // Test the resulting XML file (with a basic test) by loading it
149  // with the session reader, generating the MeshGraph and testing if
150  // each element is valid.
151  if (m_config["test"].beenSet)
152  {
153  vector<string> filenames(1);
154  filenames[0] = filename;
155 
160 
161  TestElmts<SpatialDomains::SegGeom>(graphShPt);
162  TestElmts<SpatialDomains::TriGeom>(graphShPt);
163  TestElmts<SpatialDomains::QuadGeom>(graphShPt);
164  TestElmts<SpatialDomains::TetGeom>(graphShPt);
165  TestElmts<SpatialDomains::PrismGeom>(graphShPt);
166  TestElmts<SpatialDomains::PyrGeom>(graphShPt);
167  TestElmts<SpatialDomains::HexGeom>(graphShPt);
168  }
169 }
170 
171 void OutputNekpp::WriteXmlNodes(TiXmlElement *pRoot)
172 {
173  bool UnCompressed = m_config["uncompress"].as<bool>();
174 
175  TiXmlElement *verTag = new TiXmlElement("VERTEX");
177 
178  std::set<NodeSharedPtr> tmp(m_mesh->m_vertexSet.begin(),
179  m_mesh->m_vertexSet.end());
180 
181  if (UnCompressed)
182  {
183  for (it = tmp.begin(); it != tmp.end(); ++it)
184  {
185  NodeSharedPtr n = *it;
186  stringstream s;
187  s << scientific << setprecision(8) << n->m_x << " " << n->m_y << " "
188  << n->m_z;
189  TiXmlElement *v = new TiXmlElement("V");
190  v->SetAttribute("ID", n->m_id);
191  v->LinkEndChild(new TiXmlText(s.str()));
192  verTag->LinkEndChild(v);
193  }
194  }
195  else
196  {
197  std::vector<LibUtilities::MeshVertex> vertInfo;
198  for (it = tmp.begin(); it != tmp.end(); ++it)
199  {
201  NodeSharedPtr n = *it;
202  v.id = n->m_id;
203  v.x = n->m_x;
204  v.y = n->m_y;
205  v.z = n->m_z;
206  vertInfo.push_back(v);
207  }
208  std::string vertStr;
210  verTag->SetAttribute("COMPRESSED",
212  verTag->SetAttribute("BITSIZE",
214 
215  verTag->LinkEndChild(new TiXmlText(vertStr));
216  }
217 
218  pRoot->LinkEndChild(verTag);
219 }
220 
221 void OutputNekpp::WriteXmlEdges(TiXmlElement *pRoot)
222 {
223  bool UnCompressed = m_config["uncompress"].as<bool>();
224 
225  if (m_mesh->m_expDim >= 2)
226  {
227  TiXmlElement *verTag = new TiXmlElement("EDGE");
228 
230  std::set<EdgeSharedPtr> tmp(m_mesh->m_edgeSet.begin(),
231  m_mesh->m_edgeSet.end());
232  if (UnCompressed)
233  {
234  for (it = tmp.begin(); it != tmp.end(); ++it)
235  {
236  EdgeSharedPtr ed = *it;
237  stringstream s;
238 
239  s << setw(5) << ed->m_n1->m_id << " " << ed->m_n2->m_id
240  << " ";
241  TiXmlElement *e = new TiXmlElement("E");
242  e->SetAttribute("ID", ed->m_id);
243  e->LinkEndChild(new TiXmlText(s.str()));
244  verTag->LinkEndChild(e);
245  }
246  }
247  else
248  {
249  std::vector<LibUtilities::MeshEdge> edgeInfo;
250  for (it = tmp.begin(); it != tmp.end(); ++it)
251  {
253  EdgeSharedPtr ed = *it;
254 
255  e.id = ed->m_id;
256  e.v0 = ed->m_n1->m_id;
257  e.v1 = ed->m_n2->m_id;
258 
259  edgeInfo.push_back(e);
260  }
261  std::string edgeStr;
263  edgeStr);
264  verTag->SetAttribute(
266  verTag->SetAttribute("BITSIZE",
268  verTag->LinkEndChild(new TiXmlText(edgeStr));
269  }
270  pRoot->LinkEndChild(verTag);
271  }
272 }
273 
274 void OutputNekpp::WriteXmlFaces(TiXmlElement *pRoot)
275 {
276  bool UnCompressed = m_config["uncompress"].as<bool>();
277 
278  if (m_mesh->m_expDim == 3)
279  {
280  TiXmlElement *verTag = new TiXmlElement("FACE");
282  std::set<FaceSharedPtr> tmp(m_mesh->m_faceSet.begin(),
283  m_mesh->m_faceSet.end());
284 
285  if (UnCompressed)
286  {
287  for (it = tmp.begin(); it != tmp.end(); ++it)
288  {
289  stringstream s;
290  FaceSharedPtr fa = *it;
291 
292  for (int j = 0; j < fa->m_edgeList.size(); ++j)
293  {
294  s << setw(10) << fa->m_edgeList[j]->m_id;
295  }
296  TiXmlElement *f;
297  switch (fa->m_vertexList.size())
298  {
299  case 3:
300  f = new TiXmlElement("T");
301  break;
302  case 4:
303  f = new TiXmlElement("Q");
304  break;
305  default:
306  abort();
307  }
308  f->SetAttribute("ID", fa->m_id);
309  f->LinkEndChild(new TiXmlText(s.str()));
310  verTag->LinkEndChild(f);
311  }
312  }
313  else
314  {
315  std::vector<LibUtilities::MeshTri> TriFaceInfo;
316  std::vector<LibUtilities::MeshQuad> QuadFaceInfo;
317 
318  for (it = tmp.begin(); it != tmp.end(); ++it)
319  {
320  FaceSharedPtr fa = *it;
321 
322  switch (fa->m_edgeList.size())
323  {
324  case 3:
325  {
327  f.id = fa->m_id;
328  for (int i = 0; i < 3; ++i)
329  {
330  f.e[i] = fa->m_edgeList[i]->m_id;
331  }
332  TriFaceInfo.push_back(f);
333  }
334  break;
335  case 4:
336  {
338  f.id = fa->m_id;
339  for (int i = 0; i < 4; ++i)
340  {
341  f.e[i] = fa->m_edgeList[i]->m_id;
342  }
343  QuadFaceInfo.push_back(f);
344  }
345  break;
346  default:
347  ASSERTL0(false, "Unkonwn face type");
348  }
349  }
350 
351  if (TriFaceInfo.size())
352  {
353  std::string vType("T");
354  TiXmlElement *x = new TiXmlElement(vType);
355  std::string faceStr;
357  faceStr);
358  x->SetAttribute(
359  "COMPRESSED",
361  x->SetAttribute("BITSIZE",
363  x->LinkEndChild(new TiXmlText(faceStr));
364  verTag->LinkEndChild(x);
365  }
366 
367  if (QuadFaceInfo.size())
368  {
369  std::string vType("Q");
370  TiXmlElement *x = new TiXmlElement(vType);
371  std::string faceStr;
373  faceStr);
374  x->SetAttribute(
375  "COMPRESSED",
377  x->SetAttribute("BITSIZE",
379  x->LinkEndChild(new TiXmlText(faceStr));
380  verTag->LinkEndChild(x);
381  }
382  }
383  pRoot->LinkEndChild(verTag);
384  }
385 }
386 
387 void OutputNekpp::WriteXmlElements(TiXmlElement *pRoot)
388 {
389  bool UnCompressed = m_config["uncompress"].as<bool>();
390 
391  TiXmlElement *verTag = new TiXmlElement("ELEMENT");
392  vector<ElementSharedPtr> &elmt = m_mesh->m_element[m_mesh->m_expDim];
393 
394  if (UnCompressed)
395  {
396  for (int i = 0; i < elmt.size(); ++i)
397  {
398  TiXmlElement *elm_tag = new TiXmlElement(elmt[i]->GetTag());
399  elm_tag->SetAttribute("ID", elmt[i]->GetId());
400  elm_tag->LinkEndChild(new TiXmlText(elmt[i]->GetXmlString()));
401  verTag->LinkEndChild(elm_tag);
402  }
403  }
404  else
405  {
406  std::vector<LibUtilities::MeshEdge> SegInfo;
407  std::vector<LibUtilities::MeshTri> TriInfo;
408  std::vector<LibUtilities::MeshQuad> QuadInfo;
409  std::vector<LibUtilities::MeshTet> TetInfo;
410  std::vector<LibUtilities::MeshPyr> PyrInfo;
411  std::vector<LibUtilities::MeshPrism> PrismInfo;
412  std::vector<LibUtilities::MeshHex> HexInfo;
413 
414  for (int i = 0; i < elmt.size(); ++i)
415  {
416  switch (elmt[i]->GetTag()[0])
417  {
418  case 'S':
419  {
421  e.id = elmt[i]->GetId();
422  e.v0 = elmt[i]->GetVertex(0)->m_id;
423  e.v1 = elmt[i]->GetVertex(1)->m_id;
424  SegInfo.push_back(e);
425  }
426  break;
427  case 'T':
428  {
430  e.id = elmt[i]->GetId();
431  for (int j = 0; j < 3; ++j)
432  {
433  e.e[j] = elmt[i]->GetEdge(j)->m_id;
434  }
435  TriInfo.push_back(e);
436  }
437  break;
438  case 'Q':
439  {
441  e.id = elmt[i]->GetId();
442  for (int j = 0; j < 4; ++j)
443  {
444  e.e[j] = elmt[i]->GetEdge(j)->m_id;
445  }
446  QuadInfo.push_back(e);
447  }
448  break;
449  case 'A':
450  {
452  e.id = elmt[i]->GetId();
453  for (int j = 0; j < 4; ++j)
454  {
455  e.f[j] = elmt[i]->GetFace(j)->m_id;
456  }
457  TetInfo.push_back(e);
458  }
459  break;
460  case 'P':
461  {
463  e.id = elmt[i]->GetId();
464  for (int j = 0; j < 5; ++j)
465  {
466  e.f[j] = elmt[i]->GetFace(j)->m_id;
467  }
468  PyrInfo.push_back(e);
469  }
470  break;
471  case 'R':
472  {
474  e.id = elmt[i]->GetId();
475  for (int j = 0; j < 5; ++j)
476  {
477  e.f[j] = elmt[i]->GetFace(j)->m_id;
478  }
479  PrismInfo.push_back(e);
480  }
481  break;
482  case 'H':
483  {
485  e.id = elmt[i]->GetId();
486  for (int j = 0; j < 6; ++j)
487  {
488  e.f[j] = elmt[i]->GetFace(j)->m_id;
489  }
490  HexInfo.push_back(e);
491  }
492  break;
493  default:
494  ASSERTL0(false, "Unknown element type");
495  }
496  }
497 
498  if (SegInfo.size())
499  {
500  std::string vType("S");
501  TiXmlElement *x = new TiXmlElement(vType);
502  std::string Str;
504  x->SetAttribute("COMPRESSED",
506  x->SetAttribute("BITSIZE",
508  x->LinkEndChild(new TiXmlText(Str));
509  verTag->LinkEndChild(x);
510  }
511 
512  if (TriInfo.size())
513  {
514  std::string vType("T");
515  TiXmlElement *x = new TiXmlElement(vType);
516  std::string Str;
518  x->SetAttribute("COMPRESSED",
520  x->SetAttribute("BITSIZE",
522  x->LinkEndChild(new TiXmlText(Str));
523  verTag->LinkEndChild(x);
524  }
525 
526  if (QuadInfo.size())
527  {
528  std::string vType("Q");
529  TiXmlElement *x = new TiXmlElement(vType);
530  std::string Str;
532  x->SetAttribute("COMPRESSED",
534  x->SetAttribute("BITSIZE",
536  x->LinkEndChild(new TiXmlText(Str));
537  verTag->LinkEndChild(x);
538  }
539 
540  if (TetInfo.size())
541  {
542  std::string vType("A");
543  TiXmlElement *x = new TiXmlElement(vType);
544  std::string Str;
546  x->SetAttribute("COMPRESSED",
548  x->SetAttribute("BITSIZE",
550  x->LinkEndChild(new TiXmlText(Str));
551  verTag->LinkEndChild(x);
552  }
553 
554  if (PyrInfo.size())
555  {
556  std::string vType("P");
557  TiXmlElement *x = new TiXmlElement(vType);
558  std::string Str;
560  x->SetAttribute("COMPRESSED",
562  x->SetAttribute("BITSIZE",
564  x->LinkEndChild(new TiXmlText(Str));
565  verTag->LinkEndChild(x);
566  }
567 
568  if (PrismInfo.size())
569  {
570  std::string vType("R");
571  TiXmlElement *x = new TiXmlElement(vType);
572  std::string Str;
574  x->SetAttribute("COMPRESSED",
576  x->SetAttribute("BITSIZE",
578  x->LinkEndChild(new TiXmlText(Str));
579  verTag->LinkEndChild(x);
580  }
581 
582  if (HexInfo.size())
583  {
584  std::string vType("H");
585  TiXmlElement *x = new TiXmlElement(vType);
586  std::string Str;
588  x->SetAttribute("COMPRESSED",
590  x->SetAttribute("BITSIZE",
592  x->LinkEndChild(new TiXmlText(Str));
593  verTag->LinkEndChild(x);
594  }
595  }
596 
597  pRoot->LinkEndChild(verTag);
598 }
599 
600 void OutputNekpp::WriteXmlCurves(TiXmlElement *pRoot)
601 {
602  bool UnCompressed = m_config["uncompress"].as<bool>();
603 
604  int edgecnt = 0;
605 
606  bool curve = false;
608  if (m_mesh->m_expDim > 1)
609  {
610  for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end();
611  ++it)
612  {
613  if ((*it)->m_edgeNodes.size() > 0)
614  {
615  curve = true;
616  break;
617  }
618  }
619  }
620  else if (m_mesh->m_expDim == 1)
621  {
622  for (int i = 0; i < m_mesh->m_element[1].size(); ++i)
623  {
624  if (m_mesh->m_element[1][i]->GetVolumeNodes().size() > 0)
625  {
626  curve = true;
627  break;
628  }
629  }
630  }
631  if (!curve)
632  return;
633 
634  TiXmlElement *curved = new TiXmlElement("CURVED");
635 
636  if (UnCompressed)
637  {
638  for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end();
639  ++it)
640  {
641  if ((*it)->m_edgeNodes.size() > 0)
642  {
643  TiXmlElement *e = new TiXmlElement("E");
644  e->SetAttribute("ID", edgecnt++);
645  e->SetAttribute("EDGEID", (*it)->m_id);
646  e->SetAttribute("NUMPOINTS", (*it)->GetNodeCount());
647  e->SetAttribute(
648  "TYPE", LibUtilities::kPointsTypeStr[(*it)->m_curveType]);
649  TiXmlText *t0 = new TiXmlText((*it)->GetXmlCurveString());
650  e->LinkEndChild(t0);
651  curved->LinkEndChild(e);
652  }
653  }
654 
655  int facecnt = 0;
656 
657  if (m_mesh->m_expDim == 1 && m_mesh->m_spaceDim > 1)
658  {
660  for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
661  it != m_mesh->m_element[m_mesh->m_expDim].end();
662  ++it)
663  {
664  // Only generate face curve if there are volume nodes
665  if ((*it)->GetVolumeNodes().size() > 0)
666  {
667  TiXmlElement *e = new TiXmlElement("E");
668  e->SetAttribute("ID", facecnt++);
669  e->SetAttribute("EDGEID", (*it)->GetId());
670  e->SetAttribute("NUMPOINTS", (*it)->GetNodeCount());
671  e->SetAttribute(
672  "TYPE",
673  LibUtilities::kPointsTypeStr[(*it)->GetCurveType()]);
674 
675  TiXmlText *t0 = new TiXmlText((*it)->GetXmlCurveString());
676  e->LinkEndChild(t0);
677  curved->LinkEndChild(e);
678  }
679  }
680  }
681  // 2D elements in 3-space, output face curvature information
682  else if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
683  {
685  for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
686  it != m_mesh->m_element[m_mesh->m_expDim].end();
687  ++it)
688  {
689  // Only generate face curve if there are volume nodes
690  if ((*it)->GetVolumeNodes().size() > 0)
691  {
692  TiXmlElement *e = new TiXmlElement("F");
693  e->SetAttribute("ID", facecnt++);
694  e->SetAttribute("FACEID", (*it)->GetId());
695  e->SetAttribute("NUMPOINTS", (*it)->GetNodeCount());
696  e->SetAttribute(
697  "TYPE",
698  LibUtilities::kPointsTypeStr[(*it)->GetCurveType()]);
699 
700  TiXmlText *t0 = new TiXmlText((*it)->GetXmlCurveString());
701  e->LinkEndChild(t0);
702  curved->LinkEndChild(e);
703  }
704  }
705  }
706  else if (m_mesh->m_expDim == 3)
707  {
708  FaceSet::iterator it2;
709  for (it2 = m_mesh->m_faceSet.begin();
710  it2 != m_mesh->m_faceSet.end();
711  ++it2)
712  {
713  if ((*it2)->m_faceNodes.size() > 0)
714  {
715  TiXmlElement *f = new TiXmlElement("F");
716  f->SetAttribute("ID", facecnt++);
717  f->SetAttribute("FACEID", (*it2)->m_id);
718  f->SetAttribute("NUMPOINTS", (*it2)->GetNodeCount());
719  f->SetAttribute(
720  "TYPE",
721  LibUtilities::kPointsTypeStr[(*it2)->m_curveType]);
722  TiXmlText *t0 = new TiXmlText((*it2)->GetXmlCurveString());
723  f->LinkEndChild(t0);
724  curved->LinkEndChild(f);
725  }
726  }
727  }
728  }
729  else
730  {
731  std::vector<LibUtilities::MeshCurvedInfo> edgeinfo;
732  std::vector<LibUtilities::MeshCurvedInfo> faceinfo;
733  LibUtilities::MeshCurvedPts curvedpts;
734  curvedpts.id = 0; // assume all points are going in here
735  int ptoffset = 0;
736  int newidx = 0;
737  NodeSet cvertlist;
738 
739  for (it = m_mesh->m_edgeSet.begin(); it != m_mesh->m_edgeSet.end();
740  ++it)
741  {
742  if ((*it)->m_edgeNodes.size() > 0)
743  {
745  cinfo.id = edgecnt++;
746  cinfo.entityid = (*it)->m_id;
747  cinfo.npoints = (*it)->m_edgeNodes.size() + 2;
748  cinfo.ptype = (*it)->m_curveType;
749  cinfo.ptid = 0; // set to just one point set
750  cinfo.ptoffset = ptoffset;
751 
752  edgeinfo.push_back(cinfo);
753 
754  std::vector<NodeSharedPtr> nodeList;
755  (*it)->GetCurvedNodes(nodeList);
756 
757  // fill in points
758  for (int i = 0; i < nodeList.size(); ++i)
759  {
760  pair<NodeSet::iterator, bool> testIns =
761  cvertlist.insert(nodeList[i]);
762 
763  if (testIns.second) // have inserted node
764  {
765  (*(testIns.first))->m_id = newidx;
766 
768  v.id = newidx;
769  v.x = nodeList[i]->m_x;
770  v.y = nodeList[i]->m_y;
771  v.z = nodeList[i]->m_z;
772  curvedpts.pts.push_back(v);
773  newidx++;
774  }
775 
776  curvedpts.index.push_back((*(testIns.first))->m_id);
777  }
778 
779  ptoffset += cinfo.npoints;
780  }
781  }
782 
783  int facecnt = 0;
784 
785  // 1D element in 2 or 3 space
786  if (m_mesh->m_expDim == 1 && m_mesh->m_spaceDim > 1)
787  {
789  for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
790  it != m_mesh->m_element[m_mesh->m_expDim].end();
791  ++it)
792  {
793  // Only generate face curve if there are volume nodes
794  if ((*it)->GetVolumeNodes().size() > 0)
795  {
797  cinfo.id = facecnt++;
798  cinfo.entityid = (*it)->GetId();
799  cinfo.npoints = (*it)->GetNodeCount();
800  cinfo.ptype = (*it)->GetCurveType();
801  cinfo.ptid = 0; // set to just one point set
802  cinfo.ptoffset = ptoffset;
803 
804  edgeinfo.push_back(cinfo);
805 
806  // fill in points
807  vector<NodeSharedPtr> tmp;
808  (*it)->GetCurvedNodes(tmp);
809 
810  for (int i = 0; i < tmp.size(); ++i)
811  {
812  pair<NodeSet::iterator, bool> testIns =
813  cvertlist.insert(tmp[i]);
814 
815  if (testIns.second) // have inserted node
816  {
817  (*(testIns.first))->m_id = newidx;
818 
820  v.id = newidx;
821  v.x = tmp[i]->m_x;
822  v.y = tmp[i]->m_y;
823  v.z = tmp[i]->m_z;
824  curvedpts.pts.push_back(v);
825  newidx++;
826  }
827  curvedpts.index.push_back((*(testIns.first))->m_id);
828  }
829  ptoffset += cinfo.npoints;
830  }
831  }
832  }
833  // 2D elements in 3-space, output face curvature information
834  else if (m_mesh->m_expDim == 2 && m_mesh->m_spaceDim == 3)
835  {
837  for (it = m_mesh->m_element[m_mesh->m_expDim].begin();
838  it != m_mesh->m_element[m_mesh->m_expDim].end();
839  ++it)
840  {
841  // Only generate face curve if there are volume nodes
842  if ((*it)->GetVolumeNodes().size() > 0)
843  {
845  cinfo.id = facecnt++;
846  cinfo.entityid = (*it)->GetId();
847  cinfo.npoints = (*it)->GetNodeCount();
848  cinfo.ptype = (*it)->GetCurveType();
849  cinfo.ptid = 0; // set to just one point set
850  cinfo.ptoffset = ptoffset;
851 
852  faceinfo.push_back(cinfo);
853 
854  // fill in points
855  vector<NodeSharedPtr> tmp;
856  (*it)->GetCurvedNodes(tmp);
857 
858  for (int i = 0; i < tmp.size(); ++i)
859  {
860  pair<NodeSet::iterator, bool> testIns =
861  cvertlist.insert(tmp[i]);
862 
863  if (testIns.second) // have inserted node
864  {
865  (*(testIns.first))->m_id = newidx;
866 
868  v.id = newidx;
869  v.x = tmp[i]->m_x;
870  v.y = tmp[i]->m_y;
871  v.z = tmp[i]->m_z;
872  curvedpts.pts.push_back(v);
873  newidx++;
874  }
875  curvedpts.index.push_back((*(testIns.first))->m_id);
876  }
877  ptoffset += cinfo.npoints;
878  }
879  }
880  }
881  else if (m_mesh->m_expDim == 3)
882  {
883  FaceSet::iterator it2;
884  for (it2 = m_mesh->m_faceSet.begin();
885  it2 != m_mesh->m_faceSet.end();
886  ++it2)
887  {
888  if ((*it2)->m_faceNodes.size() > 0)
889  {
890  vector<NodeSharedPtr> tmp;
891  (*it2)->GetCurvedNodes(tmp);
892 
894  cinfo.id = facecnt++;
895  cinfo.entityid = (*it2)->m_id;
896  cinfo.npoints = tmp.size();
897  cinfo.ptype = (*it2)->m_curveType;
898  cinfo.ptid = 0; // set to just one point set
899  cinfo.ptoffset = ptoffset;
900 
901  faceinfo.push_back(cinfo);
902 
903  for (int i = 0; i < tmp.size(); ++i)
904  {
905  pair<NodeSet::iterator, bool> testIns =
906  cvertlist.insert(tmp[i]);
907 
908  if (testIns.second) // have inserted node
909  {
910  (*(testIns.first))->m_id = newidx;
911 
913  v.id = newidx;
914  v.x = tmp[i]->m_x;
915  v.y = tmp[i]->m_y;
916  v.z = tmp[i]->m_z;
917  curvedpts.pts.push_back(v);
918  newidx++;
919  }
920  curvedpts.index.push_back((*(testIns.first))->m_id);
921  }
922  ptoffset += cinfo.npoints;
923  }
924  }
925  }
926 
927  // add xml information
928  if (edgeinfo.size())
929  {
930  curved->SetAttribute(
932  curved->SetAttribute("BITSIZE",
934 
935  TiXmlElement *x = new TiXmlElement("E");
936  std::string dataStr;
938  dataStr);
939  x->LinkEndChild(new TiXmlText(dataStr));
940  curved->LinkEndChild(x);
941  }
942 
943  if (faceinfo.size())
944  {
945  curved->SetAttribute(
947  curved->SetAttribute("BITSIZE",
949 
950  TiXmlElement *x = new TiXmlElement("F");
951  std::string dataStr;
953  dataStr);
954  x->LinkEndChild(new TiXmlText(dataStr));
955  curved->LinkEndChild(x);
956  }
957 
958  if (edgeinfo.size() || faceinfo.size())
959  {
960  TiXmlElement *x = new TiXmlElement("DATAPOINTS");
961  x->SetAttribute("ID", curvedpts.id);
962 
963  TiXmlElement *subx = new TiXmlElement("INDEX");
964  std::string dataStr;
966  dataStr);
967  subx->LinkEndChild(new TiXmlText(dataStr));
968  x->LinkEndChild(subx);
969 
970  subx = new TiXmlElement("POINTS");
972  dataStr);
973  subx->LinkEndChild(new TiXmlText(dataStr));
974  x->LinkEndChild(subx);
975 
976  curved->LinkEndChild(x);
977  }
978  }
979  pRoot->LinkEndChild(curved);
980 }
981 
982 void OutputNekpp::WriteXmlComposites(TiXmlElement *pRoot)
983 {
984  TiXmlElement *verTag = new TiXmlElement("COMPOSITE");
987  int j = 0;
988 
989  for (it = m_mesh->m_composite.begin(); it != m_mesh->m_composite.end();
990  ++it, ++j)
991  {
992  if (it->second->m_items.size() > 0)
993  {
994  TiXmlElement *comp_tag = new TiXmlElement("C"); // Composite
995  bool doSort = true;
996 
997  // Ensure that this composite is not used for periodic BCs!
998  for (it2 = m_mesh->m_condition.begin();
999  it2 != m_mesh->m_condition.end();
1000  ++it2)
1001  {
1002  ConditionSharedPtr c = it2->second;
1003 
1004  // Ignore non-periodic boundary conditions.
1005  if (find(c->type.begin(), c->type.end(), ePeriodic) ==
1006  c->type.end())
1007  {
1008  continue;
1009  }
1010 
1011  for (int i = 0; i < c->m_composite.size(); ++i)
1012  {
1013  if (c->m_composite[i] == j)
1014  {
1015  doSort = false;
1016  }
1017  }
1018  }
1019 
1020  doSort = doSort && it->second->m_reorder;
1021  comp_tag->SetAttribute("ID", it->second->m_id);
1022  if (it->second->m_label.size())
1023  {
1024  comp_tag->SetAttribute("LABEL", it->second->m_label);
1025  }
1026  comp_tag->LinkEndChild(
1027  new TiXmlText(it->second->GetXmlString(doSort)));
1028  verTag->LinkEndChild(comp_tag);
1029  }
1030  else
1031  {
1032  cout << "Composite " << it->second->m_id << " "
1033  << "contains nothing." << endl;
1034  }
1035  }
1036 
1037  pRoot->LinkEndChild(verTag);
1038 }
1039 
1040 void OutputNekpp::WriteXmlDomain(TiXmlElement *pRoot)
1041 {
1042  // Write the <DOMAIN> subsection.
1043  TiXmlElement *domain = new TiXmlElement("DOMAIN");
1044  std::string list;
1046 
1047  for (it = m_mesh->m_composite.begin(); it != m_mesh->m_composite.end();
1048  ++it)
1049  {
1050  if (it->second->m_items[0]->GetDim() == m_mesh->m_expDim)
1051  {
1052  if (list.length() > 0)
1053  {
1054  list += ",";
1055  }
1056  list += boost::lexical_cast<std::string>(it->second->m_id);
1057  }
1058  }
1059  domain->LinkEndChild(new TiXmlText(" C[" + list + "] "));
1060  pRoot->LinkEndChild(domain);
1061 }
1062 
1063 void OutputNekpp::WriteXmlExpansions(TiXmlElement *pRoot)
1064 {
1065  // Write a default <EXPANSIONS> section.
1066  TiXmlElement *expansions = new TiXmlElement("EXPANSIONS");
1068 
1069  for (it = m_mesh->m_composite.begin(); it != m_mesh->m_composite.end();
1070  ++it)
1071  {
1072  if (it->second->m_items[0]->GetDim() == m_mesh->m_expDim)
1073  {
1074  TiXmlElement *exp = new TiXmlElement("E");
1075  exp->SetAttribute(
1076  "COMPOSITE",
1077  "C[" + boost::lexical_cast<std::string>(it->second->m_id) +
1078  "]");
1079  exp->SetAttribute("NUMMODES", 4);
1080  exp->SetAttribute("TYPE", "MODIFIED");
1081 
1082  if (m_mesh->m_fields.size() == 0)
1083  {
1084  exp->SetAttribute("FIELDS", "u");
1085  }
1086  else
1087  {
1088  string fstr;
1089  for (int i = 0; i < m_mesh->m_fields.size(); ++i)
1090  {
1091  fstr += m_mesh->m_fields[i] + ",";
1092  }
1093  fstr = fstr.substr(0, fstr.length() - 1);
1094  exp->SetAttribute("FIELDS", fstr);
1095  }
1096 
1097  expansions->LinkEndChild(exp);
1098  }
1099  }
1100  pRoot->LinkEndChild(expansions);
1101 }
1102 
1103 void OutputNekpp::WriteXmlConditions(TiXmlElement *pRoot)
1104 {
1105  TiXmlElement *conditions = new TiXmlElement("CONDITIONS");
1106  TiXmlElement *boundaryregions = new TiXmlElement("BOUNDARYREGIONS");
1107  TiXmlElement *boundaryconditions = new TiXmlElement("BOUNDARYCONDITIONS");
1108  TiXmlElement *variables = new TiXmlElement("VARIABLES");
1110 
1111  for (it = m_mesh->m_condition.begin(); it != m_mesh->m_condition.end();
1112  ++it)
1113  {
1114  ConditionSharedPtr c = it->second;
1115  string tmp;
1116 
1117  // First set up boundary regions.
1118  TiXmlElement *b = new TiXmlElement("B");
1119  b->SetAttribute("ID", boost::lexical_cast<string>(it->first));
1120 
1121  for (int i = 0; i < c->m_composite.size(); ++i)
1122  {
1123  tmp += boost::lexical_cast<string>(c->m_composite[i]) + ",";
1124  }
1125 
1126  tmp = tmp.substr(0, tmp.length() - 1);
1127 
1128  TiXmlText *t0 = new TiXmlText("C[" + tmp + "]");
1129  b->LinkEndChild(t0);
1130  boundaryregions->LinkEndChild(b);
1131 
1132  TiXmlElement *region = new TiXmlElement("REGION");
1133  region->SetAttribute("REF", boost::lexical_cast<string>(it->first));
1134 
1135  for (int i = 0; i < c->type.size(); ++i)
1136  {
1137  string tagId;
1138 
1139  switch (c->type[i])
1140  {
1141  case eDirichlet:
1142  tagId = "D";
1143  break;
1144  case eNeumann:
1145  tagId = "N";
1146  break;
1147  case ePeriodic:
1148  tagId = "P";
1149  break;
1150  case eHOPCondition:
1151  tagId = "N";
1152  break;
1153  default:
1154  break;
1155  }
1156 
1157  TiXmlElement *tag = new TiXmlElement(tagId);
1158  tag->SetAttribute("VAR", c->field[i]);
1159  tag->SetAttribute("VALUE", c->value[i]);
1160 
1161  if (c->type[i] == eHOPCondition)
1162  {
1163  tag->SetAttribute("USERDEFINEDTYPE", "H");
1164  }
1165 
1166  region->LinkEndChild(tag);
1167  }
1168 
1169  boundaryconditions->LinkEndChild(region);
1170  }
1171 
1172  for (int i = 0; i < m_mesh->m_fields.size(); ++i)
1173  {
1174  TiXmlElement *v = new TiXmlElement("V");
1175  v->SetAttribute("ID", boost::lexical_cast<std::string>(i));
1176  TiXmlText *t0 = new TiXmlText(m_mesh->m_fields[i]);
1177  v->LinkEndChild(t0);
1178  variables->LinkEndChild(v);
1179  }
1180 
1181  if (m_mesh->m_fields.size() > 0)
1182  {
1183  conditions->LinkEndChild(variables);
1184  }
1185 
1186  if (m_mesh->m_condition.size() > 0)
1187  {
1188  conditions->LinkEndChild(boundaryregions);
1189  conditions->LinkEndChild(boundaryconditions);
1190  }
1191 
1192  pRoot->LinkEndChild(conditions);
1193 }
1194 }
1195 }
void WriteXmlCurves(TiXmlElement *pRoot)
Writes the section of the XML file if needed.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void Process()
Write mesh to output file.
Definition: OutputNekpp.cpp:94
pair< ModuleType, string > ModuleKey
static boost::shared_ptr< MeshGraph > Read(const LibUtilities::SessionReaderSharedPtr &pSession, DomainRangeShPtr &rng=NullDomainRangeShPtr)
Definition: MeshGraph.cpp:121
NekInt64 ptoffset
the id of point data map (currently always 0 since we are using just one set).
void WriteXmlFaces(TiXmlElement *pRoot)
Writes the section of the XML file if needed.
Abstract base class for output modules.
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:69
void WriteXmlNodes(TiXmlElement *pRoot)
Writes the section of the XML file.
std::vector< NekInt64 > index
id of this Point set
map< string, ConfigOption > m_config
List of configuration values.
STL namespace.
MeshSharedPtr m_mesh
Mesh object.
void TestElmts(SpatialDomains::MeshGraphSharedPtr &graph)
Definition: OutputNekpp.cpp:77
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
NekInt64 entityid
Id of this curved information.
boost::unordered_set< NodeSharedPtr, NodeHash > NodeSet
Definition: Node.h:357
static SessionReaderSharedPtr CreateInstance(int argc, char *argv[])
Creates an instance of the SessionReader class.
string GetXmlString(char tag, vector< unsigned int > &ids)
NekInt64 npoints
The entity id corresponding to the global edge/curve.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
boost::shared_ptr< Condition > ConditionSharedPtr
Definition: Mesh.h:78
int ZlibEncodeToBase64Str(std::vector< T > &in, std::string &out64)
Definition: CompressData.h:151
Represents a command-line configuration option.
NekInt64 ptype
point offset of data entry for this curve
void WriteXmlExpansions(TiXmlElement *pRoot)
Writes the section of the XML file.
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:196
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:137
void WriteXmlDomain(TiXmlElement *pRoot)
Writes the section of the XML file.
void WriteXmlComposites(TiXmlElement *pRoot)
Writes the section of the XML file.
void WriteXmlElements(TiXmlElement *pRoot)
Writes the section of the XML file.
NekInt64 ptid
The number of points in this curved entity.
std::vector< MeshVertex > pts
mapping to access pts value.
InputIterator find(InputIterator first, InputIterator last, InputIterator startingpoint, const EqualityComparable &value)
Definition: StdRegions.hpp:315
void WriteXmlConditions(TiXmlElement *pRoot)
Writes the section of the XML file.
boost::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:63
boost::shared_ptr< Face > FaceSharedPtr
Shared pointer to a face.
Definition: Face.h:378
void WriteXmlEdges(TiXmlElement *pRoot)
Writes the section of the XML file.
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215