Nektar++
MeshGraphXmlCompressed.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: BoundaryConditions.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description:
32 //
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include "MeshGraphXmlCompressed.h"
37 
42 
45 
46 // These are required for the Write(...) and Import(...) functions.
47 #include <boost/archive/iterators/base64_from_binary.hpp>
48 #include <boost/archive/iterators/binary_from_base64.hpp>
49 #include <boost/archive/iterators/transform_width.hpp>
50 #include <boost/algorithm/string/predicate.hpp>
51 #include <boost/format.hpp>
52 #include <boost/iostreams/copy.hpp>
53 #include <boost/iostreams/filter/zlib.hpp>
54 #include <boost/iostreams/filtering_stream.hpp>
55 #include <boost/make_shared.hpp>
56 
57 #include <tinyxml.h>
58 using namespace std;
59 
60 namespace Nektar
61 {
62 namespace SpatialDomains
63 {
64 
65 std::string MeshGraphXmlCompressed::className =
67  "XmlCompressed", MeshGraphXmlCompressed::create,
68  "IO with Xml geometry");
69 
70 void MeshGraphXmlCompressed::ReadVertices()
71 {
72  // Now read the vertices
73  TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
74  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
75 
76  NekDouble xscale, yscale, zscale;
77 
78  // check to see if any scaling parameters are in
79  // attributes and determine these values
80  LibUtilities::Interpreter expEvaluator;
81  const char *xscal = element->Attribute("XSCALE");
82  if (!xscal)
83  {
84  xscale = 1.0;
85  }
86  else
87  {
88  std::string xscalstr = xscal;
89  int expr_id = expEvaluator.DefineFunction("", xscalstr);
90  xscale = expEvaluator.Evaluate(expr_id);
91  }
92 
93  const char *yscal = element->Attribute("YSCALE");
94  if (!yscal)
95  {
96  yscale = 1.0;
97  }
98  else
99  {
100  std::string yscalstr = yscal;
101  int expr_id = expEvaluator.DefineFunction("", yscalstr);
102  yscale = expEvaluator.Evaluate(expr_id);
103  }
104 
105  const char *zscal = element->Attribute("ZSCALE");
106  if (!zscal)
107  {
108  zscale = 1.0;
109  }
110  else
111  {
112  std::string zscalstr = zscal;
113  int expr_id = expEvaluator.DefineFunction("", zscalstr);
114  zscale = expEvaluator.Evaluate(expr_id);
115  }
116 
117  NekDouble xmove, ymove, zmove;
118 
119  // check to see if any moving parameters are in
120  // attributes and determine these values
121 
122  const char *xmov = element->Attribute("XMOVE");
123  if (!xmov)
124  {
125  xmove = 0.0;
126  }
127  else
128  {
129  std::string xmovstr = xmov;
130  int expr_id = expEvaluator.DefineFunction("", xmovstr);
131  xmove = expEvaluator.Evaluate(expr_id);
132  }
133 
134  const char *ymov = element->Attribute("YMOVE");
135  if (!ymov)
136  {
137  ymove = 0.0;
138  }
139  else
140  {
141  std::string ymovstr = ymov;
142  int expr_id = expEvaluator.DefineFunction("", ymovstr);
143  ymove = expEvaluator.Evaluate(expr_id);
144  }
145 
146  const char *zmov = element->Attribute("ZMOVE");
147  if (!zmov)
148  {
149  zmove = 0.0;
150  }
151  else
152  {
153  std::string zmovstr = zmov;
154  int expr_id = expEvaluator.DefineFunction("", zmovstr);
155  zmove = expEvaluator.Evaluate(expr_id);
156  }
157 
158  string IsCompressed;
159  element->QueryStringAttribute("COMPRESSED", &IsCompressed);
160 
161  if (boost::iequals(IsCompressed,
163  {
164  // Extract the vertex body
165  TiXmlNode *vertexChild = element->FirstChild();
166  ASSERTL0(vertexChild, "Unable to extract the data from the compressed "
167  "vertex tag.");
168 
169  std::string vertexStr;
170  if (vertexChild->Type() == TiXmlNode::TINYXML_TEXT)
171  {
172  vertexStr += vertexChild->ToText()->ValueStr();
173  }
174 
175  std::vector<SpatialDomains::MeshVertex> vertData;
177  vertData);
178 
179  int indx;
180  NekDouble xval, yval, zval;
181  for (int i = 0; i < vertData.size(); ++i)
182  {
183  indx = vertData[i].id;
184  xval = vertData[i].x;
185  yval = vertData[i].y;
186  zval = vertData[i].z;
187 
188  xval = xval * xscale + xmove;
189  yval = yval * yscale + ymove;
190  zval = zval * zscale + zmove;
191 
193  m_spaceDimension, indx, xval, yval, zval));
194 
195  vert->SetGlobalID(indx);
196  m_vertSet[indx] = vert;
197  }
198  }
199  else
200  {
201  ASSERTL0(false, "Compressed formats do not match. Expected :" +
203  " but got " + std::string(IsCompressed));
204  }
205 }
206 
207 void MeshGraphXmlCompressed::ReadCurves()
208 {
209  // check to see if any scaling parameters are in
210  // attributes and determine these values
211  TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
212  ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
213 
214  NekDouble xscale, yscale, zscale;
215 
216  LibUtilities::Interpreter expEvaluator;
217  const char *xscal = element->Attribute("XSCALE");
218  if (!xscal)
219  {
220  xscale = 1.0;
221  }
222  else
223  {
224  std::string xscalstr = xscal;
225  int expr_id = expEvaluator.DefineFunction("", xscalstr);
226  xscale = expEvaluator.Evaluate(expr_id);
227  }
228 
229  const char *yscal = element->Attribute("YSCALE");
230  if (!yscal)
231  {
232  yscale = 1.0;
233  }
234  else
235  {
236  std::string yscalstr = yscal;
237  int expr_id = expEvaluator.DefineFunction("", yscalstr);
238  yscale = expEvaluator.Evaluate(expr_id);
239  }
240 
241  const char *zscal = element->Attribute("ZSCALE");
242  if (!zscal)
243  {
244  zscale = 1.0;
245  }
246  else
247  {
248  std::string zscalstr = zscal;
249  int expr_id = expEvaluator.DefineFunction("", zscalstr);
250  zscale = expEvaluator.Evaluate(expr_id);
251  }
252 
253  NekDouble xmove, ymove, zmove;
254 
255  // check to see if any moving parameters are in
256  // attributes and determine these values
257 
258  const char *xmov = element->Attribute("XMOVE");
259  if (!xmov)
260  {
261  xmove = 0.0;
262  }
263  else
264  {
265  std::string xmovstr = xmov;
266  int expr_id = expEvaluator.DefineFunction("", xmovstr);
267  xmove = expEvaluator.Evaluate(expr_id);
268  }
269 
270  const char *ymov = element->Attribute("YMOVE");
271  if (!ymov)
272  {
273  ymove = 0.0;
274  }
275  else
276  {
277  std::string ymovstr = ymov;
278  int expr_id = expEvaluator.DefineFunction("", ymovstr);
279  ymove = expEvaluator.Evaluate(expr_id);
280  }
281 
282  const char *zmov = element->Attribute("ZMOVE");
283  if (!zmov)
284  {
285  zmove = 0.0;
286  }
287  else
288  {
289  std::string zmovstr = zmov;
290  int expr_id = expEvaluator.DefineFunction("", zmovstr);
291  zmove = expEvaluator.Evaluate(expr_id);
292  }
293 
294  /// Look for elements in CURVE block.
295  TiXmlElement *field = m_xmlGeom->FirstChildElement("CURVED");
296 
297  if (!field) // return if no curved entities
298  {
299  return;
300  }
301 
302  string IsCompressed;
303  field->QueryStringAttribute("COMPRESSED", &IsCompressed);
304 
305  if (IsCompressed.size() == 0)
306  {
307  // this could be that the curved tag is empty
308  // in this case we dont want to read it
309  return;
310  }
311 
312  ASSERTL0(boost::iequals(IsCompressed,
314  "Compressed formats do not match. Expected :" +
316  boost::lexical_cast<std::string>(IsCompressed));
317 
318  std::vector<SpatialDomains::MeshCurvedInfo> edginfo;
319  std::vector<SpatialDomains::MeshCurvedInfo> facinfo;
321 
322  // read edge, face info and curved poitns.
323  TiXmlElement *x = field->FirstChildElement();
324  while (x)
325  {
326  const char *entitytype = x->Value();
327  // read in edge or face info
328  if (boost::iequals(entitytype, "E"))
329  {
330  // read in data
331  std::string elmtStr;
332  TiXmlNode *child = x->FirstChild();
333 
334  if (child->Type() == TiXmlNode::TINYXML_TEXT)
335  {
336  elmtStr += child->ToText()->ValueStr();
337  }
338 
340  edginfo);
341  }
342  else if (boost::iequals(entitytype, "F"))
343  {
344  // read in data
345  std::string elmtStr;
346  TiXmlNode *child = x->FirstChild();
347 
348  if (child->Type() == TiXmlNode::TINYXML_TEXT)
349  {
350  elmtStr += child->ToText()->ValueStr();
351  }
352 
354  facinfo);
355  }
356  else if (boost::iequals(entitytype, "DATAPOINTS"))
357  {
358  NekInt id;
359  ASSERTL0(x->Attribute("ID", &id),
360  "Failed to get ID from PTS section");
361  cpts.id = id;
362 
363  // read in data
364  std::string elmtStr;
365 
366  TiXmlElement *DataIdx = x->FirstChildElement("INDEX");
367  ASSERTL0(DataIdx, "Cannot read data index tag in compressed "
368  "curved section");
369 
370  TiXmlNode *child = DataIdx->FirstChild();
371  if (child->Type() == TiXmlNode::TINYXML_TEXT)
372  {
373  elmtStr = child->ToText()->ValueStr();
374  }
375 
377  cpts.index);
378 
379  TiXmlElement *DataPts = x->FirstChildElement("POINTS");
380  ASSERTL0(DataPts, "Cannot read data pts tag in compressed "
381  "curved section");
382 
383  child = DataPts->FirstChild();
384  if (child->Type() == TiXmlNode::TINYXML_TEXT)
385  {
386  elmtStr = child->ToText()->ValueStr();
387  }
388 
390  cpts.pts);
391  }
392  else
393  {
394  ASSERTL0(false, "Unknown tag in curved section");
395  }
396  x = x->NextSiblingElement();
397  }
398 
399  // rescale (x,y,z) points;
400  for (int i = 0; i < cpts.pts.size(); ++i)
401  {
402  cpts.pts[i].x = xscale * cpts.pts[i].x + xmove;
403  cpts.pts[i].y = yscale * cpts.pts[i].y + ymove;
404  cpts.pts[i].z = zscale * cpts.pts[i].z + zmove;
405  }
406 
407  for (int i = 0; i < edginfo.size(); ++i)
408  {
409  int edgeid = edginfo[i].entityid;
411 
413  edgeid, ptype = (LibUtilities::PointsType)edginfo[i].ptype));
414 
415  // load points
416  int offset = edginfo[i].ptoffset;
417  for (int j = 0; j < edginfo[i].npoints; ++j)
418  {
419  int idx = cpts.index[offset + j];
420 
422  m_spaceDimension, edginfo[i].id, cpts.pts[idx].x,
423  cpts.pts[idx].y, cpts.pts[idx].z));
424  curve->m_points.push_back(vert);
425  }
426 
427  m_curvedEdges[edgeid] = curve;
428  }
429 
430  for (int i = 0; i < facinfo.size(); ++i)
431  {
432  int faceid = facinfo[i].entityid;
434 
436  faceid, ptype = (LibUtilities::PointsType)facinfo[i].ptype));
437 
438  int offset = facinfo[i].ptoffset;
439  for (int j = 0; j < facinfo[i].npoints; ++j)
440  {
441  int idx = cpts.index[offset + j];
442 
444  m_spaceDimension, facinfo[i].id, cpts.pts[idx].x,
445  cpts.pts[idx].y, cpts.pts[idx].z));
446  curve->m_points.push_back(vert);
447  }
448 
449  m_curvedFaces[faceid] = curve;
450  }
451 }
452 
453 void MeshGraphXmlCompressed::ReadEdges()
454 {
455  CurveMap::iterator it;
456 
457  /// Look for elements in ELEMENT block.
458  TiXmlElement *field = m_xmlGeom->FirstChildElement("EDGE");
459 
460  ASSERTL0(field, "Unable to find EDGE tag in file.");
461 
462  string IsCompressed;
463  field->QueryStringAttribute("COMPRESSED", &IsCompressed);
464 
465  ASSERTL0(boost::iequals(IsCompressed,
467  "Compressed formats do not match. Expected :" +
469  std::string(IsCompressed));
470  // Extract the edge body
471  TiXmlNode *edgeChild = field->FirstChild();
472  ASSERTL0(edgeChild, "Unable to extract the data from "
473  "the compressed edge tag.");
474 
475  std::string edgeStr;
476  if (edgeChild->Type() == TiXmlNode::TINYXML_TEXT)
477  {
478  edgeStr += edgeChild->ToText()->ValueStr();
479  }
480 
481  std::vector<SpatialDomains::MeshEdge> edgeData;
483 
484  int indx;
485  for (int i = 0; i < edgeData.size(); ++i)
486  {
487  indx = edgeData[i].id;
488  PointGeomSharedPtr vertices[2] = {GetVertex(edgeData[i].v0),
489  GetVertex(edgeData[i].v1)};
490  SegGeomSharedPtr edge;
491 
492  it = m_curvedEdges.find(indx);
493  if (it == m_curvedEdges.end())
494  {
496  indx, m_spaceDimension, vertices);
497  }
498  else
499  {
501  indx, m_spaceDimension, vertices, it->second);
502  }
503  m_segGeoms[indx] = edge;
504  }
505 }
506 
507 void MeshGraphXmlCompressed::ReadFaces()
508 {
509  /// Look for elements in FACE block.
510  TiXmlElement *field = m_xmlGeom->FirstChildElement("FACE");
511 
512  ASSERTL0(field, "Unable to find FACE tag in file.");
513 
514  /// All faces are of the form: "<? ID="#"> ... </?>", with
515  /// ? being an element type (either Q or T).
516  /// They might be in compressed format and so then need upacking.
517 
518  TiXmlElement *element = field->FirstChildElement();
519  CurveMap::iterator it;
520 
521  while (element)
522  {
523  std::string elementType(element->ValueStr());
524 
525  ASSERTL0(elementType == "Q" || elementType == "T",
526  (std::string("Unknown 3D face type: ") + elementType).c_str());
527 
528  string IsCompressed;
529  element->QueryStringAttribute("COMPRESSED", &IsCompressed);
530 
531  ASSERTL0(
532  boost::iequals(IsCompressed,
534  "Compressed formats do not match. Expected :" +
536  std::string(IsCompressed));
537 
538  // Extract the face body
539  TiXmlNode *faceChild = element->FirstChild();
540  ASSERTL0(faceChild, "Unable to extract the data from "
541  "the compressed face tag.");
542 
543  std::string faceStr;
544  if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
545  {
546  faceStr += faceChild->ToText()->ValueStr();
547  }
548 
549  int indx;
550  if (elementType == "T")
551  {
552  std::vector<SpatialDomains::MeshTri> faceData;
554  faceData);
555 
556  for (int i = 0; i < faceData.size(); ++i)
557  {
558  indx = faceData[i].id;
559 
560  /// See if this face has curves.
561  it = m_curvedFaces.find(indx);
562 
563  /// Create a TriGeom to hold the new definition.
564  SegGeomSharedPtr edges[TriGeom::kNedges] = {
565  GetSegGeom(faceData[i].e[0]), GetSegGeom(faceData[i].e[1]),
566  GetSegGeom(faceData[i].e[2])};
567 
568  TriGeomSharedPtr trigeom;
569  if (it == m_curvedFaces.end())
570  {
571  trigeom =
573  }
574  else
575  {
577  indx, edges, it->second);
578  }
579  trigeom->SetGlobalID(indx);
580  m_triGeoms[indx] = trigeom;
581  }
582  }
583  else if (elementType == "Q")
584  {
585  std::vector<SpatialDomains::MeshQuad> faceData;
587  faceData);
588 
589  for (int i = 0; i < faceData.size(); ++i)
590  {
591  indx = faceData[i].id;
592 
593  /// See if this face has curves.
594  it = m_curvedFaces.find(indx);
595 
596  /// Create a QuadGeom to hold the new definition.
597  SegGeomSharedPtr edges[QuadGeom::kNedges] = {
598  GetSegGeom(faceData[i].e[0]), GetSegGeom(faceData[i].e[1]),
599  GetSegGeom(faceData[i].e[2]), GetSegGeom(faceData[i].e[3])};
600 
601  QuadGeomSharedPtr quadgeom;
602  if (it == m_curvedFaces.end())
603  {
604  quadgeom =
606  }
607  else
608  {
610  indx, edges, it->second);
611  }
612  quadgeom->SetGlobalID(indx);
613  m_quadGeoms[indx] = quadgeom;
614  }
615  }
616  /// Keep looking
617  element = element->NextSiblingElement();
618  }
619 }
620 
621 void MeshGraphXmlCompressed::ReadElements1D()
622 {
623  TiXmlElement *field = NULL;
624 
625  /// Look for elements in ELEMENT block.
626  field = m_xmlGeom->FirstChildElement("ELEMENT");
627 
628  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
629 
630  /// All elements are of the form: "<S ID = n> ... </S>", with
631  /// ? being the element type.
632 
633  TiXmlElement *segment = field->FirstChildElement("S");
634  CurveMap::iterator it;
635 
636  while (segment)
637  {
638  string IsCompressed;
639  segment->QueryStringAttribute("COMPRESSED", &IsCompressed);
640  ASSERTL0(
641  boost::iequals(IsCompressed,
643  "Compressed formats do not match. Expected :" +
645  std::string(IsCompressed));
646 
647  // Extract the face body
648  TiXmlNode *child = segment->FirstChild();
649  ASSERTL0(child, "Unable to extract the data from "
650  "the compressed face tag.");
651 
652  std::string str;
653  if (child->Type() == TiXmlNode::TINYXML_TEXT)
654  {
655  str += child->ToText()->ValueStr();
656  }
657 
658  int indx;
659 
660  std::vector<SpatialDomains::MeshEdge> data;
662 
663  for (int i = 0; i < data.size(); ++i)
664  {
665  indx = data[i].id;
666 
667  /// See if this face has curves.
668  it = m_curvedEdges.find(indx);
669 
670  PointGeomSharedPtr vertices[2] = {GetVertex(data[i].v0),
671  GetVertex(data[i].v1)};
672  SegGeomSharedPtr seg;
673 
674  if (it == m_curvedEdges.end())
675  {
677  indx, m_spaceDimension, vertices);
678  seg->SetGlobalID(indx); // Set global mesh id
679  }
680  else
681  {
683  indx, m_spaceDimension, vertices, it->second);
684  seg->SetGlobalID(indx); // Set global mesh id
685  }
686  seg->SetGlobalID(indx);
687  m_segGeoms[indx] = seg;
688  }
689  /// Keep looking for additional segments
690  segment = segment->NextSiblingElement("S");
691  }
692 }
693 
694 void MeshGraphXmlCompressed::ReadElements2D()
695 {
696  /// Look for elements in ELEMENT block.
697  TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
698 
699  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
700 
701  // Set up curve map for curved elements on an embedded manifold.
702  CurveMap::iterator it;
703 
704  /// All elements are of the form: "<? ID="#"> ... </?>", with
705  /// ? being the element type.
706 
707  TiXmlElement *element = field->FirstChildElement();
708 
709  while (element)
710  {
711  std::string elementType(element->ValueStr());
712 
713  ASSERTL0(
714  elementType == "Q" || elementType == "T",
715  (std::string("Unknown 2D element type: ") + elementType).c_str());
716 
717  string IsCompressed;
718  element->QueryStringAttribute("COMPRESSED", &IsCompressed);
719 
720  ASSERTL0(
721  boost::iequals(IsCompressed,
723  "Compressed formats do not match. Expected :" +
725  std::string(IsCompressed));
726 
727  // Extract the face body
728  TiXmlNode *faceChild = element->FirstChild();
729  ASSERTL0(faceChild, "Unable to extract the data from "
730  "the compressed face tag.");
731 
732  std::string faceStr;
733  if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
734  {
735  faceStr += faceChild->ToText()->ValueStr();
736  }
737 
738  int indx;
739  if (elementType == "T")
740  {
741  std::vector<SpatialDomains::MeshTri> faceData;
743  faceData);
744 
745  for (int i = 0; i < faceData.size(); ++i)
746  {
747  indx = faceData[i].id;
748 
749  /// See if this face has curves.
750  it = m_curvedFaces.find(indx);
751 
752  /// Create a TriGeom to hold the new definition.
753  SegGeomSharedPtr edges[TriGeom::kNedges] = {
754  GetSegGeom(faceData[i].e[0]), GetSegGeom(faceData[i].e[1]),
755  GetSegGeom(faceData[i].e[2])};
756 
757  TriGeomSharedPtr trigeom;
758  if (it == m_curvedFaces.end())
759  {
760  trigeom =
762  }
763  else
764  {
766  indx, edges, it->second);
767  }
768  trigeom->SetGlobalID(indx);
769  m_triGeoms[indx] = trigeom;
770  }
771  }
772  else if (elementType == "Q")
773  {
774  std::vector<SpatialDomains::MeshQuad> faceData;
776  faceData);
777 
778  for (int i = 0; i < faceData.size(); ++i)
779  {
780  indx = faceData[i].id;
781 
782  /// See if this face has curves.
783  it = m_curvedFaces.find(indx);
784 
785  /// Create a QuadGeom to hold the new definition.
786  SegGeomSharedPtr edges[QuadGeom::kNedges] = {
787  GetSegGeom(faceData[i].e[0]), GetSegGeom(faceData[i].e[1]),
788  GetSegGeom(faceData[i].e[2]), GetSegGeom(faceData[i].e[3])};
789 
790  QuadGeomSharedPtr quadgeom;
791  if (it == m_curvedFaces.end())
792  {
793  quadgeom =
795  }
796  else
797  {
799  indx, edges, it->second);
800  }
801  quadgeom->SetGlobalID(indx);
802  m_quadGeoms[indx] = quadgeom;
803  }
804  }
805  /// Keep looking
806  element = element->NextSiblingElement();
807  }
808 }
809 
810 void MeshGraphXmlCompressed::ReadElements3D()
811 {
812  /// Look for elements in ELEMENT block.
813  TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
814 
815  ASSERTL0(field, "Unable to find ELEMENT tag in file.");
816 
817  /// All elements are of the form: "<? ID="#"> ... </?>", with
818  /// ? being the element type.
819 
820  TiXmlElement *element = field->FirstChildElement();
821 
822  while (element)
823  {
824  std::string elementType(element->ValueStr());
825 
826  // A - tet, P - pyramid, R - prism, H - hex
827  ASSERTL0(
828  elementType == "A" || elementType == "P" || elementType == "R" ||
829  elementType == "H",
830  (std::string("Unknown 3D element type: ") + elementType).c_str());
831 
832  string IsCompressed;
833  element->QueryStringAttribute("COMPRESSED", &IsCompressed);
834 
835  ASSERTL0(
836  boost::iequals(IsCompressed,
838  "Compressed formats do not match. Expected :" +
840  std::string(IsCompressed));
841 
842  // Extract the face body
843  TiXmlNode *child = element->FirstChild();
844  ASSERTL0(child, "Unable to extract the data from "
845  "the compressed face tag.");
846 
847  std::string str;
848  if (child->Type() == TiXmlNode::TINYXML_TEXT)
849  {
850  str += child->ToText()->ValueStr();
851  }
852 
853  int indx;
854  if (elementType == "A")
855  {
856  std::vector<SpatialDomains::MeshTet> data;
858  TriGeomSharedPtr tfaces[4];
859  for (int i = 0; i < data.size(); ++i)
860  {
861  indx = data[i].id;
862  for (int j = 0; j < 4; ++j)
863  {
864  Geometry2DSharedPtr face = GetGeometry2D(data[i].f[j]);
865  tfaces[j] = static_pointer_cast<TriGeom>(face);
866  }
867 
868  TetGeomSharedPtr tetgeom(
870  m_tetGeoms[indx] = tetgeom;
871  PopulateFaceToElMap(tetgeom, 4);
872  }
873  }
874  else if (elementType == "P")
875  {
876  std::vector<SpatialDomains::MeshPyr> data;
878  Geometry2DSharedPtr faces[5];
879  for (int i = 0; i < data.size(); ++i)
880  {
881  indx = data[i].id;
882  int Ntfaces = 0;
883  int Nqfaces = 0;
884  for (int j = 0; j < 5; ++j)
885  {
886  Geometry2DSharedPtr face = GetGeometry2D(data[i].f[j]);
887 
888  if (face == Geometry2DSharedPtr() ||
889  (face->GetShapeType() != LibUtilities::eTriangle &&
890  face->GetShapeType() != LibUtilities::eQuadrilateral))
891  {
892  std::stringstream errorstring;
893  errorstring << "Element " << indx
894  << " has invalid face: " << j;
895  ASSERTL0(false, errorstring.str().c_str());
896  }
897  else if (face->GetShapeType() == LibUtilities::eTriangle)
898  {
899  faces[j] = static_pointer_cast<TriGeom>(face);
900  Ntfaces++;
901  }
902  else if (face->GetShapeType() ==
904  {
905  faces[j] = static_pointer_cast<QuadGeom>(face);
906  Nqfaces++;
907  }
908  }
909  ASSERTL0((Ntfaces == 4) && (Nqfaces == 1),
910  "Did not identify the correct number of "
911  "triangular and quadrilateral faces for a "
912  "pyramid");
913 
914  PyrGeomSharedPtr pyrgeom(
916 
917  m_pyrGeoms[indx] = pyrgeom;
918  PopulateFaceToElMap(pyrgeom, 5);
919  }
920  }
921  else if (elementType == "R")
922  {
923  std::vector<SpatialDomains::MeshPrism> data;
925  Geometry2DSharedPtr faces[5];
926  for (int i = 0; i < data.size(); ++i)
927  {
928  indx = data[i].id;
929  int Ntfaces = 0;
930  int Nqfaces = 0;
931  for (int j = 0; j < 5; ++j)
932  {
933  Geometry2DSharedPtr face = GetGeometry2D(data[i].f[j]);
934  if (face == Geometry2DSharedPtr() ||
935  (face->GetShapeType() != LibUtilities::eTriangle &&
936  face->GetShapeType() != LibUtilities::eQuadrilateral))
937  {
938  std::stringstream errorstring;
939  errorstring << "Element " << indx
940  << " has invalid face: " << j;
941  ASSERTL0(false, errorstring.str().c_str());
942  }
943  else if (face->GetShapeType() == LibUtilities::eTriangle)
944  {
945  faces[j] = static_pointer_cast<TriGeom>(face);
946  Ntfaces++;
947  }
948  else if (face->GetShapeType() ==
950  {
951  faces[j] = static_pointer_cast<QuadGeom>(face);
952  Nqfaces++;
953  }
954  }
955  ASSERTL0((Ntfaces == 2) && (Nqfaces == 3),
956  "Did not identify the correct number of "
957  "triangular and quadrilateral faces for a "
958  "prism");
959 
960  PrismGeomSharedPtr prismgeom(
962 
963  m_prismGeoms[indx] = prismgeom;
964  PopulateFaceToElMap(prismgeom, 5);
965  }
966  }
967  else if (elementType == "H")
968  {
969  std::vector<SpatialDomains::MeshHex> data;
971 
972  QuadGeomSharedPtr faces[6];
973  for (int i = 0; i < data.size(); ++i)
974  {
975  indx = data[i].id;
976  for (int j = 0; j < 6; ++j)
977  {
978  Geometry2DSharedPtr face = GetGeometry2D(data[i].f[j]);
979  faces[j] = static_pointer_cast<QuadGeom>(face);
980  }
981 
982  HexGeomSharedPtr hexgeom(
984  m_hexGeoms[indx] = hexgeom;
985  PopulateFaceToElMap(hexgeom, 6);
986  }
987  }
988  /// Keep looking
989  element = element->NextSiblingElement();
990  }
991 }
992 
993 void MeshGraphXmlCompressed::WriteVertices(TiXmlElement *geomTag,
994  PointGeomMap &verts)
995 {
996  if (verts.size() == 0)
997  {
998  return;
999  }
1000 
1001  TiXmlElement *vertTag = new TiXmlElement("VERTEX");
1002 
1003  vector<MeshVertex> vertInfo;
1004 
1005  for (auto &i : verts)
1006  {
1007  MeshVertex v;
1008  v.id = i.first;
1009  v.x = i.second->x();
1010  v.y = i.second->y();
1011  v.z = i.second->z();
1012  vertInfo.push_back(v);
1013  }
1014 
1015  vertTag->SetAttribute("COMPRESSED",
1017  vertTag->SetAttribute("BITSIZE",
1019 
1020  string vertStr;
1022 
1023  vertTag->LinkEndChild(new TiXmlText(vertStr));
1024 
1025  geomTag->LinkEndChild(vertTag);
1026 }
1027 
1028 void MeshGraphXmlCompressed::WriteEdges(TiXmlElement *geomTag,
1029  SegGeomMap &edges)
1030 {
1031  if (edges.size() == 0)
1032  {
1033  return;
1034  }
1035 
1036  TiXmlElement *edgeTag =
1037  new TiXmlElement(m_meshDimension == 1 ? "S" : "EDGE");
1038 
1039  vector<MeshEdge> edgeInfo;
1040 
1041  for (auto &i : edges)
1042  {
1043  MeshEdge e;
1044  e.id = i.first;
1045  e.v0 = i.second->GetVid(0);
1046  e.v1 = i.second->GetVid(1);
1047  edgeInfo.push_back(e);
1048  }
1049 
1050  string edgeStr;
1052 
1053  edgeTag->SetAttribute("COMPRESSED",
1055  edgeTag->SetAttribute("BITSIZE",
1057 
1058  edgeTag->LinkEndChild(new TiXmlText(edgeStr));
1059 
1060  if (m_meshDimension == 1)
1061  {
1062  TiXmlElement *tmp = new TiXmlElement("ELEMENT");
1063  tmp->LinkEndChild(edgeTag);
1064  geomTag->LinkEndChild(tmp);
1065  }
1066  else
1067  {
1068  geomTag->LinkEndChild(edgeTag);
1069  }
1070 }
1071 
1072 void MeshGraphXmlCompressed::WriteTris(TiXmlElement *faceTag, TriGeomMap &tris)
1073 {
1074  if (tris.size() == 0)
1075  {
1076  return;
1077  }
1078 
1079  string tag = "T";
1080 
1081  vector<MeshTri> triInfo;
1082 
1083  for (auto &i : tris)
1084  {
1085  MeshTri t;
1086  t.id = i.first;
1087  t.e[0] = i.second->GetEid(0);
1088  t.e[1] = i.second->GetEid(1);
1089  t.e[2] = i.second->GetEid(2);
1090  triInfo.push_back(t);
1091  }
1092 
1093  TiXmlElement *x = new TiXmlElement(tag);
1094  string triStr;
1096 
1097  x->SetAttribute("COMPRESSED",
1099  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1100 
1101  x->LinkEndChild(new TiXmlText(triStr));
1102 
1103  faceTag->LinkEndChild(x);
1104 }
1105 
1106 void MeshGraphXmlCompressed::WriteQuads(TiXmlElement *faceTag,
1107  QuadGeomMap &quads)
1108 {
1109  if (quads.size() == 0)
1110  {
1111  return;
1112  }
1113 
1114  string tag = "Q";
1115 
1116  vector<MeshQuad> quadInfo;
1117 
1118  for (auto &i : quads)
1119  {
1120  MeshQuad q;
1121  q.id = i.first;
1122  q.e[0] = i.second->GetEid(0);
1123  q.e[1] = i.second->GetEid(1);
1124  q.e[2] = i.second->GetEid(2);
1125  q.e[3] = i.second->GetEid(3);
1126  quadInfo.push_back(q);
1127  }
1128 
1129  TiXmlElement *x = new TiXmlElement(tag);
1130  string quadStr;
1132 
1133  x->SetAttribute("COMPRESSED",
1135  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1136 
1137  x->LinkEndChild(new TiXmlText(quadStr));
1138 
1139  faceTag->LinkEndChild(x);
1140 }
1141 
1142 void MeshGraphXmlCompressed::WriteHexs(TiXmlElement *elmtTag, HexGeomMap &hexs)
1143 {
1144  if (hexs.size() == 0)
1145  {
1146  return;
1147  }
1148 
1149  string tag = "H";
1150 
1151  vector<MeshHex> elementInfo;
1152 
1153  for (auto &i : hexs)
1154  {
1155  MeshHex e;
1156  e.id = i.first;
1157  e.f[0] = i.second->GetFid(0);
1158  e.f[1] = i.second->GetFid(1);
1159  e.f[2] = i.second->GetFid(2);
1160  e.f[3] = i.second->GetFid(3);
1161  e.f[4] = i.second->GetFid(4);
1162  e.f[5] = i.second->GetFid(5);
1163  elementInfo.push_back(e);
1164  }
1165 
1166  TiXmlElement *x = new TiXmlElement(tag);
1167  string elStr;
1169 
1170  x->SetAttribute("COMPRESSED",
1172  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1173 
1174  x->LinkEndChild(new TiXmlText(elStr));
1175 
1176  elmtTag->LinkEndChild(x);
1177 }
1178 
1179 void MeshGraphXmlCompressed::WritePrisms(TiXmlElement *elmtTag,
1180  PrismGeomMap &pris)
1181 {
1182  if (pris.size() == 0)
1183  {
1184  return;
1185  }
1186 
1187  string tag = "R";
1188 
1189  vector<MeshPrism> elementInfo;
1190 
1191  for (auto &i : pris)
1192  {
1193  MeshPrism e;
1194  e.id = i.first;
1195  e.f[0] = i.second->GetFid(0);
1196  e.f[1] = i.second->GetFid(1);
1197  e.f[2] = i.second->GetFid(2);
1198  e.f[3] = i.second->GetFid(3);
1199  e.f[4] = i.second->GetFid(4);
1200  elementInfo.push_back(e);
1201  }
1202 
1203  TiXmlElement *x = new TiXmlElement(tag);
1204  string elStr;
1206 
1207  x->SetAttribute("COMPRESSED",
1209  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1210 
1211  x->LinkEndChild(new TiXmlText(elStr));
1212 
1213  elmtTag->LinkEndChild(x);
1214 }
1215 
1216 void MeshGraphXmlCompressed::WritePyrs(TiXmlElement *elmtTag, PyrGeomMap &pyrs)
1217 {
1218  if (pyrs.size() == 0)
1219  {
1220  return;
1221  }
1222 
1223  string tag = "P";
1224 
1225  vector<MeshPyr> elementInfo;
1226 
1227  for (auto &i : pyrs)
1228  {
1229  MeshPyr e;
1230  e.id = i.first;
1231  e.f[0] = i.second->GetFid(0);
1232  e.f[1] = i.second->GetFid(1);
1233  e.f[2] = i.second->GetFid(2);
1234  e.f[3] = i.second->GetFid(3);
1235  e.f[4] = i.second->GetFid(4);
1236  elementInfo.push_back(e);
1237  }
1238 
1239  TiXmlElement *x = new TiXmlElement(tag);
1240  string elStr;
1242 
1243  x->SetAttribute("COMPRESSED",
1245  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1246 
1247  x->LinkEndChild(new TiXmlText(elStr));
1248 
1249  elmtTag->LinkEndChild(x);
1250 }
1251 
1252 void MeshGraphXmlCompressed::WriteTets(TiXmlElement *elmtTag, TetGeomMap &tets)
1253 {
1254  if (tets.size() == 0)
1255  {
1256  return;
1257  }
1258 
1259  string tag = "A";
1260 
1261  vector<MeshTet> elementInfo;
1262 
1263  for (auto &i : tets)
1264  {
1265  MeshTet e;
1266  e.id = i.first;
1267  e.f[0] = i.second->GetFid(0);
1268  e.f[1] = i.second->GetFid(1);
1269  e.f[2] = i.second->GetFid(2);
1270  e.f[3] = i.second->GetFid(3);
1271  elementInfo.push_back(e);
1272  }
1273 
1274  TiXmlElement *x = new TiXmlElement(tag);
1275  string elStr;
1277 
1278  x->SetAttribute("COMPRESSED",
1280  x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1281 
1282  x->LinkEndChild(new TiXmlText(elStr));
1283 
1284  elmtTag->LinkEndChild(x);
1285 }
1286 
1287 void MeshGraphXmlCompressed::WriteCurves(TiXmlElement *geomTag, CurveMap &edges,
1288  CurveMap &faces)
1289 {
1290  if (edges.size() == 0 && faces.size() == 0)
1291  {
1292  return;
1293  }
1294 
1295  TiXmlElement *curveTag = new TiXmlElement("CURVED");
1296 
1297  vector<MeshCurvedInfo> edgeInfo;
1298  vector<MeshCurvedInfo> faceInfo;
1299  MeshCurvedPts curvedPts;
1300  curvedPts.id = 0;
1301  int ptOffset = 0;
1302  int newIdx = 0;
1303  int edgeCnt = 0;
1304  int faceCnt = 0;
1305 
1306  for (auto &i : edges)
1307  {
1308  MeshCurvedInfo cinfo;
1309  cinfo.id = edgeCnt++;
1310  cinfo.entityid = i.first;
1311  cinfo.npoints = i.second->m_points.size();
1312  cinfo.ptype = i.second->m_ptype;
1313  cinfo.ptid = 0;
1314  cinfo.ptoffset = ptOffset;
1315 
1316  edgeInfo.push_back(cinfo);
1317 
1318  for (int j = 0; j < i.second->m_points.size(); j++)
1319  {
1320  MeshVertex v;
1321  v.id = newIdx;
1322  v.x = i.second->m_points[j]->x();
1323  v.y = i.second->m_points[j]->y();
1324  v.z = i.second->m_points[j]->z();
1325  curvedPts.pts.push_back(v);
1326  curvedPts.index.push_back(newIdx);
1327  newIdx++;
1328  }
1329  ptOffset += cinfo.npoints;
1330  }
1331 
1332  for (auto &i : faces)
1333  {
1334  MeshCurvedInfo cinfo;
1335  cinfo.id = faceCnt++;
1336  cinfo.entityid = i.first;
1337  cinfo.npoints = i.second->m_points.size();
1338  cinfo.ptype = i.second->m_ptype;
1339  cinfo.ptid = 0;
1340  cinfo.ptoffset = ptOffset;
1341 
1342  faceInfo.push_back(cinfo);
1343 
1344  for (int j = 0; j < i.second->m_points.size(); j++)
1345  {
1346  MeshVertex v;
1347  v.id = newIdx;
1348  v.x = i.second->m_points[j]->x();
1349  v.y = i.second->m_points[j]->y();
1350  v.z = i.second->m_points[j]->z();
1351  curvedPts.pts.push_back(v);
1352  curvedPts.index.push_back(newIdx);
1353  newIdx++;
1354  }
1355  ptOffset += cinfo.npoints;
1356  }
1357 
1358  curveTag->SetAttribute("COMPRESSED",
1360  curveTag->SetAttribute("BITSIZE",
1362 
1363  if (edgeInfo.size())
1364  {
1365  TiXmlElement *x = new TiXmlElement("E");
1366  string dataStr;
1368 
1369  x->LinkEndChild(new TiXmlText(dataStr));
1370  curveTag->LinkEndChild(x);
1371  }
1372 
1373  if (faceInfo.size())
1374  {
1375  TiXmlElement *x = new TiXmlElement("F");
1376  string dataStr;
1378 
1379  x->LinkEndChild(new TiXmlText(dataStr));
1380  curveTag->LinkEndChild(x);
1381  }
1382 
1383  if (edgeInfo.size() || faceInfo.size())
1384  {
1385  TiXmlElement *x = new TiXmlElement("DATAPOINTS");
1386  x->SetAttribute("ID", curvedPts.id);
1387  TiXmlElement *subx = new TiXmlElement("INDEX");
1388  string dataStr;
1390  dataStr);
1391  subx->LinkEndChild(new TiXmlText(dataStr));
1392  x->LinkEndChild(subx);
1393 
1394  subx = new TiXmlElement("POINTS");
1396  dataStr);
1397  subx->LinkEndChild(new TiXmlText(dataStr));
1398  x->LinkEndChild(subx);
1399  curveTag->LinkEndChild(x);
1400  }
1401 
1402  geomTag->LinkEndChild(curveTag);
1403 }
1404 } // namespace SpatialDomains
1405 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
std::vector< NekInt64 > index
id of this Point set
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
NekInt64 entityid
Id of this curved information.
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:88
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:62
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
STL namespace.
std::vector< MeshVertex > pts
mapping to access pts value.
std::map< int, PrismGeomSharedPtr > PrismGeomMap
Definition: PrismGeom.h:89
std::map< int, PyrGeomSharedPtr > PyrGeomMap
Definition: PyrGeom.h:81
NekInt64 ptype
point offset of data entry for this curve
std::map< int, TetGeomSharedPtr > TetGeomMap
Definition: TetGeom.h:90
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:59
NekInt64 ptoffset
the id of point data map (currently always 0 since we are using just one set).
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
std::int32_t NekInt
MeshGraphFactory & GetMeshGraphFactory()
Definition: MeshGraph.cpp:76
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:52
double NekDouble
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
NekInt64 npoints
The entity id corresponding to the global edge/curve.
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:77
int ZlibEncodeToBase64Str(std::vector< T > &in, std::string &out64)
Definition: CompressData.h:138
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:80
NekInt64 ptid
The number of points in this curved entity.
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:54
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:90
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:61
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:199
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:88
std::map< int, PointGeomSharedPtr > PointGeomMap
Definition: PointGeom.h:54
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::map< int, HexGeomSharedPtr > HexGeomMap
Definition: HexGeom.h:91
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:231