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