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