Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MeshPartition.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MeshPartition.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:
33 //
34 //
35 ////////////////////////////////////////////////////////////////////////////////
36 
37 #ifndef TIXML_USE_STL
38 #define TIXML_USE_STL
39 #endif
40 
41 #include <LibUtilities/BasicUtils/MeshPartition.h>
42 
43 #include <iostream>
44 #include <iomanip>
45 #include <vector>
46 #include <map>
47 
48 #include <tinyxml.h>
49 
50 #include <LibUtilities/BasicUtils/ParseUtils.hpp>
51 #include <LibUtilities/BasicUtils/SessionReader.h>
52 #include <LibUtilities/BasicUtils/ShapeType.hpp>
53 #include <LibUtilities/BasicUtils/FileSystem.h>
54 #include <LibUtilities/BasicUtils/CompressData.h>
55 
56 #include <LibUtilities/Foundations/Foundations.hpp>
57 
58 
59 #include <boost/algorithm/string.hpp>
60 #include <boost/graph/adjacency_list.hpp>
61 #include <boost/graph/adjacency_iterator.hpp>
62 #include <boost/graph/detail/edge.hpp>
63 #include <boost/format.hpp>
64 
65 namespace Nektar
66 {
67  namespace LibUtilities
68  {
69  MeshPartitionFactory& GetMeshPartitionFactory()
70  {
71  typedef Loki::SingletonHolder<MeshPartitionFactory,
72  Loki::CreateUsingNew,
73  Loki::NoDestroy,
74  Loki::SingleThreaded> Type;
75  return Type::Instance();
76  }
77 
78  MeshPartition::MeshPartition(const LibUtilities::SessionReaderSharedPtr& pSession) :
79  m_isCompressed(false),
80  m_numFields(0),
81  m_fieldNameToId(),
82  m_comm(pSession->GetComm()),
83  m_weightingRequired(false)
84  {
85  ReadConditions(pSession);
86  ReadGeometry(pSession);
87  ReadExpansions(pSession);
88  }
89 
90  MeshPartition::~MeshPartition()
91  {
92 
93  }
94 
95  void MeshPartition::PartitionMesh(int nParts, bool shared,
96  bool overlapping)
97  {
98  ASSERTL0(m_meshElements.size() >= nParts,
99  "Too few elements for this many processes.");
100  m_shared = shared;
101 
102  if (m_weightingRequired)
103  {
104  WeightElements();
105  }
106  CreateGraph(m_mesh);
107  PartitionGraph(m_mesh, nParts, m_localPartition, overlapping);
108  }
109 
110  void MeshPartition::WriteLocalPartition(LibUtilities::SessionReaderSharedPtr& pSession)
111  {
112  TiXmlDocument vNew;
113  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
114  vNew.LinkEndChild(decl);
115 
116  TiXmlElement* vElmtNektar;
117  vElmtNektar = new TiXmlElement("NEKTAR");
118 
119  int rank = m_comm->GetRowComm()->GetRank();
120  OutputPartition(pSession, m_localPartition[rank], vElmtNektar);
121 
122  vNew.LinkEndChild(vElmtNektar);
123 
124  std::string dirname = pSession->GetSessionName() + "_xml";
125  fs::path pdirname(dirname);
126 
127  boost::format pad("P%1$07d.xml");
128  pad % rank;
129  fs::path pFilename(pad.str());
130 
131  if(!fs::is_directory(dirname))
132  {
133  fs::create_directory(dirname);
134  }
135 
136  fs::path fullpath = pdirname / pFilename;
137  vNew.SaveFile(PortablePath(fullpath));
138  }
139 
140  void MeshPartition::WriteAllPartitions(LibUtilities::SessionReaderSharedPtr& pSession)
141  {
142  for (int i = 0; i < m_localPartition.size(); ++i)
143  {
144  TiXmlDocument vNew;
145  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
146  vNew.LinkEndChild(decl);
147 
148  TiXmlElement* vElmtNektar;
149  vElmtNektar = new TiXmlElement("NEKTAR");
150 
151  OutputPartition(pSession, m_localPartition[i], vElmtNektar);
152 
153  vNew.LinkEndChild(vElmtNektar);
154 
155  std::string dirname = pSession->GetSessionName() + "_xml";
156  fs::path pdirname(dirname);
157 
158  boost::format pad("P%1$07d.xml");
159  pad % i;
160  fs::path pFilename(pad.str());
161 
162  fs::path fullpath = pdirname / pFilename;
163 
164  if(!fs::is_directory(dirname))
165  {
166  fs::create_directory(dirname);
167  }
168 
169  vNew.SaveFile(PortablePath(fullpath));
170  }
171  }
172 
173  void MeshPartition::GetCompositeOrdering(CompositeOrdering &composites)
174  {
175  std::map<int, MeshEntity>::iterator it;
176  for (it = m_meshComposites.begin();
177  it != m_meshComposites.end(); ++it)
178  {
179  composites[it->first] = it->second.list;
180  }
181  }
182 
183  void MeshPartition::GetBndRegionOrdering(BndRegionOrdering &bndRegs)
184  {
185  bndRegs = m_bndRegOrder;
186  }
187 
188 
189  void MeshPartition::ReadExpansions(const LibUtilities::SessionReaderSharedPtr& pSession)
190  {
191  // Find the Expansions tag
192  TiXmlElement *expansionTypes = pSession->GetElement("Nektar/Expansions");
193 
194  // Find the Expansion type
195  TiXmlElement *expansion = expansionTypes->FirstChildElement();
196  std::string expType = expansion->Value();
197 
198 
199  if(expType != "E")
200  {
201  ASSERTL0(false,"Expansion type not defined or not supported at the moment");
202  }
203 
204  /// Expansiontypes will contain plenty of data,
205  /// where relevant at this stage are composite
206  /// ID(s) that this expansion type describes,
207  /// nummodes and a list of fields that this
208  /// expansion relates to. If this does not exist
209  /// the variable is only set to "DefaultVar".
210 
211  while (expansion)
212  {
213  std::vector<unsigned int> composite;
214  std::vector<unsigned int> nummodes;
215  std::vector<std::string> fieldName;
216 
217  const char *nModesStr = expansion->Attribute("NUMMODES");
218  ASSERTL0(nModesStr,"NUMMODES was not defined in EXPANSION section of input");
219  std::string numModesStr = nModesStr;
220  bool valid = ParseUtils::GenerateOrderedVector(numModesStr.c_str(), nummodes);
221  ASSERTL0(valid, "Unable to correctly parse the number of modes.");
222 
223  if (nummodes.size() == 1)
224  {
225  for (int i = 1; i < m_dim; i++)
226  {
227  nummodes.push_back( nummodes[0] );
228  }
229  }
230  ASSERTL0(nummodes.size() == m_dim,"Number of modes should match mesh dimension");
231 
232 
233  const char *fStr = expansion->Attribute("FIELDS");
234  if(fStr)
235  {
236  std::string fieldStr = fStr;
237  bool valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldName);
238  ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
239 
240  for (int i = 0; i < fieldName.size(); ++i)
241  {
242  if (m_fieldNameToId.count(fieldName[i]) == 0)
243  {
244  int k = m_fieldNameToId.size();
245  m_fieldNameToId[ fieldName[i] ] = k;
246  m_numFields++;
247  }
248  }
249  }
250  else
251  {
252  fieldName.push_back("DefaultVar");
253  int k = m_fieldNameToId.size();
254 
255  if (m_fieldNameToId.count("DefaultVar") == 0)
256  {
257  ASSERTL0(k == 0,
258  "Omitting field variables and explicitly listing " \
259  "them in different ExpansionTypes is wrong practise");
260 
261  m_fieldNameToId[ "DefaultVar" ] = k;
262  m_numFields++;
263  }
264  }
265 
266  std::string compositeStr = expansion->Attribute("COMPOSITE");
267  ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
268  int beg = compositeStr.find_first_of("[");
269  int end = compositeStr.find_first_of("]");
270  std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
271  bool parseGood = ParseUtils::GenerateSeqVector(compositeListStr.c_str(), composite);
272  ASSERTL0(parseGood && !composite.empty(),
273  (std::string("Unable to read composite index range: ") + compositeListStr).c_str());
274 
275 
276  // construct mapping (field name, CompositeID) -> nummodes
277  for (int i = 0; i < composite.size(); ++i)
278  {
279  for (int j = 0; j < fieldName.size(); j++)
280  {
281  m_expansions[composite[i]][fieldName[j]] = nummodes;
282  }
283  }
284 
285  expansion = expansion->NextSiblingElement("E");
286  }
287  }
288 
289 
290 
291 
292  void MeshPartition::ReadGeometry(const LibUtilities::SessionReaderSharedPtr& pSession)
293  {
294  TiXmlElement* x;
295  TiXmlElement *vGeometry, *vSubElement;
296 
297  vGeometry = pSession->GetElement("Nektar/Geometry");
298  m_dim = atoi(vGeometry->Attribute("DIM"));
299 
300  // Read mesh vertices
301  vSubElement = pSession->GetElement("Nektar/Geometry/Vertex");
302 
303  // Retrieve any VERTEX attributes specifying mesh transforms
304  std::string attr[] = {"XSCALE", "YSCALE", "ZSCALE",
305  "XMOVE", "YMOVE", "ZMOVE" };
306  for (int i = 0; i < 6; ++i)
307  {
308  const char *val = vSubElement->Attribute(attr[i].c_str());
309  if (val)
310  {
311  m_vertexAttributes[attr[i]] = std::string(val);
312  }
313  }
314 
315  // check to see if compressed
316  std::string IsCompressed;
317  vSubElement->QueryStringAttribute("COMPRESSED",&IsCompressed);
318 
319  if(IsCompressed.size())
320  {
321  ASSERTL0(boost::iequals(IsCompressed,
322  CompressData::GetCompressString()),
323  "Compressed formats do not match. Expected :"
324  + CompressData::GetCompressString()
325  + "but got "+ std::string(IsCompressed));
326 
327  m_isCompressed = true;
328 
329  // Extract the vertex body
330  TiXmlNode* vertexChild = vSubElement->FirstChild();
331  ASSERTL0(vertexChild, "Unable to extract the data "
332  "from the compressed vertex tag.");
333 
334  std::string vertexStr;
335  if (vertexChild->Type() == TiXmlNode::TINYXML_TEXT)
336  {
337  vertexStr += vertexChild->ToText()->ValueStr();
338  }
339 
340  std::vector<MeshVertex> vertData;
341  CompressData::ZlibDecodeFromBase64Str(vertexStr,vertData);
342 
343  for(int i = 0; i < vertData.size(); ++i)
344  {
345  m_meshVertices[vertData[i].id] = vertData[i];
346  }
347  }
348  else
349  {
350  x = vSubElement->FirstChildElement();
351 
352  while(x)
353  {
354  TiXmlAttribute* y = x->FirstAttribute();
355  ASSERTL0(y, "Failed to get attribute.");
356  MeshVertex v;
357  v.id = y->IntValue();
358  std::vector<std::string> vCoords;
359  std::string vCoordStr = x->FirstChild()->ToText()->Value();
360  boost::split(vCoords, vCoordStr, boost::is_any_of("\t "));
361  v.x = atof(vCoords[0].c_str());
362  v.y = atof(vCoords[1].c_str());
363  v.z = atof(vCoords[2].c_str());
364  m_meshVertices[v.id] = v;
365  x = x->NextSiblingElement();
366  }
367  }
368 
369  // Read mesh edges
370  if (m_dim >= 2)
371  {
372  vSubElement = pSession->GetElement("Nektar/Geometry/Edge");
373  ASSERTL0(vSubElement, "Cannot read edges");
374 
375  // check to see if compressed
376  std::string IsCompressed;
377  vSubElement->QueryStringAttribute("COMPRESSED",&IsCompressed);
378 
379  if(IsCompressed.size())
380  {
381  ASSERTL0(boost::iequals(IsCompressed,
382  CompressData::GetCompressString()),
383  "Compressed formats do not match. Expected :"
384  + CompressData::GetCompressString()
385  + " but got "
386  + boost::lexical_cast<std::string>(IsCompressed));
387 
388  m_isCompressed = true;
389 
390  // Extract the edge body
391  TiXmlNode* edgeChild = vSubElement->FirstChild();
392  ASSERTL0(edgeChild,
393  "Unable to extract the data from the compressed "
394  "edge tag.");
395 
396  std::string edgeStr;
397 
398  if (edgeChild->Type() == TiXmlNode::TINYXML_TEXT)
399  {
400  edgeStr += edgeChild->ToText()->ValueStr();
401  }
402 
403  std::vector<MeshEdge> edgeData;
404  CompressData::ZlibDecodeFromBase64Str(edgeStr,edgeData);
405 
406  for(int i = 0; i < edgeData.size(); ++i)
407  {
408  MeshEntity e;
409  e.id = edgeData[i].id;
410  e.list.push_back(edgeData[i].v0);
411  e.list.push_back(edgeData[i].v1);
412  m_meshEdges[e.id] = e;
413  }
414  }
415  else
416  {
417  x = vSubElement->FirstChildElement();
418 
419  while(x)
420  {
421  TiXmlAttribute* y = x->FirstAttribute();
422  ASSERTL0(y, "Failed to get attribute.");
423  MeshEntity e;
424  e.id = y->IntValue();
425  e.type = 'E';
426  std::vector<std::string> vVertices;
427  std::string vVerticesString = x->FirstChild()->ToText()->Value();
428  boost::split(vVertices, vVerticesString, boost::is_any_of("\t "));
429  e.list.push_back(atoi(vVertices[0].c_str()));
430  e.list.push_back(atoi(vVertices[1].c_str()));
431  m_meshEdges[e.id] = e;
432  x = x->NextSiblingElement();
433  }
434  }
435  }
436 
437  // Read mesh faces
438  if (m_dim == 3)
439  {
440  vSubElement = pSession->GetElement("Nektar/Geometry/Face");
441  ASSERTL0(vSubElement, "Cannot read faces.");
442  x = vSubElement->FirstChildElement();
443 
444  while(x)
445  {
446  // check to see if compressed
447  std::string IsCompressed;
448  x->QueryStringAttribute("COMPRESSED",&IsCompressed);
449 
450  if(IsCompressed.size())
451  {
452  ASSERTL0(boost::iequals(IsCompressed,
453  CompressData::GetCompressString()),
454  "Compressed formats do not match. Expected :"
455  + CompressData::GetCompressString()
456  + " but got "
457  + boost::lexical_cast<std::string>(
458  IsCompressed));
459  m_isCompressed = true;
460 
461  // Extract the edge body
462  TiXmlNode* faceChild = x->FirstChild();
463  ASSERTL0(faceChild, "Unable to extract the data from "
464  "the compressed edge tag.");
465 
466  std::string faceStr;
467  if (faceChild->Type() == TiXmlNode::TINYXML_TEXT)
468  {
469  faceStr += faceChild->ToText()->ValueStr();
470  }
471 
472  // uncompress and fill in values.
473  const std::string val= x->Value();
474 
475  if(boost::iequals(val,"T")) // triangle
476  {
477  std::vector<MeshTri> faceData;
478  CompressData::ZlibDecodeFromBase64Str(faceStr,
479  faceData);
480 
481  for(int i= 0; i < faceData.size(); ++i)
482  {
483  MeshEntity f;
484 
485  f.id = faceData[i].id;
486  f.type = 'T';
487  for (int j = 0; j < 3; ++j)
488  {
489  f.list.push_back(faceData[i].e[j]);
490  }
491  m_meshFaces[f.id] = f;
492  }
493  }
494  else if (boost::iequals(val,"Q"))
495  {
496  std::vector<MeshQuad> faceData;
497  CompressData::ZlibDecodeFromBase64Str(faceStr,
498  faceData);
499 
500  for(int i= 0; i < faceData.size(); ++i)
501  {
502  MeshEntity f;
503 
504  f.id = faceData[i].id;
505  f.type = 'Q';
506  for (int j = 0; j < 4; ++j)
507  {
508  f.list.push_back(faceData[i].e[j]);
509  }
510  m_meshFaces[f.id] = f;
511  }
512  }
513  else
514  {
515  ASSERTL0(false,"Unrecognised face tag");
516  }
517  }
518  else
519  {
520  TiXmlAttribute* y = x->FirstAttribute();
521  ASSERTL0(y, "Failed to get attribute.");
522  MeshEntity f;
523  f.id = y->IntValue();
524  f.type = x->Value()[0];
525  std::vector<std::string> vEdges;
526  std::string vEdgeStr = x->FirstChild()->ToText()->Value();
527  boost::split(vEdges, vEdgeStr, boost::is_any_of("\t "));
528  for (int i = 0; i < vEdges.size(); ++i)
529  {
530  f.list.push_back(atoi(vEdges[i].c_str()));
531  }
532  m_meshFaces[f.id] = f;
533  }
534  x = x->NextSiblingElement();
535  }
536  }
537 
538  // Read mesh elements
539  vSubElement = pSession->GetElement("Nektar/Geometry/Element");
540  ASSERTL0(vSubElement, "Cannot read elements.");
541  x = vSubElement->FirstChildElement();
542  while(x)
543  {
544  // check to see if compressed
545  std::string IsCompressed;
546  x->QueryStringAttribute("COMPRESSED",&IsCompressed);
547 
548  if(IsCompressed.size())
549  {
550  ASSERTL0(boost::iequals(IsCompressed,
551  CompressData::GetCompressString()),
552  "Compressed formats do not match. Expected :"
553  + CompressData::GetCompressString()
554  + " but got "
555  + boost::lexical_cast<std::string>(IsCompressed));
556  m_isCompressed = true;
557 
558  // Extract the body
559  TiXmlNode* child = x->FirstChild();
560  ASSERTL0(child, "Unable to extract the data from the "
561  "compressed element tag.");
562 
563  std::string elmtStr;
564  if (child->Type() == TiXmlNode::TINYXML_TEXT)
565  {
566  elmtStr += child->ToText()->ValueStr();
567  }
568 
569  // uncompress and fill in values.
570  const std::string val= x->Value();
571 
572  switch(val[0])
573  {
574  case 'S': // segment
575  {
576  std::vector<MeshEdge> data;
577  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
578 
579  for(int i= 0; i < data.size(); ++i)
580  {
581  MeshEntity e;
582  e.id = data[i].id;
583  e.type = 'S';
584  e.list.push_back(data[i].v0);
585  e.list.push_back(data[i].v1);
586  m_meshElements[e.id] = e;
587  }
588  }
589  break;
590  case 'T': // triangle
591  {
592  std::vector<MeshTri> data;
593  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
594 
595  for(int i= 0; i < data.size(); ++i)
596  {
597  MeshEntity f;
598  f.id = data[i].id;
599  f.type = 'T';
600  for (int j = 0; j < 3; ++j)
601  {
602  f.list.push_back(data[i].e[j]);
603  }
604  m_meshElements[f.id] = f;
605  }
606  }
607  break;
608  case 'Q': // quadrilateral
609  {
610  std::vector<MeshQuad> data;
611  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
612 
613  for(int i= 0; i < data.size(); ++i)
614  {
615  MeshEntity f;
616  f.id = data[i].id;
617  f.type = 'Q';
618  for (int j = 0; j < 4; ++j)
619  {
620  f.list.push_back(data[i].e[j]);
621  }
622  m_meshElements[f.id] = f;
623  }
624  }
625  break;
626  case 'A': // tetrahedron
627  {
628  std::vector<MeshTet> data;
629  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
630 
631  for(int i= 0; i < data.size(); ++i)
632  {
633  MeshEntity f;
634  f.id = data[i].id;
635  f.type = 'A';
636  for (int j = 0; j < 4; ++j)
637  {
638  f.list.push_back(data[i].f[j]);
639  }
640  m_meshElements[f.id] = f;
641  }
642  }
643  break;
644  case 'P': // prism
645  {
646  std::vector<MeshPyr> data;
647  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
648 
649  for(int i= 0; i < data.size(); ++i)
650  {
651  MeshEntity f;
652  f.id = data[i].id;
653  f.type = 'P';
654  for (int j = 0; j < 5; ++j)
655  {
656  f.list.push_back(data[i].f[j]);
657  }
658  m_meshElements[f.id] = f;
659  }
660  }
661  break;
662  case 'R': // pyramid
663  {
664  std::vector<MeshPrism> data;
665  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
666 
667  for(int i= 0; i < data.size(); ++i)
668  {
669  MeshEntity f;
670  f.id = data[i].id;
671  f.type = 'R';
672  for (int j = 0; j < 5; ++j)
673  {
674  f.list.push_back(data[i].f[j]);
675  }
676  m_meshElements[f.id] = f;
677  }
678  }
679  break;
680  case 'H': // hexahedron
681  {
682  std::vector<MeshHex> data;
683  CompressData::ZlibDecodeFromBase64Str(elmtStr,data);
684 
685  for(int i= 0; i < data.size(); ++i)
686  {
687  MeshEntity f;
688  f.id = data[i].id;
689  f.type = 'H';
690  for (int j = 0; j < 6; ++j)
691  {
692  f.list.push_back(data[i].f[j]);
693  }
694  m_meshElements[f.id] = f;
695  }
696  }
697  break;
698  default:
699  ASSERTL0(false,"Unrecognised element tag");
700  }
701  }
702  else
703  {
704  TiXmlAttribute* y = x->FirstAttribute();
705  ASSERTL0(y, "Failed to get attribute.");
706  MeshEntity e;
707  e.id = y->IntValue();
708  std::vector<std::string> vItems;
709  std::string vItemStr = x->FirstChild()->ToText()->Value();
710  boost::split(vItems, vItemStr, boost::is_any_of("\t "));
711  for (int i = 0; i < vItems.size(); ++i)
712  {
713  e.list.push_back(atoi(vItems[i].c_str()));
714  }
715  e.type = x->Value()[0];
716  m_meshElements[e.id] = e;
717  }
718  x = x->NextSiblingElement();
719  }
720 
721  // Read mesh curves
722  if (pSession->DefinesElement("Nektar/Geometry/Curved"))
723  {
724  vSubElement = pSession->GetElement("Nektar/Geometry/Curved");
725 
726  // check to see if compressed
727  std::string IsCompressed;
728  vSubElement->QueryStringAttribute("COMPRESSED",&IsCompressed);
729 
730  x = vSubElement->FirstChildElement();
731  while(x)
732  {
733  if(IsCompressed.size())
734  {
735  ASSERTL0(boost::iequals(IsCompressed,
736  CompressData::GetCompressString()),
737  "Compressed formats do not match. Expected :"
738  + CompressData::GetCompressString()
739  + " but got "
740  + boost::lexical_cast<std::string>(
741  IsCompressed));
742 
743  m_isCompressed = true;
744 
745  const char *entitytype = x->Value();
746  // The compressed curved information is stored
747  // in two parts: MeshCurvedInfo and
748  // MeshCurvedPts. MeshCurvedPts is just a
749  // list of MeshVertex values of unique vertex
750  // values from which we can make edges and
751  // faces.
752  //
753  // Then there is a list of NekInt64 pieces of
754  // information which make a MeshCurvedInfo
755  // struct. This contains information such as
756  // the curve id, the entity id, the number of
757  // curved points and the offset of where these
758  // points are stored in the pts vector of
759  // MeshVertex values. Finally the point type
760  // is also stored but in NekInt64 format
761  // rather than an enum for binary stride
762  // compatibility.
763  if(boost::iequals(entitytype,"E")||
764  boost::iequals(entitytype,"F"))
765  {
766  // read in data
767  std::string elmtStr;
768  TiXmlNode* child = x->FirstChild();
769 
770  if (child->Type() == TiXmlNode::TINYXML_TEXT)
771  {
772  elmtStr = child->ToText()->ValueStr();
773  }
774 
775  std::vector<MeshCurvedInfo> cinfo;
776  CompressData::ZlibDecodeFromBase64Str(elmtStr,
777  cinfo);
778 
779  // unpack list of curved edge or faces
780  for(int i = 0; i < cinfo.size(); ++i)
781  {
782  MeshCurved c;
783  c.id = cinfo[i].id;
784  c.entitytype = entitytype[0];
785  c.entityid = cinfo[i].entityid;
786  c.npoints = cinfo[i].npoints;
787  c.type = kPointsTypeStr[cinfo[i].ptype];
788  c.ptid = cinfo[i].ptid;
789  c.ptoffset = cinfo[i].ptoffset;
790  m_meshCurved[std::make_pair(c.entitytype,
791  c.id)] = c;
792  }
793  }
794  else if(boost::iequals(entitytype,"DATAPOINTS"))
795  {
796  MeshCurvedPts cpts;
797  NekInt id;
798 
799  ASSERTL0(x->Attribute("ID", &id),
800  "Failed to get ID from PTS section");
801  cpts.id = id;
802 
803  // read in data
804  std::string elmtStr;
805 
806  TiXmlElement* DataIdx =
807  x->FirstChildElement("INDEX");
808  ASSERTL0(DataIdx,
809  "Cannot read data index tag in compressed "
810  "curved section");
811 
812  TiXmlNode* child = DataIdx->FirstChild();
813  if (child->Type() == TiXmlNode::TINYXML_TEXT)
814  {
815  elmtStr = child->ToText()->ValueStr();
816  }
817 
818  CompressData::ZlibDecodeFromBase64Str(elmtStr,
819  cpts.index);
820 
821  TiXmlElement* DataPts
822  = x->FirstChildElement("POINTS");
823  ASSERTL0(DataPts,
824  "Cannot read data pts tag in compressed "
825  "curved section");
826 
827  child = DataPts->FirstChild();
828  if (child->Type() == TiXmlNode::TINYXML_TEXT)
829  {
830  elmtStr = child->ToText()->ValueStr();
831  }
832 
833  CompressData::ZlibDecodeFromBase64Str(elmtStr,
834  cpts.pts);
835 
836  m_meshCurvedPts[cpts.id] = cpts;
837  }
838  else
839  {
840  ASSERTL0(false, "Unknown tag in curved section");
841  }
842 
843  x = x->NextSiblingElement();
844  }
845  else
846  {
847  MeshCurved c;
848  ASSERTL0(x->Attribute("ID", &c.id),
849  "Failed to get attribute ID");
850  c.type = std::string(x->Attribute("TYPE"));
851  ASSERTL0(!c.type.empty(),
852  "Failed to get attribute TYPE");
853  ASSERTL0(x->Attribute("NUMPOINTS", &c.npoints),
854  "Failed to get attribute NUMPOINTS");
855  c.data = x->FirstChild()->ToText()->Value();
856  c.entitytype = x->Value()[0];
857 
858  if (c.entitytype == "E")
859  {
860  ASSERTL0(x->Attribute("EDGEID", &c.entityid),
861  "Failed to get attribute EDGEID");
862  }
863  else if (c.entitytype == "F")
864  {
865  ASSERTL0(x->Attribute("FACEID", &c.entityid),
866  "Failed to get attribute FACEID");
867  }
868  else
869  {
870  ASSERTL0(false, "Unknown curve type.");
871  }
872 
873  m_meshCurved[std::make_pair(c.entitytype, c.id)] = c;
874  x = x->NextSiblingElement();
875  }
876  }
877  }
878 
879  // Read composites
880  vSubElement = pSession->GetElement("Nektar/Geometry/Composite");
881  ASSERTL0(vSubElement, "Cannot read composites.");
882  x = vSubElement->FirstChildElement();
883  while(x)
884  {
885  TiXmlAttribute* y = x->FirstAttribute();
886  ASSERTL0(y, "Failed to get attribute.");
887  MeshEntity c;
888  c.id = y->IntValue();
889  std::string vSeqStr = x->FirstChild()->ToText()->Value();
890  c.type = vSeqStr[0];
891  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
892  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
893  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
894 
895  std::vector<unsigned int> vSeq;
896  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
897 
898  for (int i = 0; i < vSeq.size(); ++i)
899  {
900  c.list.push_back(vSeq[i]);
901  }
902  m_meshComposites[c.id] = c;
903  x = x->NextSiblingElement();
904  }
905 
906  // Read Domain
907  vSubElement = pSession->GetElement("Nektar/Geometry/Domain");
908  ASSERTL0(vSubElement, "Cannot read domain");
909  std::string vSeqStr = vSubElement->FirstChild()->ToText()->Value();
910  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
911  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
912  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
913  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), m_domain);
914  }
915 
916  void MeshPartition::PrintPartInfo(std::ostream &out)
917  {
918  int nElmt = boost::num_vertices(m_mesh);
919  int nPart = m_localPartition.size();
920 
921  out << "# Partition information:" << std::endl;
922  out << "# No. elements : " << nElmt << std::endl;
923  out << "# No. partitions: " << nPart << std::endl;
924  out << "# ID nElmt nLocDof nBndDof" << std::endl;
925 
926  BoostVertexIterator vertit, vertit_end;
927  std::vector<int> partElmtCount(nPart, 0);
928  std::vector<int> partLocCount (nPart, 0);
929  std::vector<int> partBndCount (nPart, 0);
930 
931  std::map<int, int> elmtSizes;
932  std::map<int, int> elmtBndSizes;
933 
934  for (unsigned int i = 0; i < m_domain.size(); ++i)
935  {
936  int cId = m_domain[i];
937  NummodesPerField npf = m_expansions[cId];
938 
939  for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
940  {
941  ASSERTL0(it->second.size() == m_dim,
942  " Number of directional" \
943  " modes in expansion spec for composite id = " +
944  boost::lexical_cast<std::string>(cId) +
945  " and field " +
946  boost::lexical_cast<std::string>(it->first) +
947  " does not correspond to mesh dimension");
948 
949  int na = it->second[0];
950  int nb = it->second[1];
951  int nc = 0;
952  if (m_dim == 3)
953  {
954  nc = it->second[2];
955  }
956 
957  int weight = CalculateElementWeight(
958  m_meshComposites[cId].type, false, na, nb, nc);
959  int bndWeight = CalculateElementWeight(
960  m_meshComposites[cId].type, true, na, nb, nc);
961 
962  for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
963  {
964  int elid = m_meshComposites[cId].list[j];
965  elmtSizes[elid] = weight;
966  elmtBndSizes[elid] = bndWeight;
967  }
968  }
969  }
970 
971  for (boost::tie(vertit, vertit_end) = boost::vertices(m_mesh);
972  vertit != vertit_end; ++vertit)
973  {
974  int partId = m_mesh[*vertit].partition;
975  partElmtCount[partId]++;
976  partLocCount [partId] += elmtSizes[m_mesh[*vertit].id];
977  partBndCount [partId] += elmtBndSizes[m_mesh[*vertit].id];
978  }
979 
980  for (int i = 0; i < nPart; ++i)
981  {
982  out << i << " " << partElmtCount[i] << " " << partLocCount[i] << " " << partBndCount[i] << std::endl;
983  }
984  }
985 
986  void MeshPartition::ReadConditions(const SessionReaderSharedPtr& pSession)
987  {
988  if (!pSession->DefinesElement("Nektar/Conditions/SolverInfo"))
989  {
990  // No SolverInfo = no change of default action to weight
991  // mesh graph.
992  return;
993  }
994 
995  TiXmlElement* solverInfoElement =
996  pSession->GetElement("Nektar/Conditions/SolverInfo");
997 
998  TiXmlElement* solverInfo =
999  solverInfoElement->FirstChildElement("I");
1000  ASSERTL0(solverInfo, "Cannot read SolverInfo tags");
1001 
1002  while (solverInfo)
1003  {
1004  // read the property name
1005  ASSERTL0(solverInfo->Attribute("PROPERTY"),
1006  "Missing PROPERTY attribute in solver info "
1007  "section. ");
1008  std::string solverProperty =
1009  solverInfo->Attribute("PROPERTY");
1010  ASSERTL0(!solverProperty.empty(),
1011  "Solver info properties must have a non-empty "
1012  "name. ");
1013  // make sure that solver property is capitalised
1014  std::string solverPropertyUpper =
1015  boost::to_upper_copy(solverProperty);
1016 
1017 
1018  // read the value
1019  ASSERTL0(solverInfo->Attribute("VALUE"),
1020  "Missing VALUE attribute in solver info section. ");
1021  std::string solverValue = solverInfo->Attribute("VALUE");
1022  ASSERTL0(!solverValue.empty(),
1023  "Solver info properties must have a non-empty value");
1024  // make sure that property value is capitalised
1025  std::string propertyValueUpper =
1026  boost::to_upper_copy(solverValue);
1027 
1028  if (solverPropertyUpper == "WEIGHTPARTITIONS")
1029  {
1030  if (propertyValueUpper != "UNIFORM")
1031  {
1032  m_weightingRequired = true;
1033  }
1034  return;
1035  }
1036  solverInfo = solverInfo->NextSiblingElement("I");
1037  }
1038  }
1039 
1040 
1041  /*
1042  * Calculate element weights based on
1043  * - element type (Q,T,H,P,R,A)
1044  * - nummodes in expansion which this element belongs to via composite.
1045  *
1046  * For each element we prepare two vertex weightings, one associated
1047  * with the number of matrix elements associated with it (to balance
1048  * matrix multiplication work) and another associated
1049  * with all work which scales linearly with the number of its
1050  * coefficients: communication, vector updates etc.
1051  *
1052  * \todo Refactor this code to explicitly represent performance model
1053  * and flexibly generate graph vertex weights depending on perf data.
1054  */
1055  void MeshPartition::WeightElements()
1056  {
1057  std::vector<unsigned int> weight(m_numFields, 1);
1058  std::map<int, MeshEntity>::iterator eIt;
1059  for (eIt = m_meshElements.begin(); eIt != m_meshElements.end(); ++eIt)
1060  {
1061  m_vertWeights[eIt->first] = weight;
1062  }
1063 
1064  for (unsigned int i = 0; i < m_domain.size(); ++i)
1065  {
1066  int cId = m_domain[i];
1067  NummodesPerField npf = m_expansions[cId];
1068 
1069  for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
1070  {
1071  ASSERTL0(it->second.size() == m_dim,
1072  " Number of directional" \
1073  " modes in expansion spec for composite id = " +
1074  boost::lexical_cast<std::string>(cId) +
1075  " and field " +
1076  boost::lexical_cast<std::string>(it->first) +
1077  " does not correspond to mesh dimension");
1078 
1079  int na = it->second[0];
1080  int nb = 0;
1081  int nc = 0;
1082  if (m_dim >= 2)
1083  {
1084  nb = it->second[1];
1085  }
1086  if (m_dim == 3)
1087  {
1088  nc = it->second[2];
1089  }
1090 
1091  int bndWeight = CalculateElementWeight(
1092  m_meshComposites[cId].type, true, na, nb, nc);
1093 
1094  for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
1095  {
1096  int elmtId = m_meshComposites[cId].list[j];
1097  m_vertWeights[elmtId][m_fieldNameToId[it->first]] = bndWeight;
1098  }
1099  }
1100  } // for i
1101  }
1102 
1103  void MeshPartition::CreateGraph(BoostSubGraph& pGraph)
1104  {
1105  // Maps edge/face to first mesh element id.
1106  // On locating second mesh element id, graph edge is created instead.
1107  std::map<int, int> vGraphEdges;
1108  std::map<int, MeshEntity>::iterator eIt;
1109  int vcnt = 0;
1110 
1111  for (eIt = m_meshElements.begin(); eIt != m_meshElements.end();
1112  ++eIt, ++vcnt)
1113  {
1114  BoostVertex v = boost::add_vertex(pGraph);
1115  pGraph[v].id = eIt->first;
1116  pGraph[v].partition = 0;
1117 
1118  if (m_weightingRequired)
1119  {
1120  pGraph[v].weight = m_vertWeights[eIt->first];
1121  }
1122 
1123  // Process element entries and add graph edges
1124  for (unsigned j = 0; j < eIt->second.list.size(); ++j)
1125  {
1126  int eId = eIt->second.list[j];
1127 
1128  // Look to see if we've examined this edge/face before
1129  // If so, we've got both graph vertices so add edge
1130  if (vGraphEdges.find(eId) != vGraphEdges.end())
1131  {
1132  BoostEdge e = boost::add_edge(vcnt, vGraphEdges[eId], pGraph).first;
1133  pGraph[e].id = vcnt;
1134  }
1135  else
1136  {
1137  vGraphEdges[eId] = vcnt;
1138  }
1139  }
1140  }
1141  }
1142 
1143  /**
1144  * @brief Partition the graph.
1145  *
1146  * This routine partitions the graph @p pGraph into @p nParts, producing
1147  * subgraphs that are populated in @p pLocalPartition. If the @p
1148  * overlapping option is set (which is used for post-processing
1149  * purposes), the resulting partitions are extended to cover
1150  * neighbouring elements by additional vertex on the dual graph, which
1151  * produces overlapping partitions (i.e. the intersection of two
1152  * connected partitions is non-empty).
1153  *
1154  * @param pGraph Graph to be partitioned.
1155  * @param nParts Number of partitions.
1156  * @param pLocalPartition Vector of sub-graphs representing each
1157  * partition.
1158  * @param overlapping True if resulting partitions should overlap.
1159  */
1160  void MeshPartition::PartitionGraph(BoostSubGraph& pGraph,
1161  int nParts,
1162  std::vector<BoostSubGraph>& pLocalPartition,
1163  bool overlapping)
1164  {
1165  int i;
1166  int nGraphVerts = boost::num_vertices(pGraph);
1167  int nGraphEdges = boost::num_edges(pGraph);
1168 
1169  // Convert boost graph into CSR format
1170  BoostVertexIterator vertit, vertit_end;
1171  BoostAdjacencyIterator adjvertit, adjvertit_end;
1172  Array<OneD, int> part(nGraphVerts,0);
1173 
1174  if (m_comm->GetRowComm()->TreatAsRankZero())
1175  {
1176  int acnt = 0;
1177  int vcnt = 0;
1178  int nWeight = nGraphVerts;
1179  Array<OneD, int> xadj(nGraphVerts+1,0);
1180  Array<OneD, int> adjncy(2*nGraphEdges);
1181  Array<OneD, int> vwgt(nWeight, 1);
1182  Array<OneD, int> vsize(nGraphVerts, 1);
1183 
1184  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
1185  vertit != vertit_end;
1186  ++vertit)
1187  {
1188  for ( boost::tie(adjvertit, adjvertit_end) = boost::adjacent_vertices(*vertit,pGraph);
1189  adjvertit != adjvertit_end;
1190  ++adjvertit)
1191  {
1192  adjncy[acnt++] = *adjvertit;
1193  }
1194 
1195  xadj[++vcnt] = acnt;
1196 
1197  if (m_weightingRequired)
1198  {
1199  vwgt[vcnt-1] = pGraph[*vertit].weight[0];
1200  }
1201  else
1202  {
1203  vwgt[vcnt-1] = 1;
1204  }
1205  }
1206 
1207  // Call Metis and partition graph
1208  int vol = 0;
1209 
1210  try
1211  {
1212  //////////////////////////////////////////////////////
1213  // On a cartesian communicator do mesh partiotion just on the first column
1214  // so there is no doubt the partitions are all the same in all the columns
1215  if(m_comm->GetColumnComm()->GetRank() == 0)
1216  {
1217  // Attempt partitioning using METIS.
1218  int ncon = 1;
1219  PartitionGraphImpl(nGraphVerts, ncon, xadj, adjncy, vwgt, vsize, nParts, vol, part);
1220 
1221  // Check METIS produced a valid partition and fix if not.
1222  CheckPartitions(nParts, part);
1223  if (!m_shared)
1224  {
1225  // distribute among columns
1226  for (i = 1; i < m_comm->GetColumnComm()->GetSize(); ++i)
1227  {
1228  m_comm->GetColumnComm()->Send(i, part);
1229  }
1230  }
1231  }
1232  else
1233  {
1234  m_comm->GetColumnComm()->Recv(0, part);
1235  }
1236  if (!m_shared)
1237  {
1238  m_comm->GetColumnComm()->Block();
1239 
1240  //////////////////////////////////
1241  // distribute among rows
1242  for (i = 1; i < m_comm->GetRowComm()->GetSize(); ++i)
1243  {
1244  m_comm->GetRowComm()->Send(i, part);
1245  }
1246  }
1247  }
1248  catch (...)
1249  {
1250  NEKERROR(ErrorUtil::efatal,
1251  "Error in calling metis to partition graph.");
1252  }
1253  }
1254  else
1255  {
1256  m_comm->GetRowComm()->Recv(0, part);
1257  }
1258 
1259  // Create boost subgraph for this process's partitions
1260  int nCols = nParts;
1261  pLocalPartition.resize(nCols);
1262  for (i = 0; i < nCols; ++i)
1263  {
1264  pLocalPartition[i] = pGraph.create_subgraph();
1265  }
1266 
1267  // Populate subgraph
1268  i = 0;
1269  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
1270  vertit != vertit_end;
1271  ++vertit, ++i)
1272  {
1273  pGraph[*vertit].partition = part[i];
1274  boost::add_vertex(i, pLocalPartition[part[i]]);
1275  }
1276 
1277  // If the overlapping option is set (for post-processing purposes),
1278  // add vertices that correspond to the neighbouring elements.
1279  if (overlapping)
1280  {
1281  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
1282  vertit != vertit_end;
1283  ++vertit)
1284  {
1285  for (boost::tie(adjvertit, adjvertit_end) = boost::adjacent_vertices(*vertit,pGraph);
1286  adjvertit != adjvertit_end; ++adjvertit)
1287  {
1288  if(part[*adjvertit] != part[*vertit])
1289  {
1290  boost::add_vertex(*adjvertit, pLocalPartition[part[*vertit]]);
1291  }
1292  }
1293  }
1294  }
1295  }
1296 
1297 
1298  void MeshPartition::CheckPartitions(int nParts, Array<OneD, int> &pPart)
1299  {
1300  unsigned int i = 0;
1301  unsigned int cnt = 0;
1302  bool valid = true;
1303 
1304  // Check that every process has at least one element assigned
1305  for (i = 0; i < nParts; ++i)
1306  {
1307  cnt = std::count(pPart.begin(), pPart.end(), i);
1308  if (cnt == 0)
1309  {
1310  valid = false;
1311  }
1312  }
1313 
1314  // If METIS produced an invalid partition, repartition naively.
1315  // Elements are assigned to processes in a round-robin fashion.
1316  // It is assumed that METIS failure only occurs when the number of
1317  // elements is approx. the number of processes, so this approach
1318  // should not be too inefficient communication-wise.
1319  if (!valid)
1320  {
1321  for (i = 0; i < pPart.num_elements(); ++i)
1322  {
1323  pPart[i] = i % nParts;
1324  }
1325  }
1326  }
1327 
1328 
1329  void MeshPartition::OutputPartition(
1330  LibUtilities::SessionReaderSharedPtr& pSession,
1331  BoostSubGraph& pGraph,
1332  TiXmlElement* pNektar)
1333  {
1334  // Write Geometry data
1335  std::string vDim = pSession->GetElement("Nektar/Geometry")->Attribute("DIM");
1336  std::string vSpace = pSession->GetElement("Nektar/Geometry")->Attribute("SPACE");
1337  std::string vPart = boost::lexical_cast<std::string>(pGraph[*boost::vertices(pGraph).first].partition);
1338  TiXmlElement* vElmtGeometry = new TiXmlElement("GEOMETRY");
1339  vElmtGeometry->SetAttribute("DIM", vDim);
1340  vElmtGeometry->SetAttribute("SPACE", vSpace);
1341  vElmtGeometry->SetAttribute("PARTITION", vPart);
1342 
1343  TiXmlElement *vVertex = new TiXmlElement("VERTEX");
1344  TiXmlElement *vEdge = new TiXmlElement("EDGE");
1345  TiXmlElement *vFace = new TiXmlElement("FACE");
1346  TiXmlElement *vElement = new TiXmlElement("ELEMENT");
1347  TiXmlElement *vCurved = new TiXmlElement("CURVED");
1348  TiXmlElement *vComposite = new TiXmlElement("COMPOSITE");
1349  TiXmlElement *vDomain = new TiXmlElement("DOMAIN");
1350 
1351  TiXmlElement *x;
1352  TiXmlText *y;
1353 
1354  BoostVertexIterator vertit, vertit_end;
1355  int id;
1356 
1357  std::map<int, MeshEntity> vComposites;
1358  std::map<int, MeshEntity> vElements;
1359  std::map<int, MeshEntity> vEdges;
1360  std::map<int, MeshEntity> vFaces;
1361  std::map<int, MeshVertex> vVertices;
1362  std::map<int, MeshEntity>::iterator vIt;
1363  std::map<int, MeshVertex>::iterator vVertIt;
1364  std::map<std::string, std::string>::iterator vAttrIt;
1365 
1366  std::vector<unsigned int> idxList;
1367 
1368  // Populate lists of elements, edges and vertices required.
1369  for (boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
1370  vertit != vertit_end;
1371  ++vertit)
1372  {
1373  id = pGraph[*vertit].id;
1374  vElements[id] = m_meshElements[pGraph[*vertit].id];
1375  }
1376 
1377  std::map<int, MeshEntity> * vNext = &vElements;
1378  switch (m_dim)
1379  {
1380  case 3:
1381  {
1382  // Compile list of faces
1383  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
1384  {
1385  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1386  {
1387  id = vIt->second.list[j];
1388  vFaces[id] = m_meshFaces[id];
1389  }
1390  }
1391  vNext = &vFaces;
1392  }
1393  case 2:
1394  {
1395  // Compile list of edges
1396  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
1397  {
1398  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1399  {
1400  id = vIt->second.list[j];
1401  vEdges[id] = m_meshEdges[id];
1402  }
1403  }
1404  vNext = &vEdges;
1405  }
1406  case 1:
1407  {
1408  // Compile list of vertices
1409  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
1410  {
1411  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1412  {
1413  id = vIt->second.list[j];
1414  vVertices[id] = m_meshVertices[id];
1415  }
1416  }
1417  }
1418  }
1419 
1420  // Generate XML data for these mesh entities
1421  if(m_isCompressed)
1422  {
1423  std::vector<MeshVertex> vertInfo;
1424  for (vVertIt = vVertices.begin();
1425  vVertIt != vVertices.end(); vVertIt++)
1426  {
1427  MeshVertex v;
1428  v.id = vVertIt->first;
1429  v.x = vVertIt->second.x;
1430  v.y = vVertIt->second.y;
1431  v.z = vVertIt->second.z;
1432  vertInfo.push_back(v);
1433  }
1434  std::string vertStr;
1435  CompressData::ZlibEncodeToBase64Str(vertInfo,vertStr);
1436  vVertex->SetAttribute("COMPRESSED",
1437  CompressData::GetCompressString());
1438  vVertex->SetAttribute("BITSIZE",
1439  CompressData::GetBitSizeStr());
1440  vVertex->LinkEndChild(new TiXmlText(vertStr));
1441  }
1442  else
1443  {
1444  for (vVertIt = vVertices.begin(); vVertIt != vVertices.end(); vVertIt++)
1445  {
1446  x = new TiXmlElement("V");
1447  x->SetAttribute("ID", vVertIt->first);
1448  std::stringstream vCoords;
1449  vCoords.precision(12);
1450  vCoords << std::setw(15) << vVertIt->second.x << " "
1451  << std::setw(15) << vVertIt->second.y << " "
1452  << std::setw(15) << vVertIt->second.z << " ";
1453  y = new TiXmlText(vCoords.str());
1454  x->LinkEndChild(y);
1455  vVertex->LinkEndChild(x);
1456  }
1457  }
1458 
1459  // Apply transformation attributes to VERTEX section
1460  for (vAttrIt = m_vertexAttributes.begin();
1461  vAttrIt != m_vertexAttributes.end();
1462  ++ vAttrIt)
1463  {
1464  vVertex->SetAttribute(vAttrIt->first, vAttrIt->second);
1465  }
1466 
1467  if (m_dim >= 2)
1468  {
1469  if(m_isCompressed)
1470  {
1471  std::vector<MeshEdge> edgeInfo;
1472  for (vIt = vEdges.begin(); vIt != vEdges.end(); vIt++)
1473  {
1474  MeshEdge e;
1475  e.id = vIt->first;
1476  e.v0 = vIt->second.list[0];
1477  e.v1 = vIt->second.list[1];
1478  edgeInfo.push_back(e);
1479  }
1480  std::string edgeStr;
1481  CompressData::ZlibEncodeToBase64Str(edgeInfo,edgeStr);
1482  vEdge->SetAttribute("COMPRESSED",
1483  CompressData::GetCompressString());
1484  vEdge->SetAttribute("BITSIZE",
1485  CompressData::GetBitSizeStr());
1486  vEdge->LinkEndChild(new TiXmlText(edgeStr));
1487  }
1488  else
1489  {
1490  for (vIt = vEdges.begin(); vIt != vEdges.end(); vIt++)
1491  {
1492  x = new TiXmlElement("E");
1493  x->SetAttribute("ID", vIt->first);
1494  std::stringstream vVertices;
1495  vVertices << std::setw(10) << vIt->second.list[0]
1496  << std::setw(10) << vIt->second.list[1]
1497  << " ";
1498  y = new TiXmlText(vVertices.str());
1499  x->LinkEndChild(y);
1500  vEdge->LinkEndChild(x);
1501  }
1502  }
1503  }
1504 
1505  if (m_dim >= 3)
1506  {
1507 
1508  if(m_isCompressed)
1509  {
1510  std::vector<MeshTri> TriFaceInfo;
1511  std::vector<MeshQuad> QuadFaceInfo;
1512 
1513  for (vIt = vFaces.begin(); vIt != vFaces.end(); vIt++)
1514  {
1515  switch(vIt->second.list.size())
1516  {
1517  case 3:
1518  {
1519  MeshTri f;
1520  f.id = vIt->first;
1521  for(int i = 0; i < 3; ++i)
1522  {
1523  f.e[i] = vIt->second.list[i];
1524  }
1525  TriFaceInfo.push_back(f);
1526  }
1527  break;
1528  case 4:
1529  {
1530  MeshQuad f;
1531  f.id = vIt->first;
1532  for(int i = 0; i < 4; ++i)
1533  {
1534  f.e[i] = vIt->second.list[i];
1535  }
1536  QuadFaceInfo.push_back(f);
1537  }
1538  break;
1539  default:
1540  ASSERTL0(false,"Unknown face type.");
1541  }
1542  }
1543 
1544  if(TriFaceInfo.size())
1545  {
1546  std::string vType("T");
1547  x = new TiXmlElement(vType);
1548 
1549  std::string faceStr;
1550  CompressData::ZlibEncodeToBase64Str(TriFaceInfo,
1551  faceStr);
1552  x->SetAttribute("COMPRESSED",
1553  CompressData::GetCompressString());
1554  x->SetAttribute("BITSIZE",
1555  CompressData::GetBitSizeStr());
1556  x->LinkEndChild(new TiXmlText(faceStr));
1557  vFace->LinkEndChild(x);
1558  }
1559 
1560  if(QuadFaceInfo.size())
1561  {
1562  std::string vType("Q");
1563  x = new TiXmlElement(vType);
1564  std::string faceStr;
1565  CompressData::ZlibEncodeToBase64Str(QuadFaceInfo,
1566  faceStr);
1567  x->SetAttribute("COMPRESSED",
1568  CompressData::GetCompressString());
1569  x->SetAttribute("BITSIZE",
1570  CompressData::GetBitSizeStr());
1571  x->LinkEndChild(new TiXmlText(faceStr));
1572  vFace->LinkEndChild(x);
1573  }
1574  }
1575  else
1576  {
1577  for (vIt = vFaces.begin(); vIt != vFaces.end(); vIt++)
1578  {
1579  std::string vType("F");
1580  vType[0] = vIt->second.type;
1581  x = new TiXmlElement(vType);
1582  x->SetAttribute("ID", vIt->first);
1583  std::stringstream vListStr;
1584  for (unsigned int i = 0; i < vIt->second.list.size(); ++i)
1585  {
1586  vListStr << std::setw(10) << vIt->second.list[i];
1587  }
1588  vListStr << " ";
1589  y = new TiXmlText(vListStr.str());
1590  x->LinkEndChild(y);
1591  vFace->LinkEndChild(x);
1592  }
1593  }
1594  }
1595 
1596 
1597  if(m_isCompressed)
1598  {
1599  std::vector<MeshEdge> SegInfo;
1600  std::vector<MeshTri> TriInfo;
1601  std::vector<MeshQuad> QuadInfo;
1602  std::vector<MeshTet> TetInfo;
1603  std::vector<MeshPyr> PyrInfo;
1604  std::vector<MeshPrism> PrismInfo;
1605  std::vector<MeshHex> HexInfo;
1606 
1607  //gather elemements in groups.
1608  for (vIt = vElements.begin(); vIt != vElements.end(); vIt++)
1609  {
1610  switch(vIt->second.type)
1611  {
1612  case 'S':
1613  {
1614  MeshEdge e;
1615  e.id = vIt->first;
1616  e.v0 = vIt->second.list[0];
1617  e.v1 = vIt->second.list[1];
1618  SegInfo.push_back(e);
1619  }
1620  break;
1621  case 'T':
1622  {
1623  MeshTri f;
1624  f.id = vIt->first;
1625  for(int i = 0; i < 3; ++i)
1626  {
1627  f.e[i] = vIt->second.list[i];
1628  }
1629  TriInfo.push_back(f);
1630  }
1631  break;
1632  case 'Q':
1633  {
1634  MeshQuad f;
1635  f.id = vIt->first;
1636  for(int i = 0; i < 4; ++i)
1637  {
1638  f.e[i] = vIt->second.list[i];
1639  }
1640  QuadInfo.push_back(f);
1641  }
1642  break;
1643  case 'A':
1644  {
1645  MeshTet vol;
1646  vol.id = vIt->first;
1647  for(int i = 0; i < 4; ++i)
1648  {
1649  vol.f[i] = vIt->second.list[i];
1650  }
1651  TetInfo.push_back(vol);
1652  }
1653  break;
1654  case 'P':
1655  {
1656  MeshPyr vol;
1657  vol.id = vIt->first;
1658  for(int i = 0; i < 5; ++i)
1659  {
1660  vol.f[i] = vIt->second.list[i];
1661  }
1662  PyrInfo.push_back(vol);
1663  }
1664  break;
1665  case 'R':
1666  {
1667  MeshPrism vol;
1668  vol.id = vIt->first;
1669  for(int i = 0; i < 5; ++i)
1670  {
1671  vol.f[i] = vIt->second.list[i];
1672  }
1673  PrismInfo.push_back(vol);
1674  }
1675  break;
1676  case 'H':
1677  {
1678  MeshHex vol;
1679  vol.id = vIt->first;
1680  for(int i = 0; i < 6; ++i)
1681  {
1682  vol.f[i] = vIt->second.list[i];
1683  }
1684  HexInfo.push_back(vol);
1685  }
1686  break;
1687  default:
1688  ASSERTL0(false,"Unknown element type");
1689  }
1690  }
1691 
1692  if(SegInfo.size())
1693  {
1694  std::string vType("S");
1695  x = new TiXmlElement(vType);
1696 
1697  std::string segStr;
1698  CompressData::ZlibEncodeToBase64Str(SegInfo,segStr);
1699  x->SetAttribute("COMPRESSED",
1700  CompressData::GetCompressString());
1701  x->SetAttribute("BITSIZE",
1702  CompressData::GetBitSizeStr());
1703  x->LinkEndChild(new TiXmlText(segStr));
1704  vElement->LinkEndChild(x);
1705  }
1706 
1707  if(TriInfo.size())
1708  {
1709  std::string vType("T");
1710  x = new TiXmlElement(vType);
1711 
1712  std::string faceStr;
1713  CompressData::ZlibEncodeToBase64Str(TriInfo,faceStr);
1714  x->SetAttribute("COMPRESSED",
1715  CompressData::GetCompressString());
1716  x->SetAttribute("BITSIZE",
1717  CompressData::GetBitSizeStr());
1718  x->LinkEndChild(new TiXmlText(faceStr));
1719  vElement->LinkEndChild(x);
1720  }
1721 
1722  if(QuadInfo.size())
1723  {
1724  std::string vType("Q");
1725  x = new TiXmlElement(vType);
1726  std::string faceStr;
1727  CompressData::ZlibEncodeToBase64Str(QuadInfo,faceStr);
1728  x->SetAttribute("COMPRESSED",
1729  CompressData::GetCompressString());
1730  x->SetAttribute("BITSIZE",
1731  CompressData::GetBitSizeStr());
1732  x->LinkEndChild(new TiXmlText(faceStr));
1733  vElement->LinkEndChild(x);
1734  }
1735 
1736  if(TetInfo.size())
1737  {
1738  std::string vType("A");
1739  x = new TiXmlElement(vType);
1740  std::string volStr;
1741  CompressData::ZlibEncodeToBase64Str(TetInfo,volStr);
1742  x->SetAttribute("COMPRESSED",
1743  CompressData::GetCompressString());
1744  x->SetAttribute("BITSIZE",
1745  CompressData::GetBitSizeStr());
1746  x->LinkEndChild(new TiXmlText(volStr));
1747  vElement->LinkEndChild(x);
1748  }
1749 
1750  if(PyrInfo.size())
1751  {
1752  std::string vType("P");
1753  x = new TiXmlElement(vType);
1754  std::string volStr;
1755  CompressData::ZlibEncodeToBase64Str(PyrInfo,volStr);
1756  x->SetAttribute("COMPRESSED",
1757  CompressData::GetCompressString());
1758  x->SetAttribute("BITSIZE",
1759  CompressData::GetBitSizeStr());
1760  x->LinkEndChild(new TiXmlText(volStr));
1761  vElement->LinkEndChild(x);
1762  }
1763 
1764  if(PrismInfo.size())
1765  {
1766  std::string vType("R");
1767  x = new TiXmlElement(vType);
1768  std::string volStr;
1769  CompressData::ZlibEncodeToBase64Str(PrismInfo,volStr);
1770  x->SetAttribute("COMPRESSED",
1771  CompressData::GetCompressString());
1772  x->SetAttribute("BITSIZE",
1773  CompressData::GetBitSizeStr());
1774  x->LinkEndChild(new TiXmlText(volStr));
1775  vElement->LinkEndChild(x);
1776  }
1777 
1778  if(HexInfo.size())
1779  {
1780  std::string vType("H");
1781  x = new TiXmlElement(vType);
1782  std::string volStr;
1783  CompressData::ZlibEncodeToBase64Str(HexInfo,volStr);
1784  x->SetAttribute("COMPRESSED",
1785  CompressData::GetCompressString());
1786  x->SetAttribute("BITSIZE",
1787  CompressData::GetBitSizeStr());
1788  x->LinkEndChild(new TiXmlText(volStr));
1789  vElement->LinkEndChild(x);
1790  }
1791 
1792  }
1793  else
1794  {
1795  for (vIt = vElements.begin(); vIt != vElements.end(); vIt++)
1796  {
1797  std::string vType("T");
1798  vType[0] = vIt->second.type;
1799  x = new TiXmlElement(vType.c_str());
1800  x->SetAttribute("ID", vIt->first);
1801  std::stringstream vEdges;
1802  for (unsigned i = 0; i < vIt->second.list.size(); ++i)
1803  {
1804  vEdges << std::setw(10) << vIt->second.list[i];
1805  }
1806  vEdges << " ";
1807  y = new TiXmlText(vEdges.str());
1808  x->LinkEndChild(y);
1809  vElement->LinkEndChild(x);
1810  }
1811  }
1812 
1813  if (m_dim >= 2)
1814  {
1815  std::map<MeshCurvedKey, MeshCurved>::const_iterator vItCurve;
1816 
1817  if(m_isCompressed)
1818  {
1819  std::vector<MeshCurvedInfo> edgeinfo;
1820  std::vector<MeshCurvedInfo> faceinfo;
1821  MeshCurvedPts curvedpts;
1822  curvedpts.id = 0; // assume all points are going in here
1823  int ptoffset = 0;
1824  int newidx = 0;
1825  std::map<int,int> idxmap;
1826 
1827  for (vItCurve = m_meshCurved.begin();
1828  vItCurve != m_meshCurved.end();
1829  ++vItCurve)
1830  {
1831  MeshCurved c = vItCurve->second;
1832 
1833  bool IsEdge = boost::iequals(c.entitytype,"E");
1834  bool IsFace = boost::iequals(c.entitytype,"F");
1835 
1836  if((IsEdge&&vEdges.find(c.entityid) != vEdges.end())||
1837  (IsFace&&vFaces.find(c.entityid) != vFaces.end()))
1838  {
1839  MeshCurvedInfo cinfo;
1840  // add in
1841  cinfo.id = c.id;
1842  cinfo.entityid = c.entityid;
1843  cinfo.npoints = c.npoints;
1844  for(int i = 0; i < SIZE_PointsType; ++i)
1845  {
1846  if(c.type.compare(kPointsTypeStr[i]) == 0)
1847  {
1848  cinfo.ptype = (PointsType) i;
1849  break;
1850  }
1851  }
1852  cinfo.ptid = 0; // set to just one point set
1853  cinfo.ptoffset = ptoffset;
1854  ptoffset += c.npoints;
1855 
1856  if (IsEdge)
1857  {
1858  edgeinfo.push_back(cinfo);
1859  }
1860  else
1861  {
1862  faceinfo.push_back(cinfo);
1863  }
1864 
1865  // fill in points to list.
1866  for(int i =0; i < c.npoints; ++i)
1867  {
1868  // get index from full list;
1869  int idx = m_meshCurvedPts[c.ptid]
1870  .index[c.ptoffset+i];
1871 
1872  // if index is not already in curved
1873  // points add it or set index to location
1874  if(idxmap.count(idx) == 0)
1875  {
1876  idxmap[idx] = newidx;
1877  curvedpts.index.push_back(newidx);
1878  curvedpts.pts.push_back(
1879  m_meshCurvedPts[c.ptid].pts[idx]);
1880  newidx++;
1881  }
1882  else
1883  {
1884  curvedpts.index.push_back(idxmap[idx]);
1885  }
1886  }
1887  }
1888  }
1889 
1890  // add xml information
1891  if(edgeinfo.size())
1892  {
1893  vCurved->SetAttribute("COMPRESSED",
1894  CompressData::GetCompressString());
1895  vCurved->SetAttribute("BITSIZE",
1896  CompressData::GetBitSizeStr());
1897 
1898  x = new TiXmlElement("E");
1899  std::string dataStr;
1900  CompressData::ZlibEncodeToBase64Str(edgeinfo,dataStr);
1901  x->LinkEndChild(new TiXmlText(dataStr));
1902  vCurved->LinkEndChild(x);
1903  }
1904 
1905  if(faceinfo.size())
1906  {
1907  vCurved->SetAttribute("COMPRESSED",
1908  CompressData::GetCompressString());
1909  vCurved->SetAttribute("BITSIZE",
1910  CompressData::GetBitSizeStr());
1911 
1912  x = new TiXmlElement("F");
1913  std::string dataStr;
1914  CompressData::ZlibEncodeToBase64Str(faceinfo,dataStr);
1915  x->LinkEndChild(new TiXmlText(dataStr));
1916  vCurved->LinkEndChild(x);
1917  }
1918 
1919  if(edgeinfo.size()||faceinfo.size())
1920  {
1921  x = new TiXmlElement("DATAPOINTS");
1922  x->SetAttribute("ID", curvedpts.id);
1923 
1924  TiXmlElement *subx = new TiXmlElement("INDEX");
1925  std::string dataStr;
1926  CompressData::ZlibEncodeToBase64Str(curvedpts.index,
1927  dataStr);
1928  subx->LinkEndChild(new TiXmlText(dataStr));
1929  x->LinkEndChild(subx);
1930 
1931  subx = new TiXmlElement("POINTS");
1932  CompressData::ZlibEncodeToBase64Str(curvedpts.pts,
1933  dataStr);
1934  subx->LinkEndChild(new TiXmlText(dataStr));
1935  x->LinkEndChild(subx);
1936 
1937  vCurved->LinkEndChild(x);
1938  }
1939  }
1940  else
1941  {
1942  for (vItCurve = m_meshCurved.begin();
1943  vItCurve != m_meshCurved.end();
1944  ++vItCurve)
1945  {
1946  MeshCurved c = vItCurve->second;
1947 
1948  bool IsEdge = boost::iequals(c.entitytype,"E");
1949  bool IsFace = boost::iequals(c.entitytype,"F");
1950 
1951  if((IsEdge&&vEdges.find(c.entityid) != vEdges.end())||
1952  (IsFace&&vFaces.find(c.entityid) != vFaces.end()))
1953  {
1954  x = new TiXmlElement(c.entitytype);
1955  x->SetAttribute("ID", c.id);
1956  if (IsEdge)
1957  {
1958  x->SetAttribute("EDGEID", c.entityid);
1959  }
1960  else
1961  {
1962  x->SetAttribute("FACEID", c.entityid);
1963  }
1964  x->SetAttribute("TYPE", c.type);
1965  x->SetAttribute("NUMPOINTS", c.npoints);
1966  y = new TiXmlText(c.data);
1967  x->LinkEndChild(y);
1968  vCurved->LinkEndChild(x);
1969  }
1970  }
1971  }
1972  }
1973 
1974  // Generate composites section comprising only those mesh entities
1975  // which belong to this partition.
1976  for (vIt = m_meshComposites.begin(); vIt != m_meshComposites.end(); ++vIt)
1977  {
1978  idxList.clear();
1979 
1980  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1981  {
1982  // Based on entity type, check if in this partition
1983  switch (vIt->second.type)
1984  {
1985  case 'V':
1986  if (vVertices.find(vIt->second.list[j]) == vVertices.end())
1987  {
1988  continue;
1989  }
1990  break;
1991  case 'E':
1992  if (vEdges.find(vIt->second.list[j]) == vEdges.end())
1993  {
1994  continue;
1995  }
1996  break;
1997  case 'F':
1998  if (vFaces.find(vIt->second.list[j]) == vFaces.end())
1999  {
2000  continue;
2001  }
2002  break;
2003  default:
2004  if (vElements.find(vIt->second.list[j]) == vElements.end())
2005  {
2006  continue;
2007  }
2008  break;
2009  }
2010 
2011  idxList.push_back(vIt->second.list[j]);
2012  }
2013 
2014  std::string vCompositeStr = ParseUtils::GenerateSeqString(idxList);
2015 
2016  if (vCompositeStr.length() > 0)
2017  {
2018  vComposites[vIt->first] = vIt->second;
2019  x = new TiXmlElement("C");
2020  x->SetAttribute("ID", vIt->first);
2021  vCompositeStr = "X[" + vCompositeStr + "]";
2022  vCompositeStr[0] = vIt->second.type;
2023  y = new TiXmlText(vCompositeStr.c_str());
2024  x->LinkEndChild(y);
2025  vComposite->LinkEndChild(x);
2026  }
2027  }
2028 
2029  idxList.clear();
2030  std::string vDomainListStr;
2031  for (unsigned int i = 0; i < m_domain.size(); ++i)
2032  {
2033  if (vComposites.find(m_domain[i]) != vComposites.end())
2034  {
2035  idxList.push_back(m_domain[i]);
2036  }
2037  }
2038  vDomainListStr = "C[" + ParseUtils::GenerateSeqString(idxList) + "]";
2039  TiXmlText* vDomainList = new TiXmlText(vDomainListStr);
2040  vDomain->LinkEndChild(vDomainList);
2041 
2042  vElmtGeometry->LinkEndChild(vVertex);
2043  if (m_dim >= 2)
2044  {
2045  vElmtGeometry->LinkEndChild(vEdge);
2046  }
2047  if (m_dim >= 3)
2048  {
2049  vElmtGeometry->LinkEndChild(vFace);
2050  }
2051  vElmtGeometry->LinkEndChild(vElement);
2052  if (m_dim >= 2)
2053  {
2054  vElmtGeometry->LinkEndChild(vCurved);
2055  }
2056  vElmtGeometry->LinkEndChild(vComposite);
2057  vElmtGeometry->LinkEndChild(vDomain);
2058 
2059  pNektar->LinkEndChild(vElmtGeometry);
2060 
2061  if (pSession->DefinesElement("Nektar/Conditions"))
2062  {
2063  std::set<int> vBndRegionIdList;
2064  TiXmlElement* vConditions = new TiXmlElement(*pSession->GetElement("Nektar/Conditions"));
2065  TiXmlElement* vBndRegions = vConditions->FirstChildElement("BOUNDARYREGIONS");
2066  TiXmlElement* vBndConditions = vConditions->FirstChildElement("BOUNDARYCONDITIONS");
2067  TiXmlElement* vItem;
2068 
2069  if (vBndRegions)
2070  {
2071  TiXmlElement* vNewBndRegions = new TiXmlElement("BOUNDARYREGIONS");
2072  vItem = vBndRegions->FirstChildElement();
2073  while (vItem)
2074  {
2075  std::string vSeqStr = vItem->FirstChild()->ToText()->Value();
2076  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
2077  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
2078  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
2079  std::vector<unsigned int> vSeq;
2080  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
2081 
2082  idxList.clear();
2083 
2084  for (unsigned int i = 0; i < vSeq.size(); ++i)
2085  {
2086  if (vComposites.find(vSeq[i]) != vComposites.end())
2087  {
2088  idxList.push_back(vSeq[i]);
2089  }
2090  }
2091  int p = atoi(vItem->Attribute("ID"));
2092 
2093  std::string vListStr = ParseUtils::GenerateSeqString(idxList);
2094 
2095  if (vListStr.length() == 0)
2096  {
2097  TiXmlElement* tmp = vItem;
2098  vItem = vItem->NextSiblingElement();
2099  vBndRegions->RemoveChild(tmp);
2100  }
2101  else
2102  {
2103  vListStr = "C[" + vListStr + "]";
2104  TiXmlText* vList = new TiXmlText(vListStr);
2105  TiXmlElement* vNewElement = new TiXmlElement("B");
2106  vNewElement->SetAttribute("ID", p);
2107  vNewElement->LinkEndChild(vList);
2108  vNewBndRegions->LinkEndChild(vNewElement);
2109  vBndRegionIdList.insert(p);
2110  vItem = vItem->NextSiblingElement();
2111  }
2112 
2113  // Store original order of boundary region.
2114  m_bndRegOrder[p] = vSeq;
2115 
2116  }
2117  vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
2118  }
2119 
2120  if (vBndConditions)
2121  {
2122  vItem = vBndConditions->FirstChildElement();
2123  while (vItem)
2124  {
2125  std::set<int>::iterator x;
2126  if ((x = vBndRegionIdList.find(atoi(vItem->Attribute("REF")))) != vBndRegionIdList.end())
2127  {
2128  vItem->SetAttribute("REF", *x);
2129  vItem = vItem->NextSiblingElement();
2130  }
2131  else
2132  {
2133  TiXmlElement* tmp = vItem;
2134  vItem = vItem->NextSiblingElement();
2135  vBndConditions->RemoveChild(tmp);
2136  }
2137  }
2138  }
2139  pNektar->LinkEndChild(vConditions);
2140  }
2141 
2142  // Distribute other sections of the XML to each process as is.
2143  TiXmlElement* vSrc = pSession->GetElement("Nektar")
2144  ->FirstChildElement();
2145  while (vSrc)
2146  {
2147  std::string vName = boost::to_upper_copy(vSrc->ValueStr());
2148  if (vName != "GEOMETRY" && vName != "CONDITIONS")
2149  {
2150  pNektar->LinkEndChild(new TiXmlElement(*vSrc));
2151  }
2152  vSrc = vSrc->NextSiblingElement();
2153  }
2154  }
2155 
2156  void MeshPartition::GetElementIDs(const int procid, std::vector<unsigned int> &elmtid)
2157  {
2158  BoostVertexIterator vertit, vertit_end;
2159 
2160  ASSERTL0(procid < m_localPartition.size(),"procid is less than the number of partitions");
2161 
2162  // Populate lists of elements, edges and vertices required.
2163  for ( boost::tie(vertit, vertit_end) = boost::vertices(m_localPartition[procid]);
2164  vertit != vertit_end;
2165  ++vertit)
2166  {
2167  elmtid.push_back(m_meshElements[m_localPartition[procid][*vertit].id].id);
2168  }
2169  }
2170 
2171  int MeshPartition::CalculateElementWeight(
2172  char elmtType,
2173  bool bndWeight,
2174  int na,
2175  int nb,
2176  int nc)
2177  {
2178  int weight = 0;
2179 
2180  switch (elmtType)
2181  {
2182  case 'A':
2183  weight = bndWeight ?
2184  StdTetData ::getNumberOfBndCoefficients(na, nb, nc) :
2185  StdTetData ::getNumberOfCoefficients (na, nb, nc);
2186  break;
2187  case 'R':
2188  weight = bndWeight ?
2189  StdPrismData::getNumberOfBndCoefficients(na, nb, nc) :
2190  StdPrismData::getNumberOfCoefficients (na, nb, nc);
2191  break;
2192  case 'H':
2193  weight = bndWeight ?
2194  StdHexData ::getNumberOfBndCoefficients(na, nb, nc) :
2195  StdHexData ::getNumberOfCoefficients (na, nb, nc);
2196  break;
2197  case 'P':
2198  weight = bndWeight ?
2199  StdPyrData ::getNumberOfBndCoefficients(na, nb, nc) :
2200  StdPyrData ::getNumberOfCoefficients (na, nb, nc);
2201  break;
2202  case 'Q':
2203  weight = bndWeight ?
2204  StdQuadData ::getNumberOfBndCoefficients(na, nb) :
2205  StdQuadData ::getNumberOfCoefficients (na, nb);
2206  break;
2207  case 'T':
2208  weight = bndWeight ?
2209  StdTriData ::getNumberOfBndCoefficients(na, nb) :
2210  StdTriData ::getNumberOfCoefficients (na, nb);
2211  break;
2212  case 'S':
2213  weight = bndWeight ?
2214  StdSegData ::getNumberOfBndCoefficients(na) :
2215  StdSegData ::getNumberOfCoefficients (na);
2216  break;
2217  case 'V':
2218  weight = 1;
2219  break;
2220  default:
2221  break;
2222  }
2223 
2224  return weight;
2225  }
2226  }
2227 }
Nektar::LibUtilities::MeshPartition::GetCompositeOrdering
void GetCompositeOrdering(CompositeOrdering &composites)
Definition: MeshPartition.cpp:173
Nektar::LibUtilities::MeshPartition::CalculateElementWeight
int CalculateElementWeight(char elmtType, bool bndWeight, int na, int nb, int nc)
Definition: MeshPartition.cpp:2171
Nektar::NekInt
boost::int32_t NekInt
Definition: NektarUnivTypeDefs.hpp:46
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
Nektar::ParseUtils::GenerateOrderedStringVector
static bool GenerateOrderedStringVector(const char *const str, std::vector< std::string > &vec)
Definition: ParseUtils.hpp:143
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
Nektar::LibUtilities::MeshPartition::m_vertexAttributes
std::map< std::string, std::string > m_vertexAttributes
Definition: MeshPartition.h:190
Nektar::LibUtilities::StdTriData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:120
Nektar::LibUtilities::StdPrismData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:286
Nektar::ParseUtils::GenerateOrderedVector
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:97
Nektar::LibUtilities::MeshCurvedInfo::ptoffset
NekInt64 ptoffset
the id of point data map (currently always 0 since we are using just one set).
Definition: MeshEntities.hpp:111
Nektar::LibUtilities::MeshPartition::m_meshComposites
std::map< int, MeshEntity > m_meshComposites
Definition: MeshPartition.h:188
Nektar::LibUtilities::MeshPartition::WriteLocalPartition
void WriteLocalPartition(SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:110
Nektar::LibUtilities::kPointsTypeStr
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:69
Foundations.hpp
Nektar::LibUtilities::StdHexData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:159
Nektar::LibUtilities::MeshPartition::m_meshVertices
std::map< int, MeshVertex > m_meshVertices
Definition: MeshPartition.h:182
FileSystem.h
Nektar::LibUtilities::MeshCurvedPts::index
std::vector< NekInt64 > index
id of this Point set
Definition: MeshEntities.hpp:130
Nektar::LibUtilities::MeshPartition::CheckPartitions
void CheckPartitions(int nParts, Array< OneD, int > &pPart)
Definition: MeshPartition.cpp:1298
Nektar::LibUtilities::MeshPartition::ReadConditions
void ReadConditions(const SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:986
Nektar::ParseUtils::GenerateSeqString
static std::string GenerateSeqString(const std::vector< unsigned int > &elmtids)
Definition: ParseUtils.hpp:159
Nektar::LibUtilities::StdPyrData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:266
Nektar::LibUtilities::MeshPartition::m_meshCurved
std::map< MeshCurvedKey, MeshCurved > m_meshCurved
Definition: MeshPartition.h:186
Nektar::LibUtilities::BndRegionOrdering
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshPartition.h:54
Nektar::LibUtilities::MeshPartition::PartitionGraphImpl
virtual void PartitionGraphImpl(int &nVerts, int &nVertConds, Nektar::Array< Nektar::OneD, int > &xadj, Nektar::Array< Nektar::OneD, int > &adjcy, Nektar::Array< Nektar::OneD, int > &vertWgt, Nektar::Array< Nektar::OneD, int > &vertSize, int &nparts, int &volume, Nektar::Array< Nektar::OneD, int > &part)=0
Nektar::LibUtilities::MeshPartition::BoostVertex
boost::graph_traits< BoostGraph >::vertex_descriptor BoostVertex
Definition: MeshPartition.h:162
Nektar::LibUtilities::MeshPartition::m_vertWeights
std::map< int, MultiWeight > m_vertWeights
Definition: MeshPartition.h:197
Nektar::LibUtilities::CompositeOrdering
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshPartition.h:53
Nektar::LibUtilities::StdTetData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:186
Nektar::LibUtilities::SessionReaderSharedPtr
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Nektar::ParseUtils::GenerateSeqVector
static bool GenerateSeqVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:79
Nektar::LibUtilities::MeshCurvedInfo::entityid
NekInt64 entityid
Id of this curved information.
Definition: MeshEntities.hpp:108
Nektar::LibUtilities::MeshPartition::NummodesPerField
std::map< std::string, NumModes > NummodesPerField
Definition: MeshPartition.h:177
Nektar::LibUtilities::MeshPartition::m_fieldNameToId
std::map< std::string, int > m_fieldNameToId
Definition: MeshPartition.h:196
Nektar::LibUtilities::MeshPartition::ReadExpansions
void ReadExpansions(const SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:189
Nektar::LibUtilities::MeshPartition::m_expansions
std::map< int, NummodesPerField > m_expansions
Definition: MeshPartition.h:194
Nektar::LibUtilities::MeshPartition::m_meshElements
std::map< int, MeshEntity > m_meshElements
Definition: MeshPartition.h:185
Nektar::LibUtilities::MeshCurvedInfo::npoints
NekInt64 npoints
The entity id corresponding to the global edge/curve.
Definition: MeshEntities.hpp:109
Nektar::LibUtilities::MeshPartitionFactory
LibUtilities::NekFactory< std::string, MeshPartition, const SessionReaderSharedPtr & > MeshPartitionFactory
Datatype of the NekFactory used to instantiate classes derived from the EquationSystem class...
Definition: MeshPartition.h:58
Nektar::LibUtilities::MeshPartition::CreateGraph
void CreateGraph(BoostSubGraph &pGraph)
Definition: MeshPartition.cpp:1103
Nektar::LibUtilities::MeshPartition::GetBndRegionOrdering
void GetBndRegionOrdering(BndRegionOrdering &composites)
Definition: MeshPartition.cpp:183
Nektar::LibUtilities::StdTriData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
Nektar::LibUtilities::MeshPartition::GetElementIDs
void GetElementIDs(const int procid, std::vector< unsigned int > &tmp)
Definition: MeshPartition.cpp:2156
Nektar::LibUtilities::MeshPartition::m_domain
std::vector< unsigned int > m_domain
Definition: MeshPartition.h:189
Nektar::LibUtilities::MeshPartition::m_localPartition
std::vector< BoostSubGraph > m_localPartition
Definition: MeshPartition.h:202
Nektar::LibUtilities::PortablePath
std::string PortablePath(const boost::filesystem::path &path)
create portable path on different platforms for boost::filesystem path
Definition: FileSystem.cpp:41
MeshPartition.h
Nektar::LibUtilities::MeshPartition::ReadGeometry
void ReadGeometry(const SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:292
Nektar::LibUtilities::CompressData::ZlibEncodeToBase64Str
int ZlibEncodeToBase64Str(std::vector< T > &in, std::string &out64)
Definition: CompressData.h:151
Nektar::LibUtilities::StdQuadData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:139
Nektar::LibUtilities::MeshCurvedInfo::ptype
NekInt64 ptype
point offset of data entry for this curve
Definition: MeshEntities.hpp:117
Nektar::LibUtilities::StdPyrData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:232
Nektar::iterator
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
Definition: MatrixOperations.cpp:96
Nektar::LibUtilities::MeshPartition::BoostAdjacencyIterator
boost::graph_traits< BoostGraph >::adjacency_iterator BoostAdjacencyIterator
Definition: MeshPartition.h:174
Nektar::LibUtilities::MeshPartition::m_meshEdges
std::map< int, MeshEntity > m_meshEdges
Definition: MeshPartition.h:183
Nektar::LibUtilities::StdPrismData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:297
Nektar::LibUtilities::MeshCurvedInfo::ptid
NekInt64 ptid
The number of points in this curved entity.
Definition: MeshEntities.hpp:110
Nektar::LibUtilities::StdQuadData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
Nektar::LibUtilities::MeshPartition::PartitionMesh
void PartitionMesh(int nParts, bool shared=false, bool overlapping=false)
Definition: MeshPartition.cpp:95
Nektar::LibUtilities::MeshPartition::BoostEdge
boost::graph_traits< BoostGraph >::edge_descriptor BoostEdge
Definition: MeshPartition.h:165
Nektar::LibUtilities::MeshPartition::PrintPartInfo
void PrintPartInfo(std::ostream &out)
Definition: MeshPartition.cpp:916
Nektar::LibUtilities::MeshPartition::m_meshFaces
std::map< int, MeshEntity > m_meshFaces
Definition: MeshPartition.h:184
Nektar::LibUtilities::MeshPartition::WriteAllPartitions
void WriteAllPartitions(SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:140
Nektar::LibUtilities::MeshPartition::BoostVertexIterator
boost::graph_traits< BoostGraph >::vertex_iterator BoostVertexIterator
Definition: MeshPartition.h:171
Nektar::LibUtilities::MeshCurvedPts::pts
std::vector< MeshVertex > pts
mapping to access pts value.
Definition: MeshEntities.hpp:134
SessionReader.h
Nektar::LibUtilities::MeshPartition::OutputPartition
void OutputPartition(SessionReaderSharedPtr &pSession, BoostSubGraph &pGraph, TiXmlElement *pGeometry)
Definition: MeshPartition.cpp:1329
Nektar::LibUtilities::MeshPartition::PartitionGraph
void PartitionGraph(BoostSubGraph &pGraph, int nParts, std::vector< BoostSubGraph > &pLocalPartition, bool overlapping=false)
Partition the graph.
Definition: MeshPartition.cpp:1160
Nektar::LibUtilities::StdHexData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:151
Nektar::LibUtilities::MeshPartition::MeshPartition
MeshPartition(const SessionReaderSharedPtr &pSession)
Definition: MeshPartition.cpp:78
ShapeType.hpp
Nektar::LibUtilities::MeshPartition::m_meshCurvedPts
std::map< int, MeshCurvedPts > m_meshCurvedPts
Definition: MeshPartition.h:187
ParseUtils.hpp
CompressData.h
Nektar::LibUtilities::GetMeshPartitionFactory
MeshPartitionFactory & GetMeshPartitionFactory()
Definition: MeshPartition.cpp:69
Nektar::LibUtilities::StdTetData::getNumberOfBndCoefficients
int getNumberOfBndCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:209
Nektar::LibUtilities::CompressData::ZlibDecodeFromBase64Str
int ZlibDecodeFromBase64Str(std::string &in64, std::vector< T > &out)
Definition: CompressData.h:243
Nektar::LibUtilities::NekFactory
Provides a generic Factory class.
Definition: NekFactory.hpp:116
Nektar::LibUtilities::MeshPartition::BoostSubGraph
boost::subgraph< BoostGraph > BoostSubGraph
Definition: MeshPartition.h:158
Generated on Wed May 4 2016 17:54:14 for Nektar++ by   doxygen 1.8.8