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 
55 #include <boost/algorithm/string.hpp>
56 #include <boost/graph/adjacency_list.hpp>
57 #include <boost/graph/adjacency_iterator.hpp>
58 #include <boost/graph/detail/edge.hpp>
59 #include <boost/format.hpp>
60 
61 namespace Nektar
62 {
63  namespace LibUtilities
64  {
65  MeshPartitionFactory& GetMeshPartitionFactory()
66  {
67  typedef Loki::SingletonHolder<MeshPartitionFactory,
68  Loki::CreateUsingNew,
69  Loki::NoDestroy > Type;
70  return Type::Instance();
71  }
72 
73  MeshPartition::MeshPartition(const LibUtilities::SessionReaderSharedPtr& pSession) :
74  m_numFields(0),
75  m_fieldNameToId(),
76  m_comm(pSession->GetComm()),
77  m_weightingRequired(false)
78  {
79  ReadConditions(pSession);
80  ReadGeometry(pSession);
81  ReadExpansions(pSession);
82  }
83 
84  MeshPartition::~MeshPartition()
85  {
86 
87  }
88 
89  void MeshPartition::PartitionMesh(bool shared)
90  {
91  ASSERTL0(m_comm->GetRowComm()->GetSize() > 1,
92  "Partitioning only necessary in parallel case.");
93  ASSERTL0(m_meshElements.size() >= m_comm->GetRowComm()->GetSize(),
94  "Too few elements for this many processes.");
95  m_shared = shared;
96 
97  if (m_weightingRequired) WeightElements();
98  CreateGraph(m_mesh);
99  PartitionGraph(m_mesh, m_localPartition);
100  }
101 
102  void MeshPartition::WriteLocalPartition(LibUtilities::SessionReaderSharedPtr& pSession)
103  {
104  TiXmlDocument vNew;
105  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
106  vNew.LinkEndChild(decl);
107 
108  TiXmlElement* vElmtNektar;
109  vElmtNektar = new TiXmlElement("NEKTAR");
110 
111  int rank = m_comm->GetRowComm()->GetRank();
112  OutputPartition(pSession, m_localPartition[rank], vElmtNektar);
113 
114  vNew.LinkEndChild(vElmtNektar);
115 
116  std::string dirname = pSession->GetSessionName() + "_xml";
117  fs::path pdirname(dirname);
118 
119  boost::format pad("P%1$07d.xml");
120  pad % rank;
121  fs::path pFilename(pad.str());
122 
123  if(!fs::is_directory(dirname))
124  {
125  fs::create_directory(dirname);
126  }
127 
128  fs::path fullpath = pdirname / pFilename;
129  vNew.SaveFile(PortablePath(fullpath));
130  }
131 
132  void MeshPartition::WriteAllPartitions(LibUtilities::SessionReaderSharedPtr& pSession)
133  {
134  for (int i = 0; i < m_comm->GetRowComm()->GetSize(); ++i)
135  {
136  TiXmlDocument vNew;
137  TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
138  vNew.LinkEndChild(decl);
139 
140  TiXmlElement* vElmtNektar;
141  vElmtNektar = new TiXmlElement("NEKTAR");
142 
143  OutputPartition(pSession, m_localPartition[i], vElmtNektar);
144 
145  vNew.LinkEndChild(vElmtNektar);
146 
147  std::string dirname = pSession->GetSessionName() + "_xml";
148  fs::path pdirname(dirname);
149 
150  boost::format pad("P%1$07d.xml");
151  pad % i;
152  fs::path pFilename(pad.str());
153 
154  fs::path fullpath = pdirname / pFilename;
155 
156  if(!fs::is_directory(dirname))
157  {
158  fs::create_directory(dirname);
159  }
160 
161  vNew.SaveFile(PortablePath(fullpath));
162  }
163  }
164 
165  void MeshPartition::GetCompositeOrdering(CompositeOrdering &composites)
166  {
167  std::map<int, MeshEntity>::iterator it;
168  for (it = m_meshComposites.begin();
169  it != m_meshComposites.end(); ++it)
170  {
171  composites[it->first] = it->second.list;
172  }
173  }
174 
175  void MeshPartition::GetBndRegionOrdering(BndRegionOrdering &bndRegs)
176  {
177  bndRegs = m_bndRegOrder;
178  }
179 
180 
181  void MeshPartition::ReadExpansions(const LibUtilities::SessionReaderSharedPtr& pSession)
182  {
183  // Find the Expansions tag
184  TiXmlElement *expansionTypes = pSession->GetElement("Nektar/Expansions");
185 
186  // Find the Expansion type
187  TiXmlElement *expansion = expansionTypes->FirstChildElement();
188  std::string expType = expansion->Value();
189 
190 
191  if(expType != "E")
192  {
193  ASSERTL0(false,"Expansion type not defined or not supported at the moment");
194  }
195 
196  /// Expansiontypes will contain plenty of data,
197  /// where relevant at this stage are composite
198  /// ID(s) that this expansion type describes,
199  /// nummodes and a list of fields that this
200  /// expansion relates to. If this does not exist
201  /// the variable is only set to "DefaultVar".
202 
203  while (expansion)
204  {
205  std::vector<unsigned int> composite;
206  std::vector<unsigned int> nummodes;
207  std::vector<std::string> fieldName;
208 
209  const char *nModesStr = expansion->Attribute("NUMMODES");
210  ASSERTL0(nModesStr,"NUMMODES was not defined in EXPANSION section of input");
211  std::string numModesStr = nModesStr;
212  bool valid = ParseUtils::GenerateOrderedVector(numModesStr.c_str(), nummodes);
213  ASSERTL0(valid, "Unable to correctly parse the number of modes.");
214 
215  if (nummodes.size() == 1)
216  {
217  for (int i = 1; i < m_dim; i++)
218  {
219  nummodes.push_back( nummodes[0] );
220  }
221  }
222  ASSERTL0(nummodes.size() == m_dim,"Number of modes should match mesh dimension");
223 
224 
225  const char *fStr = expansion->Attribute("FIELDS");
226  if(fStr)
227  {
228  std::string fieldStr = fStr;
229  bool valid = ParseUtils::GenerateOrderedStringVector(fieldStr.c_str(),fieldName);
230  ASSERTL0(valid,"Unable to correctly parse the field string in ExpansionTypes.");
231 
232  for (int i = 0; i < fieldName.size(); ++i)
233  {
234  if (m_fieldNameToId.count(fieldName[i]) == 0)
235  {
236  int k = m_fieldNameToId.size();
237  m_fieldNameToId[ fieldName[i] ] = k;
238  m_numFields++;
239  }
240  }
241  }
242  else
243  {
244  fieldName.push_back("DefaultVar");
245  int k = m_fieldNameToId.size();
246 
247  if (m_fieldNameToId.count("DefaultVar") == 0)
248  {
249  ASSERTL0(k == 0,
250  "Omitting field variables and explicitly listing " \
251  "them in different ExpansionTypes is wrong practise");
252 
253  m_fieldNameToId[ "DefaultVar" ] = k;
254  m_numFields++;
255  }
256  }
257 
258  std::string compositeStr = expansion->Attribute("COMPOSITE");
259  ASSERTL0(compositeStr.length() > 3, "COMPOSITE must be specified in expansion definition");
260  int beg = compositeStr.find_first_of("[");
261  int end = compositeStr.find_first_of("]");
262  std::string compositeListStr = compositeStr.substr(beg+1,end-beg-1);
263  bool parseGood = ParseUtils::GenerateSeqVector(compositeListStr.c_str(), composite);
264  ASSERTL0(parseGood && !composite.empty(),
265  (std::string("Unable to read composite index range: ") + compositeListStr).c_str());
266 
267 
268  // construct mapping (field name, CompositeID) -> nummodes
269  for (int i = 0; i < composite.size(); ++i)
270  {
271  for (int j = 0; j < fieldName.size(); j++)
272  {
273  m_expansions[composite[i]][fieldName[j]] = nummodes;
274  }
275  }
276 
277  expansion = expansion->NextSiblingElement("E");
278  }
279  }
280 
281 
282 
283 
284  void MeshPartition::ReadGeometry(const LibUtilities::SessionReaderSharedPtr& pSession)
285  {
286  TiXmlElement* x;
287  TiXmlElement *vGeometry, *vSubElement;
288  int i;
289 
290  vGeometry = pSession->GetElement("Nektar/Geometry");
291  m_dim = atoi(vGeometry->Attribute("DIM"));
292 
293  // Read mesh vertices
294  vSubElement = pSession->GetElement("Nektar/Geometry/Vertex");
295 
296  // Retrieve any VERTEX attributes specifying mesh transforms
297  std::string attr[] = {"XSCALE", "YSCALE", "ZSCALE",
298  "XMOVE", "YMOVE", "ZMOVE" };
299  for (i = 0; i < 6; ++i)
300  {
301  const char *val = vSubElement->Attribute(attr[i].c_str());
302  if (val)
303  {
304  m_vertexAttributes[attr[i]] = std::string(val);
305  }
306  }
307 
308  x = vSubElement->FirstChildElement();
309  i = 0;
310  if (x->FirstAttribute())
311  {
312  i = x->FirstAttribute()->IntValue();
313  }
314  while(x)
315  {
316  TiXmlAttribute* y = x->FirstAttribute();
317  ASSERTL0(y, "Failed to get attribute.");
318  MeshVertex v;
319  v.id = y->IntValue();
320  ASSERTL0(v.id == i++, "Vertex IDs not sequential.");
321  std::vector<std::string> vCoords;
322  std::string vCoordStr = x->FirstChild()->ToText()->Value();
323  boost::split(vCoords, vCoordStr, boost::is_any_of("\t "));
324  v.x = atof(vCoords[0].c_str());
325  v.y = atof(vCoords[1].c_str());
326  v.z = atof(vCoords[2].c_str());
327  m_meshVertices[v.id] = v;
328  x = x->NextSiblingElement();
329  }
330 
331  // Read mesh edges
332  if (m_dim >= 2)
333  {
334  vSubElement = pSession->GetElement("Nektar/Geometry/Edge");
335  ASSERTL0(vSubElement, "Cannot read edges");
336  x = vSubElement->FirstChildElement();
337  i = 0;
338  while(x)
339  {
340  TiXmlAttribute* y = x->FirstAttribute();
341  ASSERTL0(y, "Failed to get attribute.");
342  MeshEntity e;
343  e.id = y->IntValue();
344  e.type = 'E';
345  ASSERTL0(e.id == i++, "Edge IDs not sequential.");
346  std::vector<std::string> vVertices;
347  std::string vVerticesString = x->FirstChild()->ToText()->Value();
348  boost::split(vVertices, vVerticesString, boost::is_any_of("\t "));
349  e.list.push_back(atoi(vVertices[0].c_str()));
350  e.list.push_back(atoi(vVertices[1].c_str()));
351  m_meshEdges[e.id] = e;
352  x = x->NextSiblingElement();
353  }
354  }
355 
356  // Read mesh faces
357  if (m_dim == 3)
358  {
359  vSubElement = pSession->GetElement("Nektar/Geometry/Face");
360  ASSERTL0(vSubElement, "Cannot read faces.");
361  x = vSubElement->FirstChildElement();
362  i = 0;
363  while(x)
364  {
365  TiXmlAttribute* y = x->FirstAttribute();
366  ASSERTL0(y, "Failed to get attribute.");
367  MeshEntity f;
368  f.id = y->IntValue();
369  f.type = x->Value()[0];
370  ASSERTL0(f.id == i++, "Face IDs not sequential.");
371  std::vector<std::string> vEdges;
372  std::string vEdgeStr = x->FirstChild()->ToText()->Value();
373  boost::split(vEdges, vEdgeStr, boost::is_any_of("\t "));
374  for (int i = 0; i < vEdges.size(); ++i)
375  {
376  f.list.push_back(atoi(vEdges[i].c_str()));
377  }
378  m_meshFaces[f.id] = f;
379  x = x->NextSiblingElement();
380  }
381  }
382 
383  // Read mesh elements
384  vSubElement = pSession->GetElement("Nektar/Geometry/Element");
385  ASSERTL0(vSubElement, "Cannot read elements.");
386  x = vSubElement->FirstChildElement();
387  i = 0;
388  while(x)
389  {
390  TiXmlAttribute* y = x->FirstAttribute();
391  ASSERTL0(y, "Failed to get attribute.");
392  MeshEntity e;
393  e.id = y->IntValue();
394  ASSERTL0(e.id == i++, "Element IDs not sequential.");
395  std::vector<std::string> vItems;
396  std::string vItemStr = x->FirstChild()->ToText()->Value();
397  boost::split(vItems, vItemStr, boost::is_any_of("\t "));
398  for (int i = 0; i < vItems.size(); ++i)
399  {
400  e.list.push_back(atoi(vItems[i].c_str()));
401  }
402  e.type = x->Value()[0];
403  m_meshElements[e.id] = e;
404  x = x->NextSiblingElement();
405  }
406 
407  // Read mesh curves
408  if (pSession->DefinesElement("Nektar/Geometry/Curved"))
409  {
410  vSubElement = pSession->GetElement("Nektar/Geometry/Curved");
411  x = vSubElement->FirstChildElement();
412  i = 0;
413  while(x)
414  {
415  MeshCurved c;
416  ASSERTL0(x->Attribute("ID", &c.id),
417  "Failed to get attribute ID");
418  c.type = std::string(x->Attribute("TYPE"));
419  ASSERTL0(!c.type.empty(),
420  "Failed to get attribute TYPE");
421  ASSERTL0(x->Attribute("NUMPOINTS", &c.npoints),
422  "Failed to get attribute NUMPOINTS");
423  c.data = x->FirstChild()->ToText()->Value();
424  c.entitytype = x->Value()[0];
425  if (c.entitytype == "E")
426  {
427  ASSERTL0(x->Attribute("EDGEID", &c.entityid),
428  "Failed to get attribute EDGEID");
429  }
430  else if (c.entitytype == "F")
431  {
432  ASSERTL0(x->Attribute("FACEID", &c.entityid),
433  "Failed to get attribute FACEID");
434  }
435  else
436  {
437  ASSERTL0(false, "Unknown curve type.");
438  }
439  m_meshCurved[std::make_pair(c.entitytype, c.id)] = c;
440  x = x->NextSiblingElement();
441  }
442  }
443 
444  // Read composites
445  vSubElement = pSession->GetElement("Nektar/Geometry/Composite");
446  ASSERTL0(vSubElement, "Cannot read composites.");
447  x = vSubElement->FirstChildElement();
448  i = 0;
449  while(x)
450  {
451  TiXmlAttribute* y = x->FirstAttribute();
452  ASSERTL0(y, "Failed to get attribute.");
453  MeshEntity c;
454  c.id = y->IntValue();
455  std::string vSeqStr = x->FirstChild()->ToText()->Value();
456  c.type = vSeqStr[0];
457  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
458  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
459  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
460 
461  std::vector<unsigned int> vSeq;
462  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
463 
464  for (int i = 0; i < vSeq.size(); ++i)
465  {
466  c.list.push_back(vSeq[i]);
467  }
468  m_meshComposites[c.id] = c;
469  x = x->NextSiblingElement();
470  }
471 
472  // Read Domain
473  vSubElement = pSession->GetElement("Nektar/Geometry/Domain");
474  ASSERTL0(vSubElement, "Cannot read domain");
475  std::string vSeqStr = vSubElement->FirstChild()->ToText()->Value();
476  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
477  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
478  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
479  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), m_domain);
480  }
481 
482  void MeshPartition::PrintPartInfo(std::ostream &out)
483  {
484  int nElmt = boost::num_vertices(m_mesh);
485  int nPart = m_comm->GetRowComm()->GetSize();
486 
487  out << "# Partition information:" << std::endl;
488  out << "# No. elements : " << nElmt << std::endl;
489  out << "# No. partitions: " << nPart << std::endl;
490 
491  BoostVertexIterator vertit, vertit_end;
492  std::vector<int> partElmtCount(nPart, 0);
493  std::vector<int> partLocCount (nPart, 0);
494  std::vector<int> partBndCount (nPart, 0);
495 
496  std::map<int, int> elmtSizes;
497  std::map<int, int> elmtBndSizes;
498 
499  for (unsigned int i = 0; i < m_domain.size(); ++i)
500  {
501  int cId = m_domain[i];
502  NummodesPerField npf = m_expansions[cId];
503 
504  for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
505  {
506  ASSERTL0(it->second.size() == m_dim,
507  " Number of directional" \
508  " modes in expansion spec for composite id = " +
509  boost::lexical_cast<std::string>(cId) +
510  " and field " +
511  boost::lexical_cast<std::string>(it->first) +
512  " does not correspond to mesh dimension");
513 
514  int na = it->second[0];
515  int nb = it->second[1];
516  int nc = 0;
517  if (m_dim == 3)
518  {
519  nc = it->second[2];
520  }
521 
522  int weight = 0;
523  int bndWeight = 0;
524  switch (m_meshComposites[cId].type)
525  {
526  case 'A':
527  weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
528  bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
529  break;
530  case 'R':
531  weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
532  bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
533  break;
534  case 'H':
535  weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
536  bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
537  break;
538  case 'P':
539  weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
540  bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
541  break;
542  case 'Q':
543  weight = StdQuadData::getNumberOfCoefficients(na, nb);
544  bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
545  break;
546  case 'T':
547  weight = StdTriData::getNumberOfCoefficients(na, nb);
548  bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
549  break;
550  case 'S':
551  weight = StdSegData::getNumberOfCoefficients(na);
552  bndWeight = StdSegData::getNumberOfBndCoefficients(na);
553  break;
554  case 'V':
555  weight = 1;
556  bndWeight = 1;
557  break;
558  default:
559  break;
560  }
561 
562  for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
563  {
564  int elid = m_meshComposites[cId].list[j];
565  elmtSizes[elid] = weight;
566  elmtBndSizes[elid] = bndWeight;
567  }
568  }
569  }
570 
571  for (boost::tie(vertit, vertit_end) = boost::vertices(m_mesh);
572  vertit != vertit_end; ++vertit)
573  {
574  int partId = m_mesh[*vertit].partition;
575  partElmtCount[partId]++;
576  partLocCount [partId] += elmtSizes[m_mesh[*vertit].id];
577  partBndCount [partId] += elmtBndSizes[m_mesh[*vertit].id];
578  }
579 
580  for (int i = 0; i < nPart; ++i)
581  {
582  out << i << " " << partElmtCount[i] << " " << partLocCount[i] << " " << partBndCount[i] << std::endl;
583  }
584  }
585 
586  void MeshPartition::ReadConditions(const SessionReaderSharedPtr& pSession)
587  {
588  if (!pSession->DefinesElement("Nektar/Conditions/SolverInfo"))
589  {
590  // No SolverInfo = no change of default action to weight
591  // mesh graph.
592  return;
593  }
594 
595  TiXmlElement* solverInfoElement =
596  pSession->GetElement("Nektar/Conditions/SolverInfo");
597 
598  TiXmlElement* solverInfo =
599  solverInfoElement->FirstChildElement("I");
600  ASSERTL0(solverInfo, "Cannot read SolverInfo tags");
601 
602  while (solverInfo)
603  {
604  // read the property name
605  ASSERTL0(solverInfo->Attribute("PROPERTY"),
606  "Missing PROPERTY attribute in solver info "
607  "section. ");
608  std::string solverProperty =
609  solverInfo->Attribute("PROPERTY");
610  ASSERTL0(!solverProperty.empty(),
611  "Solver info properties must have a non-empty "
612  "name. ");
613  // make sure that solver property is capitalised
614  std::string solverPropertyUpper =
615  boost::to_upper_copy(solverProperty);
616 
617 
618  // read the value
619  ASSERTL0(solverInfo->Attribute("VALUE"),
620  "Missing VALUE attribute in solver info section. ");
621  std::string solverValue = solverInfo->Attribute("VALUE");
622  ASSERTL0(!solverValue.empty(),
623  "Solver info properties must have a non-empty value");
624  // make sure that property value is capitalised
625  std::string propertyValueUpper =
626  boost::to_upper_copy(solverValue);
627 
628  if (solverPropertyUpper == "WEIGHTPARTITIONS")
629  {
630  if (propertyValueUpper != "UNIFORM")
631  {
632  m_weightingRequired = true;
633  }
634  return;
635  }
636  solverInfo = solverInfo->NextSiblingElement("I");
637  }
638  }
639 
640 
641  /*
642  * Calculate element weights based on
643  * - element type (Q,T,H,P,R,A)
644  * - nummodes in expansion which this element belongs to via composite.
645  *
646  * For each element we prepare two vertex weightings, one associated
647  * with the number of matrix elements associated with it (to balance
648  * matrix multiplication work) and another associated
649  * with all work which scales linearly with the number of its
650  * coefficients: communication, vector updates etc.
651  *
652  * \todo Refactor this code to explicitly represent performance model
653  * and flexibly generate graph vertex weights depending on perf data.
654  */
655  void MeshPartition::WeightElements()
656  {
657  std::vector<unsigned int> weight(m_numFields, 1);
658  for (int i = 0; i < m_meshElements.size(); ++i)
659  {
660  m_vertWeights.push_back( weight );
661  }
662 
663  for (unsigned int i = 0; i < m_domain.size(); ++i)
664  {
665  int cId = m_domain[i];
666  NummodesPerField npf = m_expansions[cId];
667 
668  for (NummodesPerField::iterator it = npf.begin(); it != npf.end(); ++it)
669  {
670  ASSERTL0(it->second.size() == m_dim,
671  " Number of directional" \
672  " modes in expansion spec for composite id = " +
673  boost::lexical_cast<std::string>(cId) +
674  " and field " +
675  boost::lexical_cast<std::string>(it->first) +
676  " does not correspond to mesh dimension");
677 
678  int na = it->second[0];
679  int nb = 0;
680  int nc = 0;
681  if (m_dim >= 2)
682  {
683  nb = it->second[1];
684  }
685  if (m_dim == 3)
686  {
687  nc = it->second[2];
688  }
689 
690  int weight = 0;
691  int bndWeight = 0;
692  switch (m_meshComposites[cId].type)
693  {
694  case 'A':
695  weight = StdTetData::getNumberOfCoefficients(na, nb, nc);
696  bndWeight = StdTetData::getNumberOfBndCoefficients(na, nb, nc);
697  break;
698  case 'R':
699  weight = StdPrismData::getNumberOfCoefficients(na, nb, nc);
700  bndWeight = StdPrismData::getNumberOfBndCoefficients(na, nb, nc);
701  break;
702  case 'H':
703  weight = StdHexData::getNumberOfCoefficients(na, nb, nc);
704  bndWeight = StdHexData::getNumberOfBndCoefficients(na, nb, nc);
705  break;
706  case 'P':
707  weight = StdPyrData::getNumberOfCoefficients(na, nb, nc);
708  bndWeight = StdPyrData::getNumberOfBndCoefficients(na, nb, nc);
709  break;
710  case 'Q':
711  weight = StdQuadData::getNumberOfCoefficients(na, nb);
712  bndWeight = StdQuadData::getNumberOfBndCoefficients(na, nb);
713  break;
714  case 'T':
715  weight = StdTriData::getNumberOfCoefficients(na, nb);
716  bndWeight = StdTriData::getNumberOfBndCoefficients(na, nb);
717  break;
718  case 'S':
719  weight = StdSegData::getNumberOfCoefficients(na);
720  bndWeight = StdSegData::getNumberOfBndCoefficients(na);
721  break;
722  case 'V':
723  weight = 1;
724  bndWeight = 1;
725  break;
726  default:
727  break;
728  }
729 
730  for (unsigned int j = 0; j < m_meshComposites[cId].list.size(); ++j)
731  {
732  int elmtId = m_meshComposites[cId].list[j];
733  m_vertWeights[elmtId][ m_fieldNameToId[it->first]] = bndWeight;
734  //m_vertWeights[elmtId][ m_fieldNameToId[ it->first ] * 2 + 1 ] = weight*weight;
735  }
736  }
737  } // for i
738  }
739 
740  void MeshPartition::CreateGraph(BoostSubGraph& pGraph)
741  {
742  // Maps edge/face to first mesh element id.
743  // On locating second mesh element id, graph edge is created instead.
744  std::map<int, int> vGraphEdges;
745 
746  for (unsigned int i = 0; i < m_meshElements.size(); ++i)
747  {
748  int p = m_meshElements[i].id;
749  BoostVertex v = boost::add_vertex(pGraph);
750  pGraph[v].id = p;
751  pGraph[v].partition = 0;
752  if (m_weightingRequired)
753  {
754  pGraph[v].weight = m_vertWeights[i];
755  }
756 
757  // Process element entries and add graph edges
758  for (unsigned j = 0; j < m_meshElements[i].list.size(); ++j)
759  {
760  int eId = m_meshElements[i].list[j];
761 
762  // Look to see if we've examined this edge/face before
763  // If so, we've got both graph vertices so add edge
764  if (vGraphEdges.find(eId) != vGraphEdges.end())
765  {
766  BoostEdge e = boost::add_edge( p, vGraphEdges[eId], pGraph).first;
767  pGraph[e].id = eId;
768  }
769  else
770  {
771  vGraphEdges[eId] = p;
772  }
773  }
774  }
775  }
776 
777  void MeshPartition::PartitionGraph(BoostSubGraph& pGraph,
778  std::vector<BoostSubGraph>& pLocalPartition)
779  {
780  int i;
781  int nGraphVerts = boost::num_vertices(pGraph);
782  int nGraphEdges = boost::num_edges(pGraph);
783 
784  // Convert boost graph into CSR format
785  BoostVertexIterator vertit, vertit_end;
786  Array<OneD, int> part(nGraphVerts,0);
787 
788  if (m_comm->GetRowComm()->TreatAsRankZero())
789  {
790  int acnt = 0;
791  int vcnt = 0;
792  int nWeight = nGraphVerts;
793  BoostAdjacencyIterator adjvertit, adjvertit_end;
794  Array<OneD, int> xadj(nGraphVerts+1,0);
795  Array<OneD, int> adjncy(2*nGraphEdges);
796  Array<OneD, int> vwgt(nWeight, 1);
797  Array<OneD, int> vsize(nGraphVerts, 1);
798  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
799  vertit != vertit_end;
800  ++vertit)
801  {
802  for ( boost::tie(adjvertit, adjvertit_end) = boost::adjacent_vertices(*vertit,pGraph);
803  adjvertit != adjvertit_end;
804  ++adjvertit)
805  {
806  adjncy[acnt++] = *adjvertit;
807 
808  }
809  xadj[++vcnt] = acnt;
810 
811  if (m_weightingRequired)
812  {
813  vwgt[pGraph[*vertit].id ] = pGraph[*vertit].weight[0];
814  }
815  else
816  {
817  vwgt[pGraph[*vertit].id] = 1;
818  }
819  }
820 
821  // Call Metis and partition graph
822  int npart = m_comm->GetRowComm()->GetSize();
823  int vol = 0;
824 
825  try
826  {
827  //////////////////////////////////////////////////////
828  // On a cartesian communicator do mesh partiotion just on the first column
829  // so there is no doubt the partitions are all the same in all the columns
830  if(m_comm->GetColumnComm()->GetRank() == 0)
831  {
832  // Attempt partitioning using METIS.
833  int ncon = 1;
834  PartitionGraphImpl(nGraphVerts, ncon, xadj, adjncy, vwgt, vsize, npart, vol, part);
835 
836  // Check METIS produced a valid partition and fix if not.
837  CheckPartitions(part);
838  if (!m_shared)
839  {
840  // distribute among columns
841  for (i = 1; i < m_comm->GetColumnComm()->GetSize(); ++i)
842  {
843  m_comm->GetColumnComm()->Send(i, part);
844  }
845  }
846  }
847  else
848  {
849  m_comm->GetColumnComm()->Recv(0, part);
850  }
851  if (!m_shared)
852  {
853  m_comm->GetColumnComm()->Block();
854 
855  //////////////////////////////////
856  // distribute among rows
857  for (i = 1; i < m_comm->GetRowComm()->GetSize(); ++i)
858  {
859  m_comm->GetRowComm()->Send(i, part);
860  }
861  }
862  }
863  catch (...)
864  {
865  NEKERROR(ErrorUtil::efatal,
866  "Error in calling metis to partition graph.");
867  }
868  }
869  else
870  {
871  m_comm->GetRowComm()->Recv(0, part);
872  }
873 
874  // Create boost subgraph for this process's partitions
875  int nCols = m_comm->GetRowComm()->GetSize();
876  pLocalPartition.resize(nCols);
877  for (i = 0; i < nCols; ++i)
878  {
879  pLocalPartition[i] = pGraph.create_subgraph();
880  }
881 
882  // Populate subgraph
883  i = 0;
884  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
885  vertit != vertit_end;
886  ++vertit, ++i)
887  {
888  pGraph[*vertit].partition = part[i];
889  pGraph[*vertit].partid = boost::num_vertices(pLocalPartition[part[i]]);
890  boost::add_vertex(i, pLocalPartition[part[i]]);
891  }
892  }
893 
894 
895  void MeshPartition::CheckPartitions(Array<OneD, int> &pPart)
896  {
897  unsigned int i = 0;
898  unsigned int cnt = 0;
899  const unsigned int npart = m_comm->GetRowComm()->GetSize();
900  bool valid = true;
901 
902  // Check that every process has at least one element assigned
903  for (i = 0; i < npart; ++i)
904  {
905  cnt = std::count(pPart.begin(), pPart.end(), i);
906  if (cnt == 0)
907  {
908  valid = false;
909  }
910  }
911 
912  // If METIS produced an invalid partition, repartition naively.
913  // Elements are assigned to processes in a round-robin fashion.
914  // It is assumed that METIS failure only occurs when the number of
915  // elements is approx. the number of processes, so this approach
916  // should not be too inefficient communication-wise.
917  if (!valid)
918  {
919  for (i = 0; i < pPart.num_elements(); ++i)
920  {
921  pPart[i] = i % npart;
922  }
923  }
924  }
925 
926 
927  void MeshPartition::OutputPartition(
928  LibUtilities::SessionReaderSharedPtr& pSession,
929  BoostSubGraph& pGraph,
930  TiXmlElement* pNektar)
931  {
932  // Write Geometry data
933  std::string vDim = pSession->GetElement("Nektar/Geometry")->Attribute("DIM");
934  std::string vSpace = pSession->GetElement("Nektar/Geometry")->Attribute("SPACE");
935  std::string vPart = boost::lexical_cast<std::string>(pGraph[*boost::vertices(pGraph).first].partition);
936  TiXmlElement* vElmtGeometry = new TiXmlElement("GEOMETRY");
937  vElmtGeometry->SetAttribute("DIM", vDim);
938  vElmtGeometry->SetAttribute("SPACE", vSpace);
939  vElmtGeometry->SetAttribute("PARTITION", vPart);
940 
941  TiXmlElement *vVertex = new TiXmlElement("VERTEX");
942  TiXmlElement *vEdge = new TiXmlElement("EDGE");
943  TiXmlElement *vFace = new TiXmlElement("FACE");
944  TiXmlElement *vElement = new TiXmlElement("ELEMENT");
945  TiXmlElement *vCurved = new TiXmlElement("CURVED");
946  TiXmlElement *vComposite = new TiXmlElement("COMPOSITE");
947  TiXmlElement *vDomain = new TiXmlElement("DOMAIN");
948 
949  TiXmlElement *x;
950  TiXmlText *y;
951 
952  BoostVertexIterator vertit, vertit_end;
953  int id;
954 
955  std::map<int, MeshEntity> vComposites;
956  std::map<int, MeshEntity> vElements;
957  std::map<int, MeshEntity> vEdges;
958  std::map<int, MeshEntity> vFaces;
959  std::map<int, MeshVertex> vVertices;
960  std::map<int, MeshEntity>::iterator vIt;
961  std::map<int, MeshVertex>::iterator vVertIt;
962  std::map<std::string, std::string>::iterator vAttrIt;
963 
964  // Populate lists of elements, edges and vertices required.
965  for ( boost::tie(vertit, vertit_end) = boost::vertices(pGraph);
966  vertit != vertit_end;
967  ++vertit)
968  {
969  id = pGraph[*vertit].id;
970  vElements[id] = m_meshElements[pGraph[*vertit].id];
971  }
972 
973  std::map<int, MeshEntity> * vNext = &vElements;
974  switch (m_dim)
975  {
976  case 3:
977  {
978  // Compile list of faces
979  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
980  {
981  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
982  {
983  id = vIt->second.list[j];
984  vFaces[id] = m_meshFaces[id];
985  }
986  }
987  vNext = &vFaces;
988  }
989  case 2:
990  {
991  // Compile list of edges
992  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
993  {
994  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
995  {
996  id = vIt->second.list[j];
997  vEdges[id] = m_meshEdges[id];
998  }
999  }
1000  vNext = &vEdges;
1001  }
1002  case 1:
1003  {
1004  // Compile list of vertices
1005  for (vIt = vNext->begin(); vIt != vNext->end(); vIt++)
1006  {
1007  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1008  {
1009  id = vIt->second.list[j];
1010  vVertices[id] = m_meshVertices[id];
1011  }
1012  }
1013  }
1014  }
1015 
1016  // Generate XML data for these mesh entities
1017  for (vVertIt = vVertices.begin(); vVertIt != vVertices.end(); vVertIt++)
1018  {
1019  x = new TiXmlElement("V");
1020  x->SetAttribute("ID", vVertIt->first);
1021  std::stringstream vCoords;
1022  vCoords.precision(12);
1023  vCoords << std::setw(15) << vVertIt->second.x
1024  << std::setw(15) << vVertIt->second.y
1025  << std::setw(15) << vVertIt->second.z << " ";
1026  y = new TiXmlText(vCoords.str());
1027  x->LinkEndChild(y);
1028  vVertex->LinkEndChild(x);
1029  }
1030 
1031  // Apply transformation attributes to VERTEX section
1032  for (vAttrIt = m_vertexAttributes.begin();
1033  vAttrIt != m_vertexAttributes.end();
1034  ++ vAttrIt)
1035  {
1036  vVertex->SetAttribute(vAttrIt->first, vAttrIt->second);
1037  }
1038 
1039  if (m_dim >= 2)
1040  {
1041  for (vIt = vEdges.begin(); vIt != vEdges.end(); vIt++)
1042  {
1043  x = new TiXmlElement("E");
1044  x->SetAttribute("ID", vIt->first);
1045  std::stringstream vVertices;
1046  vVertices << std::setw(10) << vIt->second.list[0]
1047  << std::setw(10) << vIt->second.list[1] << " ";
1048  y = new TiXmlText(vVertices.str());
1049  x->LinkEndChild(y);
1050  vEdge->LinkEndChild(x);
1051  }
1052  }
1053 
1054  if (m_dim >= 3)
1055  {
1056  for (vIt = vFaces.begin(); vIt != vFaces.end(); vIt++)
1057  {
1058  std::string vType("F");
1059  vType[0] = vIt->second.type;
1060  x = new TiXmlElement(vType);
1061  x->SetAttribute("ID", vIt->first);
1062  std::stringstream vListStr;
1063  for (unsigned int i = 0; i < vIt->second.list.size(); ++i)
1064  {
1065  vListStr << std::setw(10) << vIt->second.list[i];
1066  }
1067  vListStr << " ";
1068  y = new TiXmlText(vListStr.str());
1069  x->LinkEndChild(y);
1070  vFace->LinkEndChild(x);
1071  }
1072  }
1073 
1074  for (vIt = vElements.begin(); vIt != vElements.end(); vIt++)
1075  {
1076  std::string vType("T");
1077  vType[0] = vIt->second.type;
1078  x = new TiXmlElement(vType.c_str());
1079  x->SetAttribute("ID", vIt->first);
1080  std::stringstream vEdges;
1081  for (unsigned i = 0; i < vIt->second.list.size(); ++i)
1082  {
1083  vEdges << std::setw(10) << vIt->second.list[i];
1084  }
1085  vEdges << " ";
1086  y = new TiXmlText(vEdges.str());
1087  x->LinkEndChild(y);
1088  vElement->LinkEndChild(x);
1089  }
1090 
1091  if (m_dim >= 2)
1092  {
1093  std::map<MeshCurvedKey, MeshCurved>::const_iterator vItCurve;
1094  for (vItCurve = m_meshCurved.begin();
1095  vItCurve != m_meshCurved.end();
1096  ++vItCurve)
1097  {
1098  MeshCurved c = vItCurve->second;
1099 
1100  if (vEdges.find(c.entityid) != vEdges.end() ||
1101  vFaces.find(c.entityid) != vFaces.end())
1102  {
1103  x = new TiXmlElement(c.entitytype);
1104  x->SetAttribute("ID", c.id);
1105  if (c.entitytype == "E")
1106  {
1107  x->SetAttribute("EDGEID", c.entityid);
1108  }
1109  else
1110  {
1111  x->SetAttribute("FACEID", c.entityid);
1112  }
1113  x->SetAttribute("TYPE", c.type);
1114  x->SetAttribute("NUMPOINTS", c.npoints);
1115  y = new TiXmlText(c.data);
1116  x->LinkEndChild(y);
1117  vCurved->LinkEndChild(x);
1118  }
1119  }
1120  }
1121 
1122  // Generate composites section comprising only those mesh entities
1123  // which belong to this partition.
1124  for (vIt = m_meshComposites.begin(); vIt != m_meshComposites.end(); ++vIt)
1125  {
1126  bool comma = false; // Set to true after first entity output
1127  bool range = false; // True when entity IDs form a range
1128  int last_idx = -2; // Last entity ID output
1129  std::string vCompositeStr = "";
1130  for (unsigned int j = 0; j < vIt->second.list.size(); ++j)
1131  {
1132  // Based on entity type, check if in this partition
1133  switch (vIt->second.type)
1134  {
1135  case 'V':
1136  if (vVertices.find(vIt->second.list[j]) == vVertices.end())
1137  {
1138  continue;
1139  }
1140  break;
1141  case 'E':
1142  if (vEdges.find(vIt->second.list[j]) == vEdges.end())
1143  {
1144  continue;
1145  }
1146  break;
1147  case 'F':
1148  if (vFaces.find(vIt->second.list[j]) == vFaces.end())
1149  {
1150  continue;
1151  }
1152  break;
1153  default:
1154  if (vElements.find(vIt->second.list[j]) == vElements.end())
1155  {
1156  continue;
1157  }
1158  break;
1159  }
1160 
1161  // Condense consecutive entity IDs into ranges
1162  // last_idx initially -2 to avoid error for ID=0
1163  if (last_idx + 1 == vIt->second.list[j])
1164  {
1165  last_idx++;
1166  range = true;
1167  continue;
1168  }
1169  // This entity is not in range, so close previous range with
1170  // last_idx
1171  if (range)
1172  {
1173  vCompositeStr += "-" + boost::lexical_cast<std::string>(last_idx);
1174  range = false;
1175  }
1176  // Output ID, which is either standalone, or will start a
1177  // range.
1178  vCompositeStr += comma ? "," : "";
1179  vCompositeStr += boost::lexical_cast<std::string>(vIt->second.list[j]);
1180  last_idx = vIt->second.list[j];
1181  comma = true;
1182  }
1183  // If last entity is part of a range, it must be output now
1184  if (range)
1185  {
1186  vCompositeStr += "-" + boost::lexical_cast<std::string>(last_idx);
1187  }
1188 
1189  if (vCompositeStr.length() > 0)
1190  {
1191  vComposites[vIt->first] = vIt->second;
1192  x = new TiXmlElement("C");
1193  x->SetAttribute("ID", vIt->first);
1194  vCompositeStr = "X[" + vCompositeStr + "]";
1195  vCompositeStr[0] = vIt->second.type;
1196  y = new TiXmlText(vCompositeStr.c_str());
1197  x->LinkEndChild(y);
1198  vComposite->LinkEndChild(x);
1199  }
1200  }
1201 
1202  std::string vDomainListStr;
1203  bool comma = false;
1204  for (unsigned int i = 0; i < m_domain.size(); ++i)
1205  {
1206  if (vComposites.find(m_domain[i]) != vComposites.end())
1207  {
1208  vDomainListStr += comma ? "," : "";
1209  comma = true;
1210  vDomainListStr += boost::lexical_cast<std::string>(m_domain[i]);
1211  }
1212  }
1213  vDomainListStr = "C[" + vDomainListStr + "]";
1214  TiXmlText* vDomainList = new TiXmlText(vDomainListStr);
1215  vDomain->LinkEndChild(vDomainList);
1216 
1217  vElmtGeometry->LinkEndChild(vVertex);
1218  if (m_dim >= 2)
1219  {
1220  vElmtGeometry->LinkEndChild(vEdge);
1221  }
1222  if (m_dim >= 3)
1223  {
1224  vElmtGeometry->LinkEndChild(vFace);
1225  }
1226  vElmtGeometry->LinkEndChild(vElement);
1227  if (m_dim >= 2)
1228  {
1229  vElmtGeometry->LinkEndChild(vCurved);
1230  }
1231  vElmtGeometry->LinkEndChild(vComposite);
1232  vElmtGeometry->LinkEndChild(vDomain);
1233 
1234  pNektar->LinkEndChild(vElmtGeometry);
1235 
1236  if (pSession->DefinesElement("Nektar/Conditions"))
1237  {
1238  std::set<int> vBndRegionIdList;
1239  TiXmlElement* vConditions = new TiXmlElement(*pSession->GetElement("Nektar/Conditions"));
1240  TiXmlElement* vBndRegions = vConditions->FirstChildElement("BOUNDARYREGIONS");
1241  TiXmlElement* vBndConditions = vConditions->FirstChildElement("BOUNDARYCONDITIONS");
1242  TiXmlElement* vItem;
1243 
1244  if (vBndRegions)
1245  {
1246  TiXmlElement* vNewBndRegions = new TiXmlElement("BOUNDARYREGIONS");
1247  vItem = vBndRegions->FirstChildElement();
1248  while (vItem)
1249  {
1250  std::string vSeqStr = vItem->FirstChild()->ToText()->Value();
1251  std::string::size_type indxBeg = vSeqStr.find_first_of('[') + 1;
1252  std::string::size_type indxEnd = vSeqStr.find_last_of(']') - 1;
1253  vSeqStr = vSeqStr.substr(indxBeg, indxEnd - indxBeg + 1);
1254  std::vector<unsigned int> vSeq;
1255  ParseUtils::GenerateSeqVector(vSeqStr.c_str(), vSeq);
1256  std::string vListStr;
1257  bool comma = false;
1258  for (unsigned int i = 0; i < vSeq.size(); ++i)
1259  {
1260  if (vComposites.find(vSeq[i]) != vComposites.end())
1261  {
1262  vListStr += comma ? "," : "";
1263  comma = true;
1264  vListStr += boost::lexical_cast<std::string>(vSeq[i]);
1265  }
1266  }
1267  int p = atoi(vItem->Attribute("ID"));
1268  if (vListStr.length() == 0)
1269  {
1270  vBndRegions->RemoveChild(vItem);
1271  }
1272  else
1273  {
1274  vListStr = "C[" + vListStr + "]";
1275  TiXmlText* vList = new TiXmlText(vListStr);
1276  TiXmlElement* vNewElement = new TiXmlElement("B");
1277  vNewElement->SetAttribute("ID", p);
1278  vNewElement->LinkEndChild(vList);
1279  vNewBndRegions->LinkEndChild(vNewElement);
1280  vBndRegionIdList.insert(p);
1281  }
1282 
1283  // Store original order of boundary region.
1284  m_bndRegOrder[p] = vSeq;
1285 
1286  vItem = vItem->NextSiblingElement();
1287  }
1288  vConditions->ReplaceChild(vBndRegions, *vNewBndRegions);
1289  }
1290 
1291  if (vBndConditions)
1292  {
1293  vItem = vBndConditions->FirstChildElement();
1294  while (vItem)
1295  {
1296  std::set<int>::iterator x;
1297  if ((x = vBndRegionIdList.find(atoi(vItem->Attribute("REF")))) != vBndRegionIdList.end())
1298  {
1299  vItem->SetAttribute("REF", *x);
1300  }
1301  else
1302  {
1303  vBndConditions->RemoveChild(vItem);
1304  }
1305  vItem = vItem->NextSiblingElement();
1306  }
1307  }
1308  pNektar->LinkEndChild(vConditions);
1309  }
1310 
1311  // Distribute other sections of the XML to each process as is.
1312  TiXmlElement* vSrc = pSession->GetElement("Nektar")
1313  ->FirstChildElement();
1314  while (vSrc)
1315  {
1316  std::string vName = boost::to_upper_copy(vSrc->ValueStr());
1317  if (vName != "GEOMETRY" && vName != "CONDITIONS")
1318  {
1319  pNektar->LinkEndChild(new TiXmlElement(*vSrc));
1320  }
1321  vSrc = vSrc->NextSiblingElement();
1322  }
1323  }
1324 
1325  void MeshPartition::GetElementIDs(const int procid, std::vector<unsigned int> &elmtid)
1326  {
1327  BoostVertexIterator vertit, vertit_end;
1328 
1329  ASSERTL0(procid < m_localPartition.size(),"procid is less than the number of partitions");
1330 
1331  // Populate lists of elements, edges and vertices required.
1332  for ( boost::tie(vertit, vertit_end) = boost::vertices(m_localPartition[procid]);
1333  vertit != vertit_end;
1334  ++vertit)
1335  {
1336  elmtid.push_back(m_meshElements[m_localPartition[procid][*vertit].id].id);
1337  }
1338  }
1339  }
1340 }
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2