Nektar++
SubStructuredGraph.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File SubStructuredGraph.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: a collection of classes that facilitates the implementation
33 // of the multi-level static condensation routines
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
40 
41 #include <iostream>
42 #include <algorithm>
43 #include <boost/config.hpp>
44 #include <boost/graph/adjacency_list.hpp>
45 #include <boost/graph/cuthill_mckee_ordering.hpp>
46 #include <boost/graph/properties.hpp>
47 #include <boost/graph/bandwidth.hpp>
48 
49 using std::max;
50 using std::cout;
51 using std::endl;
52 
53 namespace Nektar
54 {
55  namespace MultiRegions
56  {
58  {
59 
60  }
61 
62  PatchMap::PatchMap(const int nvals)
63  {
64  m_patchId = Array<OneD, int> (nvals);
65  m_dofId = Array<OneD, int> (nvals);
68  }
69 
71  {
72 
73  }
74 
76  const int n,
77  const int patchId,
78  const int dofId,
79  const unsigned int bndPatch,
80  const NekDouble sign)
81  {
82  m_patchId [n] = patchId;
83  m_dofId [n] = dofId;
84  m_bndPatch[n] = bndPatch;
85  m_sign [n] = sign;
86  }
87 
89  const Array<OneD, const int> sepTree) :
90  m_BndDofs (),
91  m_leftDaughterGraph (),
92  m_rightDaughterGraph()
93  {
94  static int offset = -5;
95  offset += 5;
96 
97  int recurLevel = sepTree[offset+0];
98  int nLeftIntDofs = sepTree[offset+2];
99  int nRightIntDofs = sepTree[offset+3];
100  int nBndDofs = sepTree[offset+4];
101 
102  bool daughtersConstructed[2] = {false,false};
103 
104  if ((offset + 5) < sepTree.num_elements())
105  {
106  while (sepTree[offset+5] > recurLevel)
107  {
108  switch (sepTree[offset+6])
109  {
110  case 1:
111  {
113  MultiLevelBisectedGraph>::AllocateSharedPtr(
114  sepTree);
115  daughtersConstructed[0] = true;
116  break;
117  }
118  case 2:
119  {
121  MultiLevelBisectedGraph>::AllocateSharedPtr(
122  sepTree);
123  daughtersConstructed[1] = true;
124  break;
125  }
126  default:
127  {
128  NEKERROR(ErrorUtil::efatal,"Invalid branch id");
129  }
130  }
131  if ((offset + 5) >= sepTree.num_elements())
132  {
133  break;
134  }
135  }
136  }
137 
139 
140  if (!daughtersConstructed[0] && nLeftIntDofs)
141  {
143  MultiLevelBisectedGraph>::AllocateSharedPtr(nLeftIntDofs);
144  }
145 
146  if (!daughtersConstructed[1] && nRightIntDofs)
147  {
149  MultiLevelBisectedGraph>::AllocateSharedPtr(nRightIntDofs);
150  }
151 
152  if (recurLevel == 1)
153  {
154  offset = -5;
155  }
156  }
157 
160  const int nPartition)
161  {
162  m_leftDaughterGraph = oldLevel;
164  }
165 
167  m_BndDofs(MemoryManager<SubGraph>::AllocateSharedPtr(nBndDofs)),
168  m_leftDaughterGraph(),
169  m_rightDaughterGraph()
170  {
171  }
172 
174  {
175  }
176 
178  {
179  static int nBndDofs = 0;
180  static int level = 0;
181  level++;
182 
183  int returnval;
184 
185  if(m_leftDaughterGraph.get())
186  {
187  m_leftDaughterGraph->GetTotDofs();
188  }
189  if(m_rightDaughterGraph.get())
190  {
191  m_rightDaughterGraph->GetTotDofs();
192  }
193 
194  nBndDofs += m_BndDofs->GetNverts();
195  returnval = nBndDofs;
196 
197  level--;
198  if(level == 0)
199  {
200  nBndDofs = 0;
201  }
202 
203  return returnval;
204  }
205 
207  {
208  static int level = 0;
209  static int offset = 0;
210  level++;
211 
212  if(m_leftDaughterGraph.get())
213  {
214  m_leftDaughterGraph->SetGlobalNumberingOffset();
215  }
216  if(m_rightDaughterGraph.get())
217  {
218  m_rightDaughterGraph->SetGlobalNumberingOffset();
219  }
220 
221  m_BndDofs->SetIdOffset(offset);
222  offset += m_BndDofs->GetNverts();
223 
224  level--;
225  if(level == 0)
226  {
227  offset = 0;
228  }
229  }
230 
232  {
233  static int level = 0;
234  level++;
235  cout << "LEVEL " << level << " " << m_BndDofs->GetNverts() << endl;
236 
237  if (m_leftDaughterGraph.get())
238  {
239  m_leftDaughterGraph->DumpNBndDofs();
240  }
241  if (m_rightDaughterGraph.get())
242  {
243  m_rightDaughterGraph->DumpNBndDofs();
244  }
245 
246  level--;
247  }
248 
250  std::vector<SubGraphSharedPtr>& leaves) const
251  {
252  int cnt = 0;
253 
254  if (m_leftDaughterGraph.get())
255  {
256  m_leftDaughterGraph->CollectLeaves(leaves);
257  cnt++;
258  }
259 
260  if (m_rightDaughterGraph.get())
261  {
262  m_rightDaughterGraph->CollectLeaves(leaves);
263  cnt++;
264  }
265 
266  if (cnt == 0)
267  {
269  leaves.push_back(leave);
270  }
271  }
272 
274  {
275  int cnt = 0;
276 
277  if (m_leftDaughterGraph.get())
278  {
279  cnt++;
280  }
281 
282  if (m_rightDaughterGraph.get())
283  {
284  cnt++;
285  }
286 
287  return cnt;
288  }
289 
291  {
292  int returnval;
293  static int level = 0;
294  static int nLeaves = 0;
295  level++;
296 
297  if (level == 1 &&
298  !m_leftDaughterGraph.get() && !m_rightDaughterGraph.get())
299  {
300  level = 0;
301  nLeaves = 0;
302  return 0;
303  }
304 
305  if (m_leftDaughterGraph.get())
306  {
307  if (m_leftDaughterGraph->GetNdaughterGraphs() == 0 &&
308  m_leftDaughterGraph->GetBndDofsGraph()->GetNverts() == 0)
309  {
311  nLeaves++;
312  }
313  else
314  {
315  m_leftDaughterGraph->CutEmptyLeaves();
316  }
317  }
318 
319  if(m_rightDaughterGraph.get())
320  {
321  if (m_rightDaughterGraph->GetNdaughterGraphs() == 0 &&
322  m_rightDaughterGraph->GetBndDofsGraph()->GetNverts() == 0)
323  {
325  nLeaves++;
326  }
327  else
328  {
329  m_rightDaughterGraph->CutEmptyLeaves();
330  }
331  }
332 
333  returnval = nLeaves;
334 
335  level--;
336  if(level == 0)
337  {
338  nLeaves = 0;
339  }
340 
341  return returnval;
342  }
343 
345  {
346  int returnval;
347  static int level = 0;
348  static int nLeaves = 0;
349  level++;
350 
351  if( (level == 1) &&
352  (!m_leftDaughterGraph.get()) &&
353  (!m_rightDaughterGraph.get()) )
354  {
355  level = 0;
356  nLeaves = 0;
357  return 0;
358  }
359 
360  if(m_leftDaughterGraph.get())
361  {
362  if(m_leftDaughterGraph->GetNdaughterGraphs() == 0)
363  {
365  nLeaves++;
366  }
367  else
368  {
369  m_leftDaughterGraph->CutLeaves();
370  }
371  }
372  if(m_rightDaughterGraph.get())
373  {
374  if(m_rightDaughterGraph->GetNdaughterGraphs() == 0)
375  {
377  nLeaves++;
378  }
379  else
380  {
381  m_rightDaughterGraph->CutLeaves();
382  }
383  }
384 
385  returnval = nLeaves;
386 
387  level--;
388  if(level == 0)
389  {
390  nLeaves = 0;
391  }
392 
393  return returnval;
394  }
395 
396 
398  const int nVerts) : m_IntBlocks(), m_daughterGraph()
399  {
400  // This constructor should only be used in the very special case
401  // when there is a graph consisting of nVerts vertices but without
402  // any connectivity whatsowever (i.e. no edges)
403 
404  SubGraphSharedPtr subgraph;
405  for (int i = 0 ; i < nVerts; i++)
406  {
408  m_IntBlocks.push_back(subgraph);
409  }
410  }
411 
413  const Array<OneD, const int> septree,
414  const int nPartition) :
415  m_IntBlocks(),
416  m_daughterGraph()
417  {
418  // First, create a top-down graph structure based upon the separator
419  // tree. This is easier as separation tree is also structured
420  // following a top-down approach
421  MultiLevelBisectedGraphSharedPtr topDownGraph =
423  septree);
424 
425  if (nPartition > 0)
426  {
428  AllocateSharedPtr(topDownGraph, nPartition);
429  }
430 
431  // set the global numbering of the top-down graph
432  topDownGraph->SetGlobalNumberingOffset();
433 
434  topDownGraph->CutEmptyLeaves();
435 
436  // Secondly, recursively construct the subgraphs of the bottom up
437  // point of view 1. Collect all the leaves of the topdown graph this
438  // will be the first level of the bottom up graph
439  topDownGraph->CollectLeaves(m_IntBlocks);
440 
441  // 2. Reduce the topdown graph by cutting the leaves (this will
442  // allow a recursive approach)
443  int ncuts = topDownGraph->CutLeaves();
444 
445  // 3. If there were leaves to cut, proceed recursively
446  if (ncuts)
447  {
449  AllocateSharedPtr(topDownGraph);
450  }
451  }
452 
454  const MultiLevelBisectedGraphSharedPtr& graph) :
455  m_IntBlocks (),
456  m_daughterGraph()
457  {
458  int ncuts;
459  graph->CutEmptyLeaves();
460  graph->CollectLeaves(m_IntBlocks);
461  ncuts = graph->CutLeaves();
462 
463  if(ncuts)
464  {
466  AllocateSharedPtr(graph);
467  }
468  }
469 
471  {
472  }
473 
475  {
476  static int nIntDofs = 0;
477  static int level = 0;
478  level++;
479 
480  int returnval;
481 
482  if( m_daughterGraph.get())
483  {
484  m_daughterGraph->GetTotDofs();
485  }
486 
487  for(int i = 0; i < m_IntBlocks.size(); i++)
488  {
489  nIntDofs += m_IntBlocks[i]->GetNverts();
490  }
491 
492  returnval = nIntDofs;
493 
494  level--;
495  if(level == 0)
496  {
497  nIntDofs = 0;
498  }
499 
500  return returnval;
501  }
502 
504  {
505  static int level = 0;
506  level++;
507 
508  cout << "LEVEL " << level << endl;
509  cout << "interior blocks" << endl;
510  for (int i = 0; i < m_IntBlocks.size(); i++)
511  {
512  cout << " " << i
513  << "/" << m_IntBlocks.size()-1
514  << ": " << m_IntBlocks[i]->GetNverts()
515  << ", " << m_IntBlocks[i]->GetIdOffset() << endl;
516  }
517  cout << endl;
518 
519  if (m_daughterGraph.get())
520  {
521  m_daughterGraph->Dump();
522  }
523 
524  level--;
525  }
526 
528  Array<OneD, int> &perm,
529  Array<OneD, int> &iperm) const
530  {
531  int nDofs = GetTotDofs();
532 
533  // Step 1: make a permutation array that goes from the current
534  // reordering in the bottom-up graph to an ordering where the
535  // interior dofs of the first (=bottom) level are ordered last,
536  // followed interior dofs of the second level, ...
537  Array<OneD, int> iperm1(nDofs);
538  SetBottomUpReordering(iperm1);
539 
540  // Now, based upon the input permutation array, update the
541  // permutation arrays between the original ordering of the dofs
542  // (defined somewhere outside) and the final reordering as given by
543  // iperm1
544  for (int i=0; i < nDofs; i++)
545  {
546  iperm[i] = iperm1[iperm[i]];
547  perm[iperm[i]] = i;
548  }
549  }
550 
552  Array<OneD, int>& iperm) const
553  {
554  static int offset = 0;
555  static int level = 0;
556  level++;
557 
558  if( m_daughterGraph.get() )
559  {
560  m_daughterGraph->SetBottomUpReordering(iperm);
561  }
562 
563  for(int i = 0; i < m_IntBlocks.size(); i++)
564  {
565  int GlobIdOffset = m_IntBlocks[i]->GetIdOffset();
566  m_IntBlocks[i]->SetIdOffset(offset);
567  for(int j = 0; j < m_IntBlocks[i]->GetNverts(); j++)
568  {
569  iperm[GlobIdOffset+j] = offset;
570  offset++;
571  }
572  }
573 
574  level--;
575  if(level == 0)
576  {
577  offset = 0;
578  }
579  }
580 
582  const Array<OneD, const int>& wgts)
583  {
584  static int offset = 0;
585  static int level = 0;
586  level++;
587 
588  if( m_daughterGraph.get())
589  {
590  m_daughterGraph->ExpandGraphWithVertexWeights(wgts);
591  }
592 
593  for(int i = 0; i < m_IntBlocks.size(); i++)
594  {
595  int OrigGlobIdOffset = m_IntBlocks[i]->GetIdOffset();
596  int newNverts = 0;
597 
598  m_IntBlocks[i]->SetIdOffset(offset);
599 
600  for(int j = 0; j < m_IntBlocks[i]->GetNverts(); j++)
601  {
602  newNverts += wgts[OrigGlobIdOffset+j];
603  offset += wgts[OrigGlobIdOffset+j];
604  }
605 
606  m_IntBlocks[i]->SetNverts(newNverts);
607  }
608 
609  // erase the blocks that do not have any vertices
610  m_IntBlocks.erase(
611  remove_if(m_IntBlocks.begin(),
612  m_IntBlocks.end(),
614  m_IntBlocks.end());
615 
616  // remove the current level if there are no interior blocks
617  if(m_IntBlocks.size() == 0 && m_daughterGraph.get())
618  {
619  m_IntBlocks = m_daughterGraph->GetInteriorBlocks();
620  m_daughterGraph = m_daughterGraph->GetDaughterGraph();
621  }
622 
623  level--;
624  if(level == 0)
625  {
626  offset = 0;
627  }
628  }
629 
631  const int leveltomask,
632  Array<OneD, NekDouble> &maskarray) const
633  {
634  static int level = 0;
635  level++;
636 
637  if(level < leveltomask)
638  {
639  m_daughterGraph->MaskPatches(leveltomask,maskarray);
640  }
641  else if(level == leveltomask)
642  {
643  int GlobIdOffset;
644  int nVerts;
645  for(int i = 0; i < m_IntBlocks.size(); i++)
646  {
647  GlobIdOffset = m_IntBlocks[i]->GetIdOffset();
648  nVerts = m_IntBlocks[i]->GetNverts();
649  for(int j = 0; j < nVerts; j++)
650  {
651  maskarray[GlobIdOffset+j] = (NekDouble) i;
652  }
653  }
654  }
655  else
656  {
658  "If statement should not arrive here");
659  }
660 
661  level--;
662  }
663 
665  const int whichlevel) const
666  {
667  int returnval = -1;
668  static int level = 0;
669  level++;
670 
671  if(level < whichlevel)
672  {
673  returnval = m_daughterGraph->GetNpatchesWithInterior(
674  whichlevel);
675  }
676  else if(level == whichlevel)
677  {
678  returnval = m_IntBlocks.size();
679  }
680  else
681  {
683  "If statement should not arrive here");
684  }
685 
686  level--;
687 
688  return returnval;
689  }
690 
692  const int whichlevel,
693  Array<OneD, unsigned int> &outarray) const
694  {
695  static int level = 0;
696  level++;
697 
698  if(level < whichlevel)
699  {
700  m_daughterGraph->GetNintDofsPerPatch(whichlevel,outarray);
701  }
702  else if(level == whichlevel)
703  {
704  ASSERTL1(outarray.num_elements() >= m_IntBlocks.size(),
705  "Array dimension not sufficient");
706 
707  for(int i = 0; i < m_IntBlocks.size(); i++)
708  {
709  outarray[i] = (unsigned int) m_IntBlocks[i]->GetNverts();
710  }
711  }
712  else
713  {
715  "If statement should not arrive here");
716  }
717 
718  level--;
719  }
720 
722  const int whichlevel, const int patch) const
723  {
724  int retval = -1;
725  static int level = 0;
726  level++;
727 
728  if(level < whichlevel)
729  {
730  retval = m_daughterGraph->GetInteriorOffset(whichlevel,patch);
731  }
732  else if(level == whichlevel)
733  {
734  retval = m_IntBlocks[patch]->GetIdOffset();
735  }
736  else
737  {
739  "If statement should not arrive here");
740  }
741 
742  level--;
743 
744  return retval;
745  }
746 
747  std::vector<SubGraphSharedPtr> BottomUpSubStructuredGraph::
748  GetInteriorBlocks(const int whichlevel) const
749  {
750  std::vector<SubGraphSharedPtr> returnval;
751  static int level = 0;
752  level++;
753 
754  if(level < whichlevel)
755  {
756  returnval = m_daughterGraph->GetInteriorBlocks(whichlevel);
757  }
758  else if(level == whichlevel)
759  {
760  returnval = m_IntBlocks;
761  }
762  else
763  {
765  "If statement should not arrive here");
766  }
767 
768  level--;
769  return returnval;
770  }
771 
773  const int whichlevel) const
774  {
775  int returnval = -1;
776  static int level = 0;
777  level++;
778 
779  if(level < whichlevel)
780  {
781  returnval = m_daughterGraph->GetNumGlobalDofs(whichlevel);
782  }
783  else if(level == whichlevel)
784  {
785  returnval = m_IntBlocks[m_IntBlocks.size()-1]->GetIdOffset()+
786  m_IntBlocks[m_IntBlocks.size()-1]->GetNverts();
787  }
788  else
789  {
791  "If statement should not arrive here");
792  }
793 
794  level--;
795 
796  return returnval;
797  }
798 
800  {
801  int returnval = 0;
802  static int level = 0;
803  level++;
804 
805  if( m_daughterGraph.get())
806  {
807  returnval = m_daughterGraph->GetNlevels();
808  }
809  returnval = max(returnval,level);
810 
811  level--;
812 
813  return returnval;
814 
815  }
816 
818  {
819  return g->GetNverts() == 0;
820  }
821 
822  namespace
823  {
824  // the first template parameter (=OutEdgeList) is chosen to be of
825  // type std::set as in the set up of the adjacency, a similar edge
826  // might be created multiple times. And to prevent the definition
827  // of parallell edges, we use std::set (=boost::setS) rather than
828  // std::vector (=boost::vecS)
829  typedef boost::adjacency_list<
830  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
831  typedef boost::graph_traits<BoostGraph>::vertex_descriptor
832  BoostVertex;
833  typedef boost::graph_traits<BoostGraph>::vertex_iterator
834  BoostVertexIterator;
835  typedef boost::graph_traits<BoostGraph>::adjacency_iterator
836  BoostAdjacencyIterator;
837  }
838 
839  void CuthillMckeeReordering(const BoostGraph& graph,
840  Array<OneD, int>& perm,
841  Array<OneD, int>& iperm)
842  {
843  int nGraphVerts = boost::num_vertices(graph);
844 
845  ASSERTL1(perm. num_elements() >= nGraphVerts &&
846  iperm.num_elements() >= nGraphVerts,
847  "Non-matching dimensions");
848 
849  // Call boost::cuthill_mckee_ordering to reorder the graph-vertices
850  // using the reverse Cuthill-Mckee algorithm
851  std::vector<BoostVertex> reorderedVerts(nGraphVerts);
852  boost::cuthill_mckee_ordering(graph, reorderedVerts.rbegin());
853 
854  //copy the reordering to the Arrays perm and iperm
855  for(int i = 0; i < nGraphVerts; i++)
856  {
857  perm[i] = reorderedVerts[i];
858  iperm[ reorderedVerts[i] ] = i;
859  }
860  }
861 
863  const BoostGraph &graph,
864  Array<OneD, int> &perm,
865  Array<OneD, int> &iperm,
866  BottomUpSubStructuredGraphSharedPtr &substructgraph,
867  std::set<int> partVerts,
868  int mdswitch)
869  {
870  int nGraphVerts = boost::num_vertices(graph);
871  int nGraphEdges = boost::num_edges (graph);
872 
873  ASSERTL1(perm. num_elements() >= nGraphVerts &&
874  iperm.num_elements() >= nGraphVerts,
875  "Non-matching dimensions");
876 
877  // We will now use METIS to reorder the graph. For the purpose of
878  // multi-level static condensation, we will use a METIS routine that
879  // partitions the graph recursively using a multi-level bisection
880  // algorithm. The name of this routine is METIS_NodeND and it was
881  // originally designed to reorder the DOFs in a matrix in order to
882  // minimise the fill-in when applying a factorisation technique
883  // (such as Cholesky). However, this reordering of DOFs also seems
884  // to be perfectly suited in the context of multilevel
885  // substructering. Therefore, we will use this metis routine instead
886  // of the more well-known graph-partitioning routines. However, as
887  // the standard metis implementation of METIS_NodeND only gives the
888  // resulting re-ordering as an output, we we will use an modified
889  // version of this routine that also returns information about the
890  // structure of the multi-level bisected partitioning.
891 
892  // This modified implementation has been written by W. GAO and
893  // collaborators and it additionally returns the separator tree
894  // compared to the standard implementation.
895 
896  // The name of this modified routine AS_METIS_NodeND (where AS
897  // stands for automated substructering) More information can be
898  // found in the paper:
899  //
900  // W. Gao, S. Li Xiaoye, C. Yang and Z. Bai
901  // 'An implementation and evaluation of the AMLS method
902  // for sparse eigenvalue problems'
903  // - ACM Trans. Math. Softw. 34, 4, Article 20 (July 2008)
904 
905  if(nGraphEdges)
906  {
907  // Step 1: Convert boost graph to a graph in adjncy-list format
908  // as required by METIS
909  int acnt = 0;
910  int vcnt = 0;
911  int i, cnt;
912  int nPartition = partVerts.size();
913  int nNonPartition = nGraphVerts - partVerts.size();
914  BoostVertexIterator vertit, vertit_end;
915  BoostAdjacencyIterator adjvertit, adjvertit_end;
916  Array<OneD, int> xadj(nNonPartition+1,0);
917  Array<OneD, int> adjncy(2*nGraphEdges);
918  Array<OneD, int> initial_perm(nGraphVerts);
919  Array<OneD, int> iinitial_perm(nGraphVerts);
920  Array<OneD, int> perm_tmp (nNonPartition);
921  Array<OneD, int> iperm_tmp(nNonPartition);
922 
924 
925  // First reorder vertices so that partition nodes are at the
926  // end. This allows METIS to partition the interior nodes.
927  for (i = cnt = 0; i < nGraphVerts; ++i)
928  {
929  if (partVerts.count(i) == 0)
930  {
931  initial_perm [i] = cnt;
932  iinitial_perm[cnt++] = i;
933  }
934  }
935 
936  for (i = 0; i < nGraphVerts; ++i)
937  {
938  if (partVerts.count(i) > 0)
939  {
940  initial_perm [i] = cnt;
941  iinitial_perm[cnt++] = i;
942  }
943  }
944 
945  // Apply this reordering to the graph.
946  boost::property_map<BoostGraph, boost::vertex_index_t>::type
947  index = get(boost::vertex_index, graph);
948 
949  for (boost::tie(vertit, vertit_end) = boost::vertices(graph);
950  vertit != vertit_end; ++vertit)
951  {
952  if (partVerts.count(index[*vertit]) > 0)
953  {
954  continue;
955  }
956 
957  for (boost::tie(adjvertit, adjvertit_end) =
958  boost::adjacent_vertices(*vertit,graph);
959  adjvertit != adjvertit_end;
960  ++adjvertit)
961  {
962  if (partVerts.count(index[*adjvertit]) > 0)
963  {
964  continue;
965  }
966  adjncy[acnt++] = initial_perm[*adjvertit];
967  }
968  xadj[++vcnt] = acnt;
969  }
970 
971  // Step 2: use metis to reorder the dofs. We do not know on
972  // forehand the size of the separator tree that METIS will
973  // return, so we just assume a really big value and try with
974  // that
975  int sizeSeparatorTree = nGraphVerts*10;
976  Array<OneD,int> septreeTmp(sizeSeparatorTree,-1);
977 
978  // The separator tree returned by metis has the following
979  // structure: It is a one dimensional array and information per
980  // level is contained per 5 elements:
981  //
982  // m_septree[i*5 + 0]: the level of recursion (top-level = 1)
983  // m_septree[i*5 + 1]: is this substructure a left or right
984  // branch? 1 = left branch, 2 = right branch
985  // m_septree[i*5 + 2]: the number of 'interior' DOFs in left
986  // branch
987  // m_septree[i*5 + 3]: the number of 'interior' DOFs in right
988  // branch
989  // m_septree[i*5 + 4]: the number of 'boundary' DOFs
990 
991  // Now try to call Call METIS.
992  try
993  {
995  nNonPartition,xadj,adjncy,perm_tmp,iperm_tmp,
996  septreeTmp, mdswitch);
997  }
998  catch(...)
999  {
1001  "Error in calling metis (the size of the separator"
1002  " tree might not be sufficient)");
1003  }
1004 
1005  // Change permutations from METIS to account for initial offset.
1006  for (i = 0; i < nGraphVerts; ++i)
1007  {
1008  if (partVerts.count(i) == 0)
1009  {
1010  iperm[i] = iperm_tmp[initial_perm[i]];
1011  perm[iperm[i]] = i;
1012  }
1013  }
1014 
1015  for (i = nNonPartition, it = partVerts.begin();
1016  it != partVerts.end(); ++it, ++i)
1017  {
1018  perm [i] = *it;
1019  iperm[*it] = i;
1020  }
1021 
1022  for (i = 0; i < nGraphVerts; ++i)
1023  {
1024  ASSERTL0(perm[iperm[i]] == i,
1025  "Perm error " + boost::lexical_cast<std::string>(i));
1026  }
1027 
1028  // Post-process the separator tree
1029  int trueSizeSepTree = 0;
1030  for (i = 0; septreeTmp[i] != -1; i++)
1031  {
1032  trueSizeSepTree++;
1033  }
1034  Array<OneD,int> septree(trueSizeSepTree);
1035  Vmath::Vcopy(trueSizeSepTree,septreeTmp,1,septree,1);
1036 
1037  // Based upon the separator tree, where are going to set up an
1038  // object of the class BottomUpSubStructuredGraph. The
1039  // constructor will read the separatortree and will interprete
1040  // the information from a bottom-up point of view.
1042  AllocateSharedPtr(septree, nPartition);
1043 
1044  // Important, we cannot simply use the ordering given by metis
1045  // as it does not order the different blocks as we would like
1046  // it. Therefore, we use following command to re-order them
1047  // again in the context of the bottom-up substructering. As a
1048  // result, we will now obtain an ordering where the interior
1049  // degrees of freedom of the first (=bottom) level will be
1050  // ordered last (block by block ofcoarse). The interior degrees
1051  // of freedom of the second level will be ordered second to
1052  // last, etc ... As a result, the boundary degrees of freedom of
1053  // the last level (i.e. the dofs that will have to solved
1054  // non-recursively) will be ordered first (after the Dirichlet
1055  // Dofs that is). (this way, we actually follow the same idea
1056  // and convention in the standard (non-multi-level) static
1057  // condensation approach)
1058  substructgraph->UpdateBottomUpReordering(perm,iperm);
1059  }
1060  else
1061  {
1062  // This is the very special case of a graph without any
1063  // connectivity i.e. a collection of vertices without any edges
1064  for(int i = 0; i < nGraphVerts; i++)
1065  {
1066  perm[i] = i;
1067  iperm[i] = i;
1068  }
1070  AllocateSharedPtr(nGraphVerts);
1071  }
1072  }
1073 
1074  void NoReordering(const BoostGraph& graph,
1075  Array<OneD, int>& perm,
1076  Array<OneD, int>& iperm)
1077  {
1078  int nGraphVerts = boost::num_vertices(graph);
1079 
1080  ASSERTL1(perm. num_elements() >= nGraphVerts &&
1081  iperm.num_elements() >= nGraphVerts,
1082  "Non-matching dimensions");
1083 
1084  for (int i = 0; i < nGraphVerts; i++)
1085  {
1086  perm [i] = i;
1087  iperm[i] = i;
1088  }
1089  }
1090  }
1091 }
Array< OneD, NekDouble > m_sign
int GetNpatchesWithInterior(const int whichlevel) const
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:158
void MaskPatches(const int leveltomask, Array< OneD, NekDouble > &maskarray) const
void GetNintDofsPerPatch(const int whichlevel, Array< OneD, unsigned int > &outarray) const
void MultiLevelBisectionReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm, BottomUpSubStructuredGraphSharedPtr &substructgraph, std::set< int > partVerts, int mdswitch)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:22
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
MultiLevelBisectedGraph(const Array< OneD, const int > sepTree)
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
void CuthillMckeeReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
BottomUpSubStructuredGraph(const Array< OneD, const int > septree, const int nPartition)
boost::shared_ptr< MultiLevelBisectedGraph > MultiLevelBisectedGraphSharedPtr
void CollectLeaves(std::vector< SubGraphSharedPtr > &leaves) const
boost::shared_ptr< SubGraph > SubGraphSharedPtr
static void as_onmetis(int nVerts, Nektar::Array< Nektar::OneD, int > xadj, Nektar::Array< Nektar::OneD, int > adjncy, Nektar::Array< Nektar::OneD, int > perm, Nektar::Array< Nektar::OneD, int > iperm, Nektar::Array< Nektar::OneD, int > map, int mdswitch=1)
Definition: Metis.hpp:53
Array< OneD, unsigned int > m_bndPatch
void ExpandGraphWithVertexWeights(const Array< OneD, const int > &wgts)
MultiLevelBisectedGraphSharedPtr m_rightDaughterGraph
double NekDouble
std::vector< SubGraphSharedPtr > GetInteriorBlocks() const
BottomUpSubStructuredGraphSharedPtr m_daughterGraph
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void UpdateBottomUpReordering(Array< OneD, int > &perm, Array< OneD, int > &iperm) const
bool SubGraphWithoutVerts(const SubGraphSharedPtr g)
MultiLevelBisectedGraphSharedPtr m_leftDaughterGraph
void SetBottomUpReordering(Array< OneD, int > &iperm) const
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
void SetPatchMap(const int n, const int patchId, const int dofId, const unsigned int bndPatch, const NekDouble sign)
int GetInteriorOffset(const int whichlevel, const int patch=0) const
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)