Nektar++
MeshGraphIOXmlCompressed.cpp
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////////////
2//
3// File: MeshGraphIOXmlCompressed.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
37
41
44
45// These are required for the Write(...) and Import(...) functions.
46#include <boost/algorithm/string/predicate.hpp>
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/format.hpp>
51#include <boost/iostreams/copy.hpp>
52#include <boost/iostreams/filter/zlib.hpp>
53#include <boost/iostreams/filtering_stream.hpp>
54#include <boost/make_shared.hpp>
55
56#include <tinyxml.h>
57using namespace std;
58
60{
61
64 "XmlCompressed", MeshGraphIOXmlCompressed::create,
65 "IO with compressed Xml geometry");
66
68{
69 PointGeomMap &vertSet = m_meshGraph->GetAllPointGeoms();
70 int spaceDimension = m_meshGraph->GetSpaceDimension();
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 spaceDimension, indx, xval, yval, zval));
194
195 vert->SetGlobalID(indx);
196 vertSet[indx] = vert;
197 }
198 }
199 else
200 {
201 ASSERTL0(false, "Compressed formats do not match. Expected :" +
203 " but got " + IsCompressed);
204 }
205}
206
208{
209 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
210 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
211 int spaceDimension = m_meshGraph->GetSpaceDimension();
212
213 // check to see if any scaling parameters are in
214 // attributes and determine these values
215 TiXmlElement *element = m_xmlGeom->FirstChildElement("VERTEX");
216 ASSERTL0(element, "Unable to find mesh VERTEX tag in file.");
217
218 NekDouble xscale, yscale, zscale;
219
220 LibUtilities::Interpreter expEvaluator;
221 const char *xscal = element->Attribute("XSCALE");
222 if (!xscal)
223 {
224 xscale = 1.0;
225 }
226 else
227 {
228 std::string xscalstr = xscal;
229 int expr_id = expEvaluator.DefineFunction("", xscalstr);
230 xscale = expEvaluator.Evaluate(expr_id);
231 }
232
233 const char *yscal = element->Attribute("YSCALE");
234 if (!yscal)
235 {
236 yscale = 1.0;
237 }
238 else
239 {
240 std::string yscalstr = yscal;
241 int expr_id = expEvaluator.DefineFunction("", yscalstr);
242 yscale = expEvaluator.Evaluate(expr_id);
243 }
244
245 const char *zscal = element->Attribute("ZSCALE");
246 if (!zscal)
247 {
248 zscale = 1.0;
249 }
250 else
251 {
252 std::string zscalstr = zscal;
253 int expr_id = expEvaluator.DefineFunction("", zscalstr);
254 zscale = expEvaluator.Evaluate(expr_id);
255 }
256
257 NekDouble xmove, ymove, zmove;
258
259 // check to see if any moving parameters are in
260 // attributes and determine these values
261
262 const char *xmov = element->Attribute("XMOVE");
263 if (!xmov)
264 {
265 xmove = 0.0;
266 }
267 else
268 {
269 std::string xmovstr = xmov;
270 int expr_id = expEvaluator.DefineFunction("", xmovstr);
271 xmove = expEvaluator.Evaluate(expr_id);
272 }
273
274 const char *ymov = element->Attribute("YMOVE");
275 if (!ymov)
276 {
277 ymove = 0.0;
278 }
279 else
280 {
281 std::string ymovstr = ymov;
282 int expr_id = expEvaluator.DefineFunction("", ymovstr);
283 ymove = expEvaluator.Evaluate(expr_id);
284 }
285
286 const char *zmov = element->Attribute("ZMOVE");
287 if (!zmov)
288 {
289 zmove = 0.0;
290 }
291 else
292 {
293 std::string zmovstr = zmov;
294 int expr_id = expEvaluator.DefineFunction("", zmovstr);
295 zmove = expEvaluator.Evaluate(expr_id);
296 }
297
298 /// Look for elements in CURVE block.
299 TiXmlElement *field = m_xmlGeom->FirstChildElement("CURVED");
300
301 if (!field) // return if no curved entities
302 {
303 return;
304 }
305
306 string IsCompressed;
307 field->QueryStringAttribute("COMPRESSED", &IsCompressed);
308
309 if (IsCompressed.size() == 0)
310 {
311 // this could be that the curved tag is empty
312 // in this case we dont want to read it
313 return;
314 }
315
316 ASSERTL0(boost::iequals(IsCompressed,
318 "Compressed formats do not match. Expected :" +
320 IsCompressed);
321
322 std::vector<SpatialDomains::MeshCurvedInfo> edginfo;
323 std::vector<SpatialDomains::MeshCurvedInfo> facinfo;
325
326 // read edge, face info and curved poitns.
327 TiXmlElement *x = field->FirstChildElement();
328 while (x)
329 {
330 const char *entitytype = x->Value();
331 // read in edge or face info
332 if (boost::iequals(entitytype, "E"))
333 {
334 // read in data
335 std::string elmtStr;
336 TiXmlNode *child = x->FirstChild();
337
338 if (child->Type() == TiXmlNode::TINYXML_TEXT)
339 {
340 elmtStr += child->ToText()->ValueStr();
341 }
342
344 edginfo);
345 }
346 else if (boost::iequals(entitytype, "F"))
347 {
348 // read in data
349 std::string elmtStr;
350 TiXmlNode *child = x->FirstChild();
351
352 if (child->Type() == TiXmlNode::TINYXML_TEXT)
353 {
354 elmtStr += child->ToText()->ValueStr();
355 }
356
358 facinfo);
359 }
360 else if (boost::iequals(entitytype, "DATAPOINTS"))
361 {
362 NekInt id;
363 ASSERTL0(x->Attribute("ID", &id),
364 "Failed to get ID from PTS section");
365 cpts.id = id;
366
367 // read in data
368 std::string elmtStr;
369
370 TiXmlElement *DataIdx = x->FirstChildElement("INDEX");
371 ASSERTL0(DataIdx, "Cannot read data index tag in compressed "
372 "curved section");
373
374 TiXmlNode *child = DataIdx->FirstChild();
375 if (child->Type() == TiXmlNode::TINYXML_TEXT)
376 {
377 elmtStr = child->ToText()->ValueStr();
378 }
379
381 cpts.index);
382
383 TiXmlElement *DataPts = x->FirstChildElement("POINTS");
384 ASSERTL0(DataPts, "Cannot read data pts tag in compressed "
385 "curved section");
386
387 child = DataPts->FirstChild();
388 if (child->Type() == TiXmlNode::TINYXML_TEXT)
389 {
390 elmtStr = child->ToText()->ValueStr();
391 }
392
394 cpts.pts);
395 }
396 else
397 {
398 ASSERTL0(false, "Unknown tag in curved section");
399 }
400 x = x->NextSiblingElement();
401 }
402
403 // rescale (x,y,z) points;
404 for (int i = 0; i < cpts.pts.size(); ++i)
405 {
406 cpts.pts[i].x = xscale * cpts.pts[i].x + xmove;
407 cpts.pts[i].y = yscale * cpts.pts[i].y + ymove;
408 cpts.pts[i].z = zscale * cpts.pts[i].z + zmove;
409 }
410
411 for (int i = 0; i < edginfo.size(); ++i)
412 {
413 int edgeid = edginfo[i].entityid;
415
417 edgeid, ptype = (LibUtilities::PointsType)edginfo[i].ptype));
418
419 // load points
420 int offset = edginfo[i].ptoffset;
421 for (int j = 0; j < edginfo[i].npoints; ++j)
422 {
423 int idx = cpts.index[offset + j];
424
426 spaceDimension, edginfo[i].id, cpts.pts[idx].x, cpts.pts[idx].y,
427 cpts.pts[idx].z));
428 curve->m_points.push_back(vert);
429 }
430
431 curvedEdges[edgeid] = curve;
432 }
433
434 for (int i = 0; i < facinfo.size(); ++i)
435 {
436 int faceid = facinfo[i].entityid;
438
440 faceid, ptype = (LibUtilities::PointsType)facinfo[i].ptype));
441
442 int offset = facinfo[i].ptoffset;
443 for (int j = 0; j < facinfo[i].npoints; ++j)
444 {
445 int idx = cpts.index[offset + j];
446
448 spaceDimension, facinfo[i].id, cpts.pts[idx].x, cpts.pts[idx].y,
449 cpts.pts[idx].z));
450 curve->m_points.push_back(vert);
451 }
452
453 curvedFaces[faceid] = curve;
454 }
455}
456
458{
459 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
460 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
461 int spaceDimension = m_meshGraph->GetSpaceDimension();
462
463 CurveMap::iterator it;
464
465 /// Look for elements in ELEMENT block.
466 TiXmlElement *field = m_xmlGeom->FirstChildElement("EDGE");
467
468 ASSERTL0(field, "Unable to find EDGE tag in file.");
469
470 string IsCompressed;
471 field->QueryStringAttribute("COMPRESSED", &IsCompressed);
472
473 ASSERTL0(boost::iequals(IsCompressed,
475 "Compressed formats do not match. Expected :" +
477 IsCompressed);
478 // Extract the edge body
479 TiXmlNode *edgeChild = field->FirstChild();
480 ASSERTL0(edgeChild, "Unable to extract the data from "
481 "the compressed edge tag.");
482
483 std::string edgeStr;
484 if (edgeChild->Type() == TiXmlNode::TINYXML_TEXT)
485 {
486 edgeStr += edgeChild->ToText()->ValueStr();
487 }
488
489 std::vector<SpatialDomains::MeshEdge> edgeData;
491
492 int indx;
493 for (int i = 0; i < edgeData.size(); ++i)
494 {
495 indx = edgeData[i].id;
496 PointGeomSharedPtr vertices[2] = {
497 m_meshGraph->GetVertex(edgeData[i].v0),
498 m_meshGraph->GetVertex(edgeData[i].v1)};
499 SegGeomSharedPtr edge;
500
501 it = curvedEdges.find(indx);
502 if (it == curvedEdges.end())
503 {
505 indx, spaceDimension, vertices);
506 }
507 else
508 {
510 indx, spaceDimension, vertices, it->second);
511 }
512 segGeoms[indx] = edge;
513 }
514}
515
517{
518 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
519 auto &triGeoms = m_meshGraph->GetAllTriGeoms();
520 auto &quadGeoms = m_meshGraph->GetAllQuadGeoms();
521
522 /// Look for elements in FACE block.
523 TiXmlElement *field = m_xmlGeom->FirstChildElement("FACE");
524
525 ASSERTL0(field, "Unable to find FACE tag in file.");
526
527 /// All faces are of the form: "<? ID="#"> ... </?>", with
528 /// ? being an element type (either Q or T).
529 /// They might be in compressed format and so then need upacking.
530
531 TiXmlElement *element = field->FirstChildElement();
532 CurveMap::iterator it;
533
534 while (element)
535 {
536 std::string elementType(element->ValueStr());
537
538 ASSERTL0(elementType == "Q" || elementType == "T",
539 (std::string("Unknown 3D face type: ") + elementType).c_str());
540
541 string IsCompressed;
542 element->QueryStringAttribute("COMPRESSED", &IsCompressed);
543
544 ASSERTL0(
545 boost::iequals(IsCompressed,
547 "Compressed formats do not match. Expected :" +
549 IsCompressed);
550
551 // Extract the face body
552 TiXmlNode *faceChild = element->FirstChild();
553 ASSERTL0(faceChild, "Unable to extract the data from "
554 "the compressed face tag.");
555
556 std::string faceStr;
557 if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
558 {
559 faceStr += faceChild->ToText()->ValueStr();
560 }
561
562 int indx;
563 if (elementType == "T")
564 {
565 std::vector<SpatialDomains::MeshTri> faceData;
567 faceData);
568
569 for (int i = 0; i < faceData.size(); ++i)
570 {
571 indx = faceData[i].id;
572
573 /// See if this face has curves.
574 it = curvedFaces.find(indx);
575
576 /// Create a TriGeom to hold the new definition.
578 m_meshGraph->GetSegGeom(faceData[i].e[0]),
579 m_meshGraph->GetSegGeom(faceData[i].e[1]),
580 m_meshGraph->GetSegGeom(faceData[i].e[2])};
581
582 TriGeomSharedPtr trigeom;
583 if (it == curvedFaces.end())
584 {
585 trigeom =
587 }
588 else
589 {
591 indx, edges, it->second);
592 }
593 trigeom->SetGlobalID(indx);
594 triGeoms[indx] = trigeom;
595 }
596 }
597 else if (elementType == "Q")
598 {
599 std::vector<SpatialDomains::MeshQuad> faceData;
601 faceData);
602
603 for (int i = 0; i < faceData.size(); ++i)
604 {
605 indx = faceData[i].id;
606
607 /// See if this face has curves.
608 it = curvedFaces.find(indx);
609
610 /// Create a QuadGeom to hold the new definition.
612 m_meshGraph->GetSegGeom(faceData[i].e[0]),
613 m_meshGraph->GetSegGeom(faceData[i].e[1]),
614 m_meshGraph->GetSegGeom(faceData[i].e[2]),
615 m_meshGraph->GetSegGeom(faceData[i].e[3])};
616
617 QuadGeomSharedPtr quadgeom;
618 if (it == curvedFaces.end())
619 {
620 quadgeom =
622 }
623 else
624 {
626 indx, edges, it->second);
627 }
628 quadgeom->SetGlobalID(indx);
629 quadGeoms[indx] = quadgeom;
630 }
631 }
632 /// Keep looking
633 element = element->NextSiblingElement();
634 }
635}
636
638{
639 auto &curvedEdges = m_meshGraph->GetCurvedEdges();
640 auto &segGeoms = m_meshGraph->GetAllSegGeoms();
641 int spaceDimension = m_meshGraph->GetSpaceDimension();
642
643 TiXmlElement *field = nullptr;
644
645 /// Look for elements in ELEMENT block.
646 field = m_xmlGeom->FirstChildElement("ELEMENT");
647
648 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
649
650 /// All elements are of the form: "<S ID = n> ... </S>", with
651 /// ? being the element type.
652
653 TiXmlElement *segment = field->FirstChildElement("S");
654 CurveMap::iterator it;
655
656 while (segment)
657 {
658 string IsCompressed;
659 segment->QueryStringAttribute("COMPRESSED", &IsCompressed);
660 ASSERTL0(
661 boost::iequals(IsCompressed,
663 "Compressed formats do not match. Expected :" +
665 IsCompressed);
666
667 // Extract the face body
668 TiXmlNode *child = segment->FirstChild();
669 ASSERTL0(child, "Unable to extract the data from "
670 "the compressed face tag.");
671
672 std::string str;
673 if (child->Type() == TiXmlNode::TINYXML_TEXT)
674 {
675 str += child->ToText()->ValueStr();
676 }
677
678 int indx;
679
680 std::vector<SpatialDomains::MeshEdge> data;
682
683 for (int i = 0; i < data.size(); ++i)
684 {
685 indx = data[i].id;
686
687 /// See if this face has curves.
688 it = curvedEdges.find(indx);
689
690 PointGeomSharedPtr vertices[2] = {
691 m_meshGraph->GetVertex(data[i].v0),
692 m_meshGraph->GetVertex(data[i].v1)};
694
695 if (it == curvedEdges.end())
696 {
698 indx, spaceDimension, vertices);
699 seg->SetGlobalID(indx); // Set global mesh id
700 }
701 else
702 {
704 indx, spaceDimension, vertices, it->second);
705 seg->SetGlobalID(indx); // Set global mesh id
706 }
707 seg->SetGlobalID(indx);
708 segGeoms[indx] = seg;
709 }
710 /// Keep looking for additional segments
711 segment = segment->NextSiblingElement("S");
712 }
713}
714
716{
717 auto &curvedFaces = m_meshGraph->GetCurvedFaces();
718 auto &triGeoms = m_meshGraph->GetAllTriGeoms();
719 auto &quadGeoms = m_meshGraph->GetAllQuadGeoms();
720
721 /// Look for elements in ELEMENT block.
722 TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
723
724 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
725
726 // Set up curve map for curved elements on an embedded manifold.
727 CurveMap::iterator it;
728
729 /// All elements are of the form: "<? ID="#"> ... </?>", with
730 /// ? being the element type.
731
732 TiXmlElement *element = field->FirstChildElement();
733
734 while (element)
735 {
736 std::string elementType(element->ValueStr());
737
738 ASSERTL0(
739 elementType == "Q" || elementType == "T",
740 (std::string("Unknown 2D element type: ") + elementType).c_str());
741
742 string IsCompressed;
743 element->QueryStringAttribute("COMPRESSED", &IsCompressed);
744
745 ASSERTL0(
746 boost::iequals(IsCompressed,
748 "Compressed formats do not match. Expected :" +
750 IsCompressed);
751
752 // Extract the face body
753 TiXmlNode *faceChild = element->FirstChild();
754 ASSERTL0(faceChild, "Unable to extract the data from "
755 "the compressed face tag.");
756
757 std::string faceStr;
758 if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
759 {
760 faceStr += faceChild->ToText()->ValueStr();
761 }
762
763 int indx;
764 if (elementType == "T")
765 {
766 std::vector<SpatialDomains::MeshTri> faceData;
768 faceData);
769
770 for (int i = 0; i < faceData.size(); ++i)
771 {
772 indx = faceData[i].id;
773
774 /// See if this face has curves.
775 it = curvedFaces.find(indx);
776
777 /// Create a TriGeom to hold the new definition.
779 m_meshGraph->GetSegGeom(faceData[i].e[0]),
780 m_meshGraph->GetSegGeom(faceData[i].e[1]),
781 m_meshGraph->GetSegGeom(faceData[i].e[2])};
782
783 TriGeomSharedPtr trigeom;
784 if (it == curvedFaces.end())
785 {
786 trigeom =
788 }
789 else
790 {
792 indx, edges, it->second);
793 }
794 trigeom->SetGlobalID(indx);
795 triGeoms[indx] = trigeom;
796 }
797 }
798 else if (elementType == "Q")
799 {
800 std::vector<SpatialDomains::MeshQuad> faceData;
802 faceData);
803
804 for (int i = 0; i < faceData.size(); ++i)
805 {
806 indx = faceData[i].id;
807
808 /// See if this face has curves.
809 it = curvedFaces.find(indx);
810
811 /// Create a QuadGeom to hold the new definition.
813 m_meshGraph->GetSegGeom(faceData[i].e[0]),
814 m_meshGraph->GetSegGeom(faceData[i].e[1]),
815 m_meshGraph->GetSegGeom(faceData[i].e[2]),
816 m_meshGraph->GetSegGeom(faceData[i].e[3])};
817
818 QuadGeomSharedPtr quadgeom;
819 if (it == curvedFaces.end())
820 {
821 quadgeom =
823 }
824 else
825 {
827 indx, edges, it->second);
828 }
829 quadgeom->SetGlobalID(indx);
830 quadGeoms[indx] = quadgeom;
831 }
832 }
833 /// Keep looking
834 element = element->NextSiblingElement();
835 }
836}
837
839{
840 auto &tetGeoms = m_meshGraph->GetAllTetGeoms();
841 auto &pyrGeoms = m_meshGraph->GetAllPyrGeoms();
842 auto &prismGeoms = m_meshGraph->GetAllPrismGeoms();
843 auto &hexGeoms = m_meshGraph->GetAllHexGeoms();
844
845 /// Look for elements in ELEMENT block.
846 TiXmlElement *field = m_xmlGeom->FirstChildElement("ELEMENT");
847
848 ASSERTL0(field, "Unable to find ELEMENT tag in file.");
849
850 /// All elements are of the form: "<? ID="#"> ... </?>", with
851 /// ? being the element type.
852
853 TiXmlElement *element = field->FirstChildElement();
854
855 while (element)
856 {
857 std::string elementType(element->ValueStr());
858
859 // A - tet, P - pyramid, R - prism, H - hex
860 ASSERTL0(
861 elementType == "A" || elementType == "P" || elementType == "R" ||
862 elementType == "H",
863 (std::string("Unknown 3D element type: ") + elementType).c_str());
864
865 string IsCompressed;
866 element->QueryStringAttribute("COMPRESSED", &IsCompressed);
867
868 ASSERTL0(
869 boost::iequals(IsCompressed,
871 "Compressed formats do not match. Expected :" +
873 IsCompressed);
874
875 // Extract the face body
876 TiXmlNode *child = element->FirstChild();
877 ASSERTL0(child, "Unable to extract the data from "
878 "the compressed face tag.");
879
880 std::string str;
881 if (child->Type() == TiXmlNode::TINYXML_TEXT)
882 {
883 str += child->ToText()->ValueStr();
884 }
885
886 int indx;
887 if (elementType == "A")
888 {
889 std::vector<SpatialDomains::MeshTet> data;
891 TriGeomSharedPtr tfaces[4];
892 for (int i = 0; i < data.size(); ++i)
893 {
894 indx = data[i].id;
895 for (int j = 0; j < 4; ++j)
896 {
898 m_meshGraph->GetGeometry2D(data[i].f[j]);
899 tfaces[j] = static_pointer_cast<TriGeom>(face);
900 }
901
902 TetGeomSharedPtr tetgeom(
904 tetGeoms[indx] = tetgeom;
905 m_meshGraph->PopulateFaceToElMap(tetgeom, 4);
906 }
907 }
908 else if (elementType == "P")
909 {
910 std::vector<SpatialDomains::MeshPyr> data;
912 Geometry2DSharedPtr faces[5];
913 for (int i = 0; i < data.size(); ++i)
914 {
915 indx = data[i].id;
916 int Ntfaces = 0;
917 int Nqfaces = 0;
918 for (int j = 0; j < 5; ++j)
919 {
921 m_meshGraph->GetGeometry2D(data[i].f[j]);
922
923 if (face == Geometry2DSharedPtr() ||
924 (face->GetShapeType() != LibUtilities::eTriangle &&
925 face->GetShapeType() != LibUtilities::eQuadrilateral))
926 {
927 std::stringstream errorstring;
928 errorstring << "Element " << indx
929 << " has invalid face: " << j;
930 ASSERTL0(false, errorstring.str().c_str());
931 }
932 else if (face->GetShapeType() == LibUtilities::eTriangle)
933 {
934 faces[j] = static_pointer_cast<TriGeom>(face);
935 Ntfaces++;
936 }
937 else if (face->GetShapeType() ==
939 {
940 faces[j] = static_pointer_cast<QuadGeom>(face);
941 Nqfaces++;
942 }
943 }
944 ASSERTL0((Ntfaces == 4) && (Nqfaces == 1),
945 "Did not identify the correct number of "
946 "triangular and quadrilateral faces for a "
947 "pyramid");
948
949 PyrGeomSharedPtr pyrgeom(
951
952 pyrGeoms[indx] = pyrgeom;
953 m_meshGraph->PopulateFaceToElMap(pyrgeom, 5);
954 }
955 }
956 else if (elementType == "R")
957 {
958 std::vector<SpatialDomains::MeshPrism> data;
960 Geometry2DSharedPtr faces[5];
961 for (int i = 0; i < data.size(); ++i)
962 {
963 indx = data[i].id;
964 int Ntfaces = 0;
965 int Nqfaces = 0;
966 for (int j = 0; j < 5; ++j)
967 {
969 m_meshGraph->GetGeometry2D(data[i].f[j]);
970 if (face == Geometry2DSharedPtr() ||
971 (face->GetShapeType() != LibUtilities::eTriangle &&
972 face->GetShapeType() != LibUtilities::eQuadrilateral))
973 {
974 std::stringstream errorstring;
975 errorstring << "Element " << indx
976 << " has invalid face: " << j;
977 ASSERTL0(false, errorstring.str().c_str());
978 }
979 else if (face->GetShapeType() == LibUtilities::eTriangle)
980 {
981 faces[j] = static_pointer_cast<TriGeom>(face);
982 Ntfaces++;
983 }
984 else if (face->GetShapeType() ==
986 {
987 faces[j] = static_pointer_cast<QuadGeom>(face);
988 Nqfaces++;
989 }
990 }
991 ASSERTL0((Ntfaces == 2) && (Nqfaces == 3),
992 "Did not identify the correct number of "
993 "triangular and quadrilateral faces for a "
994 "prism");
995
996 PrismGeomSharedPtr prismgeom(
998
999 prismGeoms[indx] = prismgeom;
1000 m_meshGraph->PopulateFaceToElMap(prismgeom, 5);
1001 }
1002 }
1003 else if (elementType == "H")
1004 {
1005 std::vector<SpatialDomains::MeshHex> data;
1007
1008 QuadGeomSharedPtr faces[6];
1009 for (int i = 0; i < data.size(); ++i)
1010 {
1011 indx = data[i].id;
1012 for (int j = 0; j < 6; ++j)
1013 {
1014 Geometry2DSharedPtr face =
1015 m_meshGraph->GetGeometry2D(data[i].f[j]);
1016 faces[j] = static_pointer_cast<QuadGeom>(face);
1017 }
1018
1019 HexGeomSharedPtr hexgeom(
1021 hexGeoms[indx] = hexgeom;
1022 m_meshGraph->PopulateFaceToElMap(hexgeom, 6);
1023 }
1024 }
1025 /// Keep looking
1026 element = element->NextSiblingElement();
1027 }
1028}
1029
1031 PointGeomMap &verts)
1032{
1033 if (verts.size() == 0)
1034 {
1035 return;
1036 }
1037
1038 TiXmlElement *vertTag = new TiXmlElement("VERTEX");
1039
1040 vector<MeshVertex> vertInfo;
1041
1042 for (auto &i : verts)
1043 {
1044 MeshVertex v;
1045 v.id = i.first;
1046 v.x = i.second->x();
1047 v.y = i.second->y();
1048 v.z = i.second->z();
1049 vertInfo.push_back(v);
1050 }
1051
1052 vertTag->SetAttribute("COMPRESSED",
1054 vertTag->SetAttribute("BITSIZE",
1056
1057 string vertStr;
1059
1060 vertTag->LinkEndChild(new TiXmlText(vertStr));
1061
1062 geomTag->LinkEndChild(vertTag);
1063}
1064
1066 SegGeomMap &edges)
1067{
1068 int meshDimension = m_meshGraph->GetMeshDimension();
1069
1070 if (edges.size() == 0)
1071 {
1072 return;
1073 }
1074
1075 TiXmlElement *edgeTag = new TiXmlElement(meshDimension == 1 ? "S" : "EDGE");
1076
1077 vector<MeshEdge> edgeInfo;
1078
1079 for (auto &i : edges)
1080 {
1081 MeshEdge e;
1082 e.id = i.first;
1083 e.v0 = i.second->GetVid(0);
1084 e.v1 = i.second->GetVid(1);
1085 edgeInfo.push_back(e);
1086 }
1087
1088 string edgeStr;
1090
1091 edgeTag->SetAttribute("COMPRESSED",
1093 edgeTag->SetAttribute("BITSIZE",
1095
1096 edgeTag->LinkEndChild(new TiXmlText(edgeStr));
1097
1098 if (meshDimension == 1)
1099 {
1100 TiXmlElement *tmp = new TiXmlElement("ELEMENT");
1101 tmp->LinkEndChild(edgeTag);
1102 geomTag->LinkEndChild(tmp);
1103 }
1104 else
1105 {
1106 geomTag->LinkEndChild(edgeTag);
1107 }
1108}
1109
1110void MeshGraphIOXmlCompressed::v_WriteTris(TiXmlElement *faceTag,
1111 TriGeomMap &tris)
1112{
1113 if (tris.size() == 0)
1114 {
1115 return;
1116 }
1117
1118 string tag = "T";
1119
1120 vector<MeshTri> triInfo;
1121
1122 for (auto &i : tris)
1123 {
1124 MeshTri t;
1125 t.id = i.first;
1126 t.e[0] = i.second->GetEid(0);
1127 t.e[1] = i.second->GetEid(1);
1128 t.e[2] = i.second->GetEid(2);
1129 triInfo.push_back(t);
1130 }
1131
1132 TiXmlElement *x = new TiXmlElement(tag);
1133 string triStr;
1135
1136 x->SetAttribute("COMPRESSED",
1138 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1139
1140 x->LinkEndChild(new TiXmlText(triStr));
1141
1142 faceTag->LinkEndChild(x);
1143}
1144
1146 QuadGeomMap &quads)
1147{
1148 if (quads.size() == 0)
1149 {
1150 return;
1151 }
1152
1153 string tag = "Q";
1154
1155 vector<MeshQuad> quadInfo;
1156
1157 for (auto &i : quads)
1158 {
1159 MeshQuad q;
1160 q.id = i.first;
1161 q.e[0] = i.second->GetEid(0);
1162 q.e[1] = i.second->GetEid(1);
1163 q.e[2] = i.second->GetEid(2);
1164 q.e[3] = i.second->GetEid(3);
1165 quadInfo.push_back(q);
1166 }
1167
1168 TiXmlElement *x = new TiXmlElement(tag);
1169 string quadStr;
1171
1172 x->SetAttribute("COMPRESSED",
1174 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1175
1176 x->LinkEndChild(new TiXmlText(quadStr));
1177
1178 faceTag->LinkEndChild(x);
1179}
1180
1181void MeshGraphIOXmlCompressed::v_WriteHexs(TiXmlElement *elmtTag,
1182 HexGeomMap &hexs)
1183{
1184 if (hexs.size() == 0)
1185 {
1186 return;
1187 }
1188
1189 string tag = "H";
1190
1191 vector<MeshHex> elementInfo;
1192
1193 for (auto &i : hexs)
1194 {
1195 MeshHex e;
1196 e.id = i.first;
1197 e.f[0] = i.second->GetFid(0);
1198 e.f[1] = i.second->GetFid(1);
1199 e.f[2] = i.second->GetFid(2);
1200 e.f[3] = i.second->GetFid(3);
1201 e.f[4] = i.second->GetFid(4);
1202 e.f[5] = i.second->GetFid(5);
1203 elementInfo.push_back(e);
1204 }
1205
1206 TiXmlElement *x = new TiXmlElement(tag);
1207 string elStr;
1209
1210 x->SetAttribute("COMPRESSED",
1212 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1213
1214 x->LinkEndChild(new TiXmlText(elStr));
1215
1216 elmtTag->LinkEndChild(x);
1217}
1218
1220 PrismGeomMap &pris)
1221{
1222 if (pris.size() == 0)
1223 {
1224 return;
1225 }
1226
1227 string tag = "R";
1228
1229 vector<MeshPrism> elementInfo;
1230
1231 for (auto &i : pris)
1232 {
1233 MeshPrism e;
1234 e.id = i.first;
1235 e.f[0] = i.second->GetFid(0);
1236 e.f[1] = i.second->GetFid(1);
1237 e.f[2] = i.second->GetFid(2);
1238 e.f[3] = i.second->GetFid(3);
1239 e.f[4] = i.second->GetFid(4);
1240 elementInfo.push_back(e);
1241 }
1242
1243 TiXmlElement *x = new TiXmlElement(tag);
1244 string elStr;
1246
1247 x->SetAttribute("COMPRESSED",
1249 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1250
1251 x->LinkEndChild(new TiXmlText(elStr));
1252
1253 elmtTag->LinkEndChild(x);
1254}
1255
1256void MeshGraphIOXmlCompressed::v_WritePyrs(TiXmlElement *elmtTag,
1257 PyrGeomMap &pyrs)
1258{
1259 if (pyrs.size() == 0)
1260 {
1261 return;
1262 }
1263
1264 string tag = "P";
1265
1266 vector<MeshPyr> elementInfo;
1267
1268 for (auto &i : pyrs)
1269 {
1270 MeshPyr e;
1271 e.id = i.first;
1272 e.f[0] = i.second->GetFid(0);
1273 e.f[1] = i.second->GetFid(1);
1274 e.f[2] = i.second->GetFid(2);
1275 e.f[3] = i.second->GetFid(3);
1276 e.f[4] = i.second->GetFid(4);
1277 elementInfo.push_back(e);
1278 }
1279
1280 TiXmlElement *x = new TiXmlElement(tag);
1281 string elStr;
1283
1284 x->SetAttribute("COMPRESSED",
1286 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1287
1288 x->LinkEndChild(new TiXmlText(elStr));
1289
1290 elmtTag->LinkEndChild(x);
1291}
1292
1293void MeshGraphIOXmlCompressed::v_WriteTets(TiXmlElement *elmtTag,
1294 TetGeomMap &tets)
1295{
1296 if (tets.size() == 0)
1297 {
1298 return;
1299 }
1300
1301 string tag = "A";
1302
1303 vector<MeshTet> elementInfo;
1304
1305 for (auto &i : tets)
1306 {
1307 MeshTet e;
1308 e.id = i.first;
1309 e.f[0] = i.second->GetFid(0);
1310 e.f[1] = i.second->GetFid(1);
1311 e.f[2] = i.second->GetFid(2);
1312 e.f[3] = i.second->GetFid(3);
1313 elementInfo.push_back(e);
1314 }
1315
1316 TiXmlElement *x = new TiXmlElement(tag);
1317 string elStr;
1319
1320 x->SetAttribute("COMPRESSED",
1322 x->SetAttribute("BITSIZE", LibUtilities::CompressData::GetBitSizeStr());
1323
1324 x->LinkEndChild(new TiXmlText(elStr));
1325
1326 elmtTag->LinkEndChild(x);
1327}
1328
1330 CurveMap &edges, CurveMap &faces)
1331{
1332 if (edges.size() == 0 && faces.size() == 0)
1333 {
1334 return;
1335 }
1336
1337 TiXmlElement *curveTag = new TiXmlElement("CURVED");
1338
1339 vector<MeshCurvedInfo> edgeInfo;
1340 vector<MeshCurvedInfo> faceInfo;
1341 MeshCurvedPts curvedPts;
1342 curvedPts.id = 0;
1343 int ptOffset = 0;
1344 int newIdx = 0;
1345 int edgeCnt = 0;
1346 int faceCnt = 0;
1347
1348 for (auto &i : edges)
1349 {
1350 MeshCurvedInfo cinfo;
1351 cinfo.id = edgeCnt++;
1352 cinfo.entityid = i.first;
1353 cinfo.npoints = i.second->m_points.size();
1354 cinfo.ptype = i.second->m_ptype;
1355 cinfo.ptid = 0;
1356 cinfo.ptoffset = ptOffset;
1357
1358 edgeInfo.push_back(cinfo);
1359
1360 for (int j = 0; j < i.second->m_points.size(); j++)
1361 {
1362 MeshVertex v;
1363 v.id = newIdx;
1364 v.x = i.second->m_points[j]->x();
1365 v.y = i.second->m_points[j]->y();
1366 v.z = i.second->m_points[j]->z();
1367 curvedPts.pts.push_back(v);
1368 curvedPts.index.push_back(newIdx);
1369 newIdx++;
1370 }
1371 ptOffset += cinfo.npoints;
1372 }
1373
1374 for (auto &i : faces)
1375 {
1376 MeshCurvedInfo cinfo;
1377 cinfo.id = faceCnt++;
1378 cinfo.entityid = i.first;
1379 cinfo.npoints = i.second->m_points.size();
1380 cinfo.ptype = i.second->m_ptype;
1381 cinfo.ptid = 0;
1382 cinfo.ptoffset = ptOffset;
1383
1384 faceInfo.push_back(cinfo);
1385
1386 for (int j = 0; j < i.second->m_points.size(); j++)
1387 {
1388 MeshVertex v;
1389 v.id = newIdx;
1390 v.x = i.second->m_points[j]->x();
1391 v.y = i.second->m_points[j]->y();
1392 v.z = i.second->m_points[j]->z();
1393 curvedPts.pts.push_back(v);
1394 curvedPts.index.push_back(newIdx);
1395 newIdx++;
1396 }
1397 ptOffset += cinfo.npoints;
1398 }
1399
1400 curveTag->SetAttribute("COMPRESSED",
1402 curveTag->SetAttribute("BITSIZE",
1404
1405 if (edgeInfo.size())
1406 {
1407 TiXmlElement *x = new TiXmlElement("E");
1408 string dataStr;
1410
1411 x->LinkEndChild(new TiXmlText(dataStr));
1412 curveTag->LinkEndChild(x);
1413 }
1414
1415 if (faceInfo.size())
1416 {
1417 TiXmlElement *x = new TiXmlElement("F");
1418 string dataStr;
1420
1421 x->LinkEndChild(new TiXmlText(dataStr));
1422 curveTag->LinkEndChild(x);
1423 }
1424
1425 if (edgeInfo.size() || faceInfo.size())
1426 {
1427 TiXmlElement *x = new TiXmlElement("DATAPOINTS");
1428 x->SetAttribute("ID", curvedPts.id);
1429 TiXmlElement *subx = new TiXmlElement("INDEX");
1430 string dataStr;
1432 dataStr);
1433 subx->LinkEndChild(new TiXmlText(dataStr));
1434 x->LinkEndChild(subx);
1435
1436 subx = new TiXmlElement("POINTS");
1438 dataStr);
1439 subx->LinkEndChild(new TiXmlText(dataStr));
1440 x->LinkEndChild(subx);
1441 curveTag->LinkEndChild(x);
1442 }
1443
1444 geomTag->LinkEndChild(curveTag);
1445}
1446} // namespace Nektar::SpatialDomains
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:76
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void v_WriteEdges(TiXmlElement *geomTag, SegGeomMap &edges) override
void v_WriteHexs(TiXmlElement *elmtTag, HexGeomMap &hexs) override
void v_WriteVertices(TiXmlElement *geomTag, PointGeomMap &verts) override
void v_WritePyrs(TiXmlElement *elmtTag, PyrGeomMap &pyrs) override
void v_WriteTets(TiXmlElement *elmtTag, TetGeomMap &tets) override
void v_WriteCurves(TiXmlElement *geomTag, CurveMap &edges, CurveMap &faces) override
void v_WritePrisms(TiXmlElement *elmtTag, PrismGeomMap &pris) override
void v_WriteTris(TiXmlElement *faceTag, TriGeomMap &tris) override
void v_WriteQuads(TiXmlElement *faceTag, QuadGeomMap &quads) override
static const int kNedges
Definition: QuadGeom.h:74
static const int kNedges
Get the orientation of face1.
Definition: TriGeom.h:70
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:216
int ZlibEncodeToBase64Str(std::vector< T > &in, std::string &out64)
Definition: CompressData.h:129
std::map< int, TriGeomSharedPtr > TriGeomMap
Definition: TriGeom.h:57
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:45
std::map< int, PyrGeomSharedPtr > PyrGeomMap
Definition: PyrGeom.h:76
std::map< int, QuadGeomSharedPtr > QuadGeomMap
Definition: QuadGeom.h:52
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:82
std::map< int, SegGeomSharedPtr > SegGeomMap
Definition: SegGeom.h:49
std::shared_ptr< Curve > CurveSharedPtr
Definition: Curve.hpp:58
std::unordered_map< int, CurveSharedPtr > CurveMap
Definition: Curve.hpp:59
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:59
std::map< int, TetGeomSharedPtr > TetGeomMap
Definition: TetGeom.h:84
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:75
MeshGraphIOFactory & GetMeshGraphIOFactory()
Definition: MeshGraphIO.cpp:47
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::map< int, PrismGeomSharedPtr > PrismGeomMap
Definition: PrismGeom.h:83
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:57
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:62
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:56
std::map< int, HexGeomSharedPtr > HexGeomMap
Definition: HexGeom.h:85
std::map< int, PointGeomSharedPtr > PointGeomMap
Definition: PointGeom.h:51
std::vector< double > q(NPUPPER *NPUPPER)
std::int32_t NekInt
double NekDouble
STL namespace.
NekInt64 ptid
The number of points in this curved entity.
NekInt64 npoints
The entity id corresponding to the global edge/curve.
NekInt64 ptype
point offset of data entry for this curve
NekInt64 ptoffset
the id of point data map (currently always 0 since we are using just one set).
NekInt64 entityid
Id of this curved information.
std::vector< MeshVertex > pts
mapping to access pts value.
std::vector< NekInt64 > index
id of this Point set