Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AssemblyMapCG3D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMapCG3D.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: C0-continuous assembly mappings specific to 3D
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/MultiRegions.hpp>
37 #include <MultiRegions/AssemblyMap/AssemblyMapCG3D.h>
38 #include <LocalRegions/Expansion2D.h>
39 #include <LocalRegions/Expansion3D.h>
40 
41 #include <boost/config.hpp>
42 #include <boost/graph/adjacency_list.hpp>
43 #include <boost/graph/cuthill_mckee_ordering.hpp>
44 #include <boost/graph/properties.hpp>
45 #include <boost/graph/bandwidth.hpp>
46 
47 namespace Nektar
48 {
49  namespace MultiRegions
50  {
51  /**
52  * @class AssemblyMapCG3D
53  * Mappings are created for three possible global solution types:
54  * - Direct full matrix
55  * - Direct static condensation
56  * - Direct multi-level static condensation
57  * In the latter case, mappings are created recursively for the
58  * different levels of static condensation.
59  *
60  * These mappings are used by GlobalLinSys to generate the global
61  * system.
62  */
63 
64  /**
65  *
66  */
67  AssemblyMapCG3D::AssemblyMapCG3D(
68  const LibUtilities::SessionReaderSharedPtr &pSession,
69  const std::string variable):
70  AssemblyMapCG(pSession,variable)
71  {
72  }
73 
74 
75 
76  /**
77  *
78  */
79  AssemblyMapCG3D::AssemblyMapCG3D(
80  const LibUtilities::SessionReaderSharedPtr &pSession,
81  const int numLocalCoeffs,
82  const ExpList &locExp,
83  const Array<OneD, const ExpListSharedPtr> &bndCondExp,
84  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
85  &bndConditions,
86  const PeriodicMap &periodicVerts,
87  const PeriodicMap &periodicEdges,
88  const PeriodicMap &periodicFaces,
89  const bool checkIfSystemSingular,
90  const std::string variable):
91  AssemblyMapCG(pSession,variable)
92  {
93  SetUp3DExpansionC0ContMap(numLocalCoeffs,
94  locExp,
95  bndCondExp,
96  bndConditions,
97  periodicVerts,
98  periodicEdges,
99  periodicFaces,
100  checkIfSystemSingular);
101  CalculateBndSystemBandWidth();
102  CalculateFullSystemBandWidth();
103  }
104 
105 
106  /**
107  *
108  */
109  AssemblyMapCG3D::AssemblyMapCG3D(
110  const LibUtilities::SessionReaderSharedPtr &pSession,
111  const int numLocalCoeffs,
112  const ExpList &locExp):
113  AssemblyMapCG(pSession)
114  {
115  SetUp3DExpansionC0ContMap(numLocalCoeffs, locExp);
116  CalculateBndSystemBandWidth();
117  CalculateFullSystemBandWidth();
118  }
119 
120 
121  /**
122  *
123  */
124  AssemblyMapCG3D::~AssemblyMapCG3D()
125  {
126  }
127 
128 
129 
130 
131  /**
132  * Construction of the local->global map is achieved in several stages.
133  * A mesh vertex, mesh edge and mesh face renumbering is constructed
134  * in #vertReorderedGraphVertId, #edgeReorderedGraphVertId and
135  * #faceReorderedGraphVertId
136  *
137  * The only unique identifiers of the vertices, edges and faces of the
138  * mesh are the vertex id and the mesh id (stored in their corresponding
139  * Geometry object). However, setting up a global numbering based on
140  * these id's will not lead to a suitable or optimal numbering. Mainly
141  * because:
142  * - we want the Dirichlet DOF's to be listed first
143  * - we want an optimal global numbering of the remaining DOF's
144  * (strategy still need to be defined but can for example be:
145  * minimum bandwith or minimum fill-in of the resulting global
146  * system matrix)
147  *
148  * That's why the vertices, edges and faces will be rearranged. This is
149  * done in the following way: The vertices, edges and faces of the mesh
150  * are considered as vertices of a graph (in a computer science way)
151  * (equivalently, they can also be considered as boundary degrees of
152  * freedom, whereby all boundary modes of a single edge are considered
153  * as a single DOF). We then will use algorithms to reorder these
154  * graph-vertices (or boundary DOF's).
155  *
156  * We will use a boost graph object to store this graph the first
157  * template parameter (=OutEdgeList) is chosen to be of type std::set
158  * as in the set up of the adjacency, a similar edge might be created
159  * multiple times. And to prevent the definition of parallel edges,
160  * we use std::set (=boost::setS) rather than std::vector
161  * (=boost::vecS).
162  *
163  * Two different containers are used to store the graph vertex id's of
164  * the different mesh vertices and edges. They are implemented as a STL
165  * map such that the graph vertex id can later be retrieved by the
166  * unique mesh vertex or edge id's which serve as the key of the map.
167  *
168  * Therefore, the algorithm proceeds as follows:
169  */
170  void AssemblyMapCG3D::SetUp3DExpansionC0ContMap(
171  const int numLocalCoeffs,
172  const ExpList &locExp,
173  const Array<OneD, const ExpListSharedPtr> &bndCondExp,
174  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &bndConditions,
175  const PeriodicMap &periodicVerts,
176  const PeriodicMap &periodicEdges,
177  const PeriodicMap &periodicFaces,
178  const bool checkIfSystemSingular)
179  {
180  int i,j,k,l;
181  int cnt = 0;
182  int intDofCnt;
183  int meshVertId;
184  int meshVertId2;
185  int meshEdgeId;
186  int meshEdgeId2;
187  int meshFaceId;
188  int globalId;
189  int nEdgeInteriorCoeffs;
190  int nFaceInteriorCoeffs;
191  int firstNonDirGraphVertId;
192  int nLocBndCondDofs = 0;
193  int nLocDirBndCondDofs = 0;
194  int graphVertId = 0;
195  LocalRegions::Expansion3DSharedPtr locExpansion;
196  LocalRegions::Expansion2DSharedPtr bndCondFaceExp;
197  LibUtilities::BasisType bType;
198  StdRegions::Orientation edgeOrient;
199  StdRegions::Orientation faceOrient;
200  Array<OneD, unsigned int> edgeInteriorMap;
201  Array<OneD, int> edgeInteriorSign;
202  Array<OneD, unsigned int> faceInteriorMap;
203  Array<OneD, int> faceInteriorSign;
204  PeriodicMap::const_iterator pIt;
205 
206  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
207 
208  m_signChange = false;
209  //m_systemSingular = false;
210 
211  map<int,int> vertReorderedGraphVertId;
212  map<int,int> edgeReorderedGraphVertId;
213  map<int,int> faceReorderedGraphVertId;
214  map<int,int>::iterator mapIt;
215  map<int,int>::const_iterator mapConstIt;
216  map<int,pair<int, StdRegions::Orientation> >::const_iterator mapFaceIt;
217 
218  bool systemSingular = (checkIfSystemSingular)? true: false;
219 
220  /**
221  * STEP 1: Order the Dirichlet vertices and edges first
222  */
223  for(i = 0; i < bndCondExp.num_elements(); i++)
224  {
225  // Check to see if any value on face has Dirichlet value.
226  cnt = 0;
227  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
228  {
229  bndCondFaceExp = bndCondExp[i]->GetExp(j)
230  ->as<LocalRegions::Expansion2D>();
231 
232  if (bndConditions[i]->GetBoundaryConditionType() ==
233  SpatialDomains::eDirichlet)
234  {
235  for(k = 0; k < bndCondFaceExp->GetNverts(); k++)
236  {
237  meshVertId = bndCondFaceExp->GetGeom2D()->GetVid(k);
238  if(vertReorderedGraphVertId.count(meshVertId) == 0)
239  {
240  vertReorderedGraphVertId[meshVertId] = graphVertId++;
241  }
242  }
243 
244  for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
245  {
246  meshEdgeId = bndCondFaceExp->GetGeom2D()->GetEid(k);
247  if(edgeReorderedGraphVertId.count(meshEdgeId) == 0)
248  {
249  edgeReorderedGraphVertId[meshEdgeId] = graphVertId++;
250  }
251  }
252  meshFaceId = bndCondFaceExp->GetGeom2D()->GetFid();
253  faceReorderedGraphVertId[meshFaceId] = graphVertId++;
254  nLocDirBndCondDofs += bndCondFaceExp->GetNcoeffs();
255  }
256  if (bndConditions[i]->GetBoundaryConditionType() !=
257  SpatialDomains::eNeumann &&
258  checkIfSystemSingular == true)
259  {
260  systemSingular = false;
261  }
262  nLocBndCondDofs += bndCondFaceExp->GetNcoeffs();
263  }
264  }
265 
266  // Number of dirichlet edges and faces (not considering periodic BCs)
267  m_numDirEdges = edgeReorderedGraphVertId.size();
268  m_numDirFaces = faceReorderedGraphVertId.size();
269 
270  /**
271  * STEP 1.5: Exchange Dirichlet mesh vertices between processes and
272  * check for singular problems.
273  */
274 
275  // Collate information on Dirichlet vertices from all processes
276  int n = m_comm->GetSize();
277  int p = m_comm->GetRank();
278 
279  // At this point vertReorderedGraphVertId only contains
280  // information from Dirichlet boundaries. Therefore make a
281  // global list of the vert and edge information on all
282  // processors
283  Array<OneD, int> vertcounts (n, 0);
284  Array<OneD, int> vertoffsets(n, 0);
285  Array<OneD, int> edgecounts (n, 0);
286  Array<OneD, int> edgeoffsets(n, 0);
287  vertcounts[p] = vertReorderedGraphVertId.size();
288  edgecounts[p] = edgeReorderedGraphVertId.size();
289  m_comm->AllReduce(vertcounts, LibUtilities::ReduceSum);
290  m_comm->AllReduce(edgecounts, LibUtilities::ReduceSum);
291 
292  for (i = 1; i < n; ++i)
293  {
294  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
295  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
296  }
297 
298  int nTotVerts = Vmath::Vsum(n,vertcounts,1);
299  int nTotEdges = Vmath::Vsum(n,edgecounts,1);
300 
301  Array<OneD, int> vertlist(nTotVerts, 0);
302  Array<OneD, int> edgelist(nTotEdges, 0);
303  std::map<int, int>::iterator it;
304  // construct list of global ids of global vertices
305  for (it = vertReorderedGraphVertId.begin(), i = 0;
306  it != vertReorderedGraphVertId.end();
307  ++it, ++i)
308  {
309  vertlist[vertoffsets[p] + i] = it->first;
310  }
311 
312  // construct list of global ids of global edges
313  for (it = edgeReorderedGraphVertId.begin(), i = 0;
314  it != edgeReorderedGraphVertId.end();
315  ++it, ++i)
316  {
317  edgelist[edgeoffsets[p] + i] = it->first;
318  }
319  m_comm->AllReduce(vertlist, LibUtilities::ReduceSum);
320  m_comm->AllReduce(edgelist, LibUtilities::ReduceSum);
321 
322  // Now we have a list of all Dirichlet vertices and edges on all
323  // processors.
324 
325  int nExtraDirichlet = 0;
326  map<int, int> extraDirVertIds, extraDirEdgeIds;
327 
328  // Ensure Dirchlet vertices are consistently recorded between
329  // processes (e.g. Dirichlet region meets Neumann region across a
330  // partition boundary requires vertex on partition to be Dirichlet).
331  //
332  // To do this we look over all elements and vertices in local
333  // partition and see if they match the values stored in the vertlist
334  // from othe processors and if so record the meshVertId/meshEdgeId
335  // and the processor it comes from.
336  for (i = 0; i < n; ++i)
337  {
338  if (i == p)
339  {
340  continue;
341  }
342 
343  for(j = 0; j < locExpVector.size(); j++)
344  {
345  locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(j)]
346  ->as<LocalRegions::Expansion3D>();
347 
348  for(k = 0; k < locExpansion->GetNverts(); k++)
349  {
350  meshVertId = locExpansion->GetGeom3D()->GetVid(k);
351  if(vertReorderedGraphVertId.count(meshVertId) == 0)
352  {
353  for (l = 0; l < vertcounts[i]; ++l)
354  {
355  if (vertlist[vertoffsets[i]+l] == meshVertId)
356  {
357  extraDirVertIds[meshVertId] = i;
358  vertReorderedGraphVertId[meshVertId] = graphVertId++;
359  nExtraDirichlet++;
360  }
361  }
362  }
363  }
364 
365  for(k = 0; k < locExpansion->GetNedges(); k++)
366  {
367  meshEdgeId = locExpansion->GetGeom3D()->GetEid(k);
368  if(edgeReorderedGraphVertId.count(meshEdgeId) == 0)
369  {
370  for (l = 0; l < edgecounts[i]; ++l)
371  {
372  if (edgelist[edgeoffsets[i]+l] == meshEdgeId)
373  {
374  extraDirEdgeIds[meshEdgeId] = i;
375  edgeReorderedGraphVertId[meshEdgeId] = graphVertId++;
376  nExtraDirichlet += locExpansion->GetEdgeNcoeffs(k) - 2;
377  }
378  }
379  }
380  }
381  }
382  }
383 
384  // Low Energy preconditioner needs to know how many extra Dirichlet
385  // edges are on this process so store map in array.
386  int m_extradiredges = extraDirEdgeIds.size();
387  m_extraDirEdges = Array<OneD, int>(m_extradiredges, -1);
388  i = 0;
389  for (mapConstIt = extraDirEdgeIds.begin();
390  mapConstIt != extraDirEdgeIds.end(); mapConstIt++)
391  {
392  meshEdgeId = mapConstIt->first;
393  m_extraDirEdges[i++] = meshEdgeId;
394  }
395 
396  // Now we have a list of all vertices and edges that are Dirichlet
397  // and not defined on the local partition as well as which processor
398  // they are stored on.
399  //
400  // Make a full list of all such entities on all processors and which
401  // processor they belong to.
402  for (i = 0; i < n; ++i)
403  {
404  vertcounts [i] = 0;
405  vertoffsets[i] = 0;
406  edgecounts [i] = 0;
407  edgeoffsets[i] = 0;
408  }
409 
410  vertcounts[p] = extraDirVertIds.size();
411  edgecounts[p] = extraDirEdgeIds.size();
412  m_comm->AllReduce(vertcounts, LibUtilities::ReduceSum);
413  m_comm->AllReduce(edgecounts, LibUtilities::ReduceSum);
414  nTotVerts = Vmath::Vsum(n, vertcounts, 1);
415  nTotEdges = Vmath::Vsum(n, edgecounts, 1);
416 
417  vertoffsets[0] = edgeoffsets[0] = 0;
418 
419  for (i = 1; i < n; ++i)
420  {
421  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
422  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
423  }
424 
425  Array<OneD, int> vertids (nTotVerts, 0);
426  Array<OneD, int> edgeids (nTotEdges, 0);
427  Array<OneD, int> vertprocs(nTotVerts, 0);
428  Array<OneD, int> edgeprocs(nTotEdges, 0);
429 
430  for (it = extraDirVertIds.begin(), i = 0;
431  it != extraDirVertIds.end(); ++it, ++i)
432  {
433  vertids [vertoffsets[p]+i] = it->first;
434  vertprocs[vertoffsets[p]+i] = it->second;
435  }
436 
437  for (it = extraDirEdgeIds.begin(), i = 0;
438  it != extraDirEdgeIds.end(); ++it, ++i)
439  {
440  edgeids [edgeoffsets[p]+i] = it->first;
441  edgeprocs[edgeoffsets[p]+i] = it->second;
442  }
443 
444  m_comm->AllReduce(vertids, LibUtilities::ReduceSum);
445  m_comm->AllReduce(vertprocs, LibUtilities::ReduceSum);
446  m_comm->AllReduce(edgeids, LibUtilities::ReduceSum);
447  m_comm->AllReduce(edgeprocs, LibUtilities::ReduceSum);
448 
449  set<int> extraDirVerts;
450  set<int> extraDirEdges;
451 
452  // Set up list of vertices that need to be shared to other partitions
453  for (i = 0; i < nTotVerts; ++i)
454  {
455  if (m_comm->GetRank() == vertprocs[i])
456  {
457  extraDirVerts.insert(vertids[i]);
458  }
459  }
460 
461  // Set up list of edges that need to be shared to other partitions
462  for (i = 0; i < nTotEdges; ++i)
463  {
464  if (m_comm->GetRank() == edgeprocs[i])
465  {
466  extraDirEdges.insert(edgeids[i]);
467  }
468  }
469 
470  // Check between processes if the whole system is singular
471  int s = (systemSingular ? 1 : 0);
472  m_comm->AllReduce(s, LibUtilities::ReduceMin);
473  systemSingular = (s == 1 ? true : false);
474 
475  // Count the number of boundary regions on each process
476  Array<OneD, int> bccounts(n, 0);
477  bccounts[p] = bndCondExp.num_elements();
478  m_comm->AllReduce(bccounts, LibUtilities::ReduceSum);
479 
480  // Find the process rank with the maximum number of boundary regions
481  int maxBCIdx = Vmath::Imax(n, bccounts, 1);
482 
483  // If the system is singular, the process with the maximum number of
484  // BCs will set a Dirichlet vertex to make system non-singular.
485  // Note: we find the process with maximum boundary regions to ensure
486  // we do not try to set a Dirichlet vertex on a partition with no
487  // intersection with the boundary.
488  meshVertId = 0;
489  if(systemSingular == true && checkIfSystemSingular && maxBCIdx == p)
490  {
491  if(m_session->DefinesParameter("SingularElement"))
492  {
493  int s_eid;
494  m_session->LoadParameter("SingularElement", s_eid);
495 
496  ASSERTL1(s_eid < locExpVector.size(),
497  "SingularElement Parameter is too large");
498 
499  meshVertId = locExpVector[s_eid]
500  ->as<LocalRegions::Expansion2D>()
501  ->GetGeom2D()->GetVid(0);
502  }
503  else if (bndCondExp.num_elements() == 0)
504  {
505  // All boundaries are periodic.
506  meshVertId = locExpVector[0]->GetGeom()->GetVid(0);
507  }
508  else
509  {
510  // Last region i and j = 0 edge
511  bndCondFaceExp = bndCondExp[bndCondExp.num_elements()-1]
512  ->GetExp(0)->as<LocalRegions::Expansion2D>();
513 
514  // First vertex 0 of the edge
515  meshVertId = bndCondFaceExp->GetGeom2D()->GetVid(0);
516  }
517 
518  if(vertReorderedGraphVertId.count(meshVertId) == 0)
519  {
520  vertReorderedGraphVertId[meshVertId] = graphVertId++;
521  }
522  }
523 
524  m_comm->AllReduce(meshVertId, LibUtilities::ReduceSum);
525 
526  // When running in parallel, we need to ensure that the singular
527  // mesh vertex is communicated to any periodic vertices, otherwise
528  // the system may diverge.
529  if(systemSingular == true && checkIfSystemSingular)
530  {
531  // Firstly, we check that no other processors have this
532  // vertex. If they do, then we mark the vertex as also being
533  // Dirichlet.
534  if (maxBCIdx != p)
535  {
536  for (i = 0; i < locExpVector.size(); ++i)
537  {
538  for (j = 0; j < locExpVector[i]->GetNverts(); ++j)
539  {
540  if (locExpVector[i]->GetGeom()->GetVid(j) !=
541  meshVertId)
542  {
543  continue;
544  }
545 
546  if (vertReorderedGraphVertId.count(meshVertId) == 0)
547  {
548  vertReorderedGraphVertId[meshVertId] =
549  graphVertId++;
550  }
551  }
552  }
553  }
554 
555  // In the case that meshVertId is periodic with other vertices,
556  // this process and all other processes need to make sure that
557  // the periodic vertices are also marked as Dirichlet.
558  int gId;
559 
560  // At least one process (maxBCidx) will have already associated
561  // a graphVertId with meshVertId. Others won't even have any of
562  // the vertices. The logic below is designed to handle both
563  // cases.
564  if (vertReorderedGraphVertId.count(meshVertId) == 0)
565  {
566  gId = -1;
567  }
568  else
569  {
570  gId = vertReorderedGraphVertId[meshVertId];
571  }
572 
573  for (pIt = periodicVerts.begin();
574  pIt != periodicVerts.end(); ++pIt)
575  {
576  // Either the vertex is local to this processor (in which
577  // case it will be in the pIt->first position) or else
578  // meshVertId might be contained within another processor's
579  // vertex list. The if statement below covers both cases. If
580  // we find it, set as Dirichlet with the vertex id gId.
581  if (pIt->first == meshVertId)
582  {
583  vertReorderedGraphVertId[meshVertId] =
584  gId < 0 ? graphVertId++ : gId;
585 
586  for (i = 0; i < pIt->second.size(); ++i)
587  {
588  if (pIt->second[i].isLocal)
589  {
590  vertReorderedGraphVertId[pIt->second[i].id] =
591  gId;
592  }
593  }
594  }
595  else
596  {
597  bool found = false;
598  for (i = 0; i < pIt->second.size(); ++i)
599  {
600  if (pIt->second[i].id == meshVertId)
601  {
602  found = true;
603  break;
604  }
605  }
606 
607  if (found)
608  {
609  vertReorderedGraphVertId[pIt->first] =
610  gId < 0 ? graphVertId++ : gId;
611 
612  for (i = 0; i < pIt->second.size(); ++i)
613  {
614  if (pIt->second[i].isLocal)
615  {
616  vertReorderedGraphVertId[
617  pIt->second[i].id] = gId;
618  }
619  }
620  }
621  }
622  }
623  }
624 
625  m_numLocalDirBndCoeffs = nLocDirBndCondDofs + nExtraDirichlet;
626  firstNonDirGraphVertId = graphVertId;
627 
628  /**
629  * STEP 2: Now order all other vertices and edges in the graph and
630  * create a temporary numbering of domain-interior vertices and
631  * edges.
632  */
633  typedef boost::adjacency_list<
634  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
635  BoostGraph boostGraphObj;
636 
637  map<int, int> vertTempGraphVertId;
638  map<int, int> edgeTempGraphVertId;
639  map<int, int> faceTempGraphVertId;
640  map<int, int> vwgts_map;
641  Array<OneD, int> localVerts;
642  Array<OneD, int> localEdges;
643  Array<OneD, int> localFaces;
644 
645  int tempGraphVertId = 0;
646  int localVertOffset = 0;
647  int localEdgeOffset = 0;
648  int localFaceOffset = 0;
649  int nTotalVerts = 0;
650  int nTotalEdges = 0;
651  int nTotalFaces = 0;
652  int nVerts;
653  int nEdges;
654  int nFaces;
655  int vertCnt;
656  int edgeCnt;
657  int faceCnt;
658 
659  m_numNonDirVertexModes = 0;
660  m_numNonDirEdges = 0;
661  m_numNonDirFaces = 0;
662  m_numNonDirFaceModes = 0;
663  m_numNonDirFaceModes = 0;
664 
665  m_numLocalBndCoeffs = 0;
666 
667  map<int,int> EdgeSize;
668  map<int,int> FaceSize;
669 
670  /// - Count verts, edges, face and add up edges and face sizes
671  for(i = 0; i < locExpVector.size(); ++i)
672  {
673  if((locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
674  ->as<LocalRegions::Expansion3D>()))
675  {
676  nTotalVerts += locExpansion->GetNverts();
677  nTotalEdges += locExpansion->GetNedges();
678  nTotalFaces += locExpansion->GetNfaces();
679 
680  nEdges = locExpansion->GetNedges();
681  for(j = 0; j < nEdges; ++j)
682  {
683  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
684  meshEdgeId = locExpansion->GetGeom3D()->GetEid(j);
685  EdgeSize[meshEdgeId] = nEdgeInteriorCoeffs;
686  }
687 
688  nFaces = locExpansion->GetNfaces();
689  faceCnt = 0;
690  for(j = 0; j < nFaces; ++j)
691  {
692  nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
693  meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
694  FaceSize[meshFaceId] = nFaceInteriorCoeffs;
695  }
696  }
697  }
698 
699 
700  /// - Periodic vertices
701  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
702  {
703  meshVertId = pIt->first;
704 
705  // This periodic vertex is joined to a Dirichlet condition.
706  if (vertReorderedGraphVertId.count(pIt->first) != 0)
707  {
708  for (i = 0; i < pIt->second.size(); ++i)
709  {
710  meshVertId2 = pIt->second[i].id;
711  if (vertReorderedGraphVertId.count(meshVertId2) == 0 &&
712  pIt->second[i].isLocal)
713  {
714  vertReorderedGraphVertId[meshVertId2] =
715  vertReorderedGraphVertId[meshVertId];
716  }
717  }
718  continue;
719  }
720 
721  // One of the attached vertices is Dirichlet.
722  bool isDirichlet = false;
723  for (i = 0; i < pIt->second.size(); ++i)
724  {
725  if (!pIt->second[i].isLocal)
726  {
727  continue;
728  }
729 
730  meshVertId2 = pIt->second[i].id;
731  if (vertReorderedGraphVertId.count(meshVertId2) > 0)
732  {
733  isDirichlet = true;
734  break;
735  }
736  }
737 
738  if (isDirichlet)
739  {
740  vertReorderedGraphVertId[meshVertId] =
741  vertReorderedGraphVertId[pIt->second[i].id];
742 
743  for (j = 0; j < pIt->second.size(); ++j)
744  {
745  meshVertId2 = pIt->second[i].id;
746  if (j == i || !pIt->second[j].isLocal ||
747  vertReorderedGraphVertId.count(meshVertId2) > 0)
748  {
749  continue;
750  }
751 
752  vertReorderedGraphVertId[meshVertId2] =
753  vertReorderedGraphVertId[pIt->second[i].id];
754  }
755 
756  continue;
757  }
758 
759  // Otherwise, see if a vertex ID has already been set.
760  for (i = 0; i < pIt->second.size(); ++i)
761  {
762  if (!pIt->second[i].isLocal)
763  {
764  continue;
765  }
766 
767  if (vertTempGraphVertId.count(pIt->second[i].id) > 0)
768  {
769  break;
770  }
771  }
772 
773  if (i == pIt->second.size())
774  {
775  vertTempGraphVertId[meshVertId] = tempGraphVertId++;
776  m_numNonDirVertexModes += 1;
777  }
778  else
779  {
780  vertTempGraphVertId[meshVertId] =
781  vertTempGraphVertId[pIt->second[i].id];
782  }
783  }
784 
785  // Store the temporary graph vertex id's of all element edges and
786  // vertices in these 3 arrays below
787  localVerts = Array<OneD, int>(nTotalVerts,-1);
788  localEdges = Array<OneD, int>(nTotalEdges,-1);
789  localFaces = Array<OneD, int>(nTotalFaces,-1);
790 
791  // Set up vertex numbering
792  for(i = 0; i < locExpVector.size(); ++i)
793  {
794  if((locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
795  ->as<LocalRegions::Expansion3D>()))
796  {
797  vertCnt = 0;
798  nVerts = locExpansion->GetNverts();
799  for(j = 0; j < nVerts; ++j)
800  {
801  meshVertId = locExpansion->GetGeom3D()->GetVid(j);
802  if(vertReorderedGraphVertId.count(meshVertId) == 0)
803  {
804  if(vertTempGraphVertId.count(meshVertId) == 0)
805  {
806  boost::add_vertex(boostGraphObj);
807  vertTempGraphVertId[meshVertId] = tempGraphVertId++;
808  m_numNonDirVertexModes+=1;
809  }
810  localVerts[localVertOffset+vertCnt++] = vertTempGraphVertId[meshVertId];
811  vwgts_map[ vertTempGraphVertId[meshVertId] ] = 1;
812  }
813  }
814  }
815  else
816  {
817  ASSERTL0(false,"dynamic cast to a local 3D expansion failed");
818  }
819  localVertOffset+=nVerts;
820  }
821 
822  /// - Periodic edges
823  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
824  {
825  meshEdgeId = pIt->first;
826 
827  // This periodic edge is joined to a Dirichlet condition.
828  if (edgeReorderedGraphVertId.count(pIt->first) != 0)
829  {
830  for (i = 0; i < pIt->second.size(); ++i)
831  {
832  meshEdgeId2 = pIt->second[i].id;
833  if (edgeReorderedGraphVertId.count(meshEdgeId2) == 0 &&
834  pIt->second[i].isLocal)
835  {
836  edgeReorderedGraphVertId[meshEdgeId2] =
837  edgeReorderedGraphVertId[meshEdgeId];
838  }
839  }
840  continue;
841  }
842 
843  // One of the attached edges is Dirichlet.
844  bool isDirichlet = false;
845  for (i = 0; i < pIt->second.size(); ++i)
846  {
847  if (!pIt->second[i].isLocal)
848  {
849  continue;
850  }
851 
852  meshEdgeId2 = pIt->second[i].id;
853  if (edgeReorderedGraphVertId.count(meshEdgeId2) > 0)
854  {
855  isDirichlet = true;
856  break;
857  }
858  }
859 
860  if (isDirichlet)
861  {
862  edgeReorderedGraphVertId[meshEdgeId] =
863  edgeReorderedGraphVertId[pIt->second[i].id];
864 
865  for (j = 0; j < pIt->second.size(); ++j)
866  {
867  meshEdgeId2 = pIt->second[i].id;
868  if (j == i || !pIt->second[j].isLocal ||
869  edgeReorderedGraphVertId.count(meshEdgeId2) > 0)
870  {
871  continue;
872  }
873 
874  edgeReorderedGraphVertId[meshEdgeId2] =
875  edgeReorderedGraphVertId[pIt->second[i].id];
876  }
877 
878  continue;
879  }
880 
881  // Otherwise, see if a edge ID has already been set.
882  for (i = 0; i < pIt->second.size(); ++i)
883  {
884  if (!pIt->second[i].isLocal)
885  {
886  continue;
887  }
888 
889  if (edgeTempGraphVertId.count(pIt->second[i].id) > 0)
890  {
891  break;
892  }
893  }
894 
895  if (i == pIt->second.size())
896  {
897  edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
898  nEdgeInteriorCoeffs = EdgeSize[meshEdgeId];
899  m_numNonDirEdgeModes+=nEdgeInteriorCoeffs;
900  m_numNonDirEdges++;
901  }
902  else
903  {
904  edgeTempGraphVertId[meshEdgeId] = edgeTempGraphVertId[pIt->second[i].id];
905  }
906  }
907 
908  // Set up edge numbering
909  for(i = 0; i < locExpVector.size(); ++i)
910  {
911  if((locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
912  ->as<LocalRegions::Expansion3D>()))
913  {
914  edgeCnt = 0;
915  nEdges = locExpansion->GetNedges();
916 
917  for(j = 0; j < nEdges; ++j)
918  {
919  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
920  meshEdgeId = locExpansion->GetGeom3D()->GetEid(j);
921  if(edgeReorderedGraphVertId.count(meshEdgeId) == 0)
922  {
923  if(edgeTempGraphVertId.count(meshEdgeId) == 0)
924  {
925  boost::add_vertex(boostGraphObj);
926  edgeTempGraphVertId[meshEdgeId] = tempGraphVertId++;
927  m_numNonDirEdgeModes+=nEdgeInteriorCoeffs;
928 
929  m_numNonDirEdges++;
930  }
931  localEdges[localEdgeOffset+edgeCnt++] = edgeTempGraphVertId[meshEdgeId];
932  vwgts_map[ edgeTempGraphVertId[meshEdgeId] ] = nEdgeInteriorCoeffs;
933  }
934  }
935  }
936  else
937  {
938  ASSERTL0(false,"dynamic cast to a local 3D expansion failed");
939  }
940  localEdgeOffset+=nEdges;
941  }
942 
943  /// - Periodic faces
944  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
945  {
946  if (!pIt->second[0].isLocal)
947  {
948  // The face mapped to is on another process.
949  meshFaceId = pIt->first;
950  ASSERTL0(faceReorderedGraphVertId.count(meshFaceId) == 0,
951  "This periodic boundary edge has been specified before");
952  faceTempGraphVertId[meshFaceId] = tempGraphVertId++;
953  nFaceInteriorCoeffs = FaceSize[meshFaceId];
954  m_numNonDirFaceModes+=nFaceInteriorCoeffs;
955  m_numNonDirFaces++;
956  }
957  else if (pIt->first < pIt->second[0].id)
958  {
959  ASSERTL0(faceReorderedGraphVertId.count(pIt->first) == 0,
960  "This periodic boundary face has been specified before");
961  ASSERTL0(faceReorderedGraphVertId.count(pIt->second[0].id) == 0,
962  "This periodic boundary face has been specified before");
963 
964  faceTempGraphVertId[pIt->first] = tempGraphVertId;
965  faceTempGraphVertId[pIt->second[0].id] = tempGraphVertId++;
966  nFaceInteriorCoeffs = FaceSize[pIt->first];
967  m_numNonDirFaceModes+=nFaceInteriorCoeffs;
968  m_numNonDirFaces++;
969  }
970  }
971 
972  // setup face numbering
973  for(i = 0; i < locExpVector.size(); ++i)
974  {
975  if((locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
976  ->as<LocalRegions::Expansion3D>()))
977  {
978  nFaces = locExpansion->GetNfaces();
979  faceCnt = 0;
980  for(j = 0; j < nFaces; ++j)
981  {
982  nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
983  meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
984  if(faceReorderedGraphVertId.count(meshFaceId) == 0)
985  {
986  if(faceTempGraphVertId.count(meshFaceId) == 0)
987  {
988  boost::add_vertex(boostGraphObj);
989  faceTempGraphVertId[meshFaceId] = tempGraphVertId++;
990  m_numNonDirFaceModes+=nFaceInteriorCoeffs;
991 
992  m_numNonDirFaces++;
993  }
994  localFaces[localFaceOffset+faceCnt++] = faceTempGraphVertId[meshFaceId];
995  vwgts_map[ faceTempGraphVertId[meshFaceId] ] = nFaceInteriorCoeffs;
996  }
997  }
998  m_numLocalBndCoeffs += locExpansion->NumBndryCoeffs();
999  }
1000  else
1001  {
1002  ASSERTL0(false,"dynamic cast to a local 3D expansion failed");
1003  }
1004  localFaceOffset+=nFaces;
1005  }
1006 
1007  localVertOffset=0;
1008  localEdgeOffset=0;
1009  localFaceOffset=0;
1010  for(i = 0; i < locExpVector.size(); ++i)
1011  {
1012  if((locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
1013  ->as<LocalRegions::Expansion3D>()))
1014  {
1015  nVerts = locExpansion->GetNverts();
1016  nEdges = locExpansion->GetNedges();
1017  nFaces = locExpansion->GetNfaces();
1018 
1019  // Now loop over all local faces, edges and vertices of this
1020  // element and define that all other faces, edges and
1021  // verices of this element are adjacent to them.
1022 
1023  // Vertices
1024  for(j = 0; j < nVerts; j++)
1025  {
1026  if(localVerts[j+localVertOffset]==-1)
1027  {
1028  break;
1029  }
1030  // associate to other vertices
1031  for(k = 0; k < nVerts; k++)
1032  {
1033  if(localVerts[k+localVertOffset]==-1)
1034  {
1035  break;
1036  }
1037  if(k!=j)
1038  {
1039  boost::add_edge( (size_t) localVerts[j+localVertOffset],
1040  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1041  }
1042  }
1043  // associate to other edges
1044  for(k = 0; k < nEdges; k++)
1045  {
1046  if(localEdges[k+localEdgeOffset]==-1)
1047  {
1048  break;
1049  }
1050  boost::add_edge( (size_t) localVerts[j+localVertOffset],
1051  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1052  }
1053  // associate to other faces
1054  for(k = 0; k < nFaces; k++)
1055  {
1056  if(localFaces[k+localFaceOffset]==-1)
1057  {
1058  break;
1059  }
1060  boost::add_edge( (size_t) localVerts[j+localVertOffset],
1061  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1062  }
1063  }
1064 
1065  // Edges
1066  for(j = 0; j < nEdges; j++)
1067  {
1068  if(localEdges[j+localEdgeOffset]==-1)
1069  {
1070  break;
1071  }
1072  // Associate to other edges
1073  for(k = 0; k < nEdges; k++)
1074  {
1075  if(localEdges[k+localEdgeOffset]==-1)
1076  {
1077  break;
1078  }
1079  if(k!=j)
1080  {
1081  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1082  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1083  }
1084  }
1085  // Associate to vertices
1086  for(k = 0; k < nVerts; k++)
1087  {
1088  if(localVerts[k+localVertOffset]==-1)
1089  {
1090  break;
1091  }
1092  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1093  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1094  }
1095  // Associate to faces
1096  for(k = 0; k < nFaces; k++)
1097  {
1098  if(localFaces[k+localFaceOffset]==-1)
1099  {
1100  break;
1101  }
1102  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1103  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1104  }
1105  }
1106 
1107  // Faces
1108  for(j = 0; j < nFaces; j++)
1109  {
1110  if(localFaces[j+localFaceOffset]==-1)
1111  {
1112  break;
1113  }
1114  // Associate to other faces
1115  for(k = 0; k < nFaces; k++)
1116  {
1117  if(localFaces[k+localFaceOffset]==-1)
1118  {
1119  break;
1120  }
1121  if(k!=j)
1122  {
1123  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1124  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1125  }
1126  }
1127  // Associate to vertices
1128  for(k = 0; k < nVerts; k++)
1129  {
1130  if(localVerts[k+localVertOffset]==-1)
1131  {
1132  break;
1133  }
1134  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1135  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1136  }
1137  // Associate to edges
1138  for(k = 0; k < nEdges; k++)
1139  {
1140  if(localEdges[k+localEdgeOffset]==-1)
1141  {
1142  break;
1143  }
1144  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1145  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1146  }
1147  }
1148  }
1149  else
1150  {
1151  ASSERTL0(false,"dynamic cast to a local 3D expansion failed");
1152  }
1153  localVertOffset+=nVerts;
1154  localEdgeOffset+=nEdges;
1155  localFaceOffset+=nFaces;
1156  }
1157 
1158  // Container to store vertices of the graph which correspond to
1159  // degrees of freedom along the boundary.
1160  set<int> partVerts;
1161 
1162  if (m_solnType == eIterativeMultiLevelStaticCond)
1163  {
1164  vector<long> procVerts, procEdges, procFaces;
1165  set <int> foundVerts, foundEdges, foundFaces;
1166 
1167  // Loop over element and construct the procVerts and procEdges
1168  // vectors, which store the geometry IDs of mesh vertices and
1169  // edges respectively which are local to this process.
1170  for(i = cnt = 0; i < locExpVector.size(); ++i)
1171  {
1172  int elmtid = locExp.GetOffset_Elmt_Id(i);
1173  if((locExpansion = locExpVector[elmtid]
1174  ->as<LocalRegions::Expansion3D>()))
1175  {
1176  for (j = 0; j < locExpansion->GetNverts(); ++j)
1177  {
1178  int vid = locExpansion->GetGeom3D()->GetVid(j)+1;
1179 
1180  if (foundVerts.count(vid) == 0)
1181  {
1182  procVerts.push_back(vid);
1183  foundVerts.insert(vid);
1184  }
1185  }
1186 
1187  for (j = 0; j < locExpansion->GetNedges(); ++j)
1188  {
1189  int eid = locExpansion->GetGeom3D()->GetEid(j)+1;
1190 
1191  if (foundEdges.count(eid) == 0)
1192  {
1193  procEdges.push_back(eid);
1194  foundEdges.insert(eid);
1195  }
1196  }
1197 
1198  for (j = 0; j < locExpansion->GetNfaces(); ++j)
1199  {
1200  int fid = locExpansion->GetGeom3D()->GetFid(j)+1;
1201 
1202  if (foundFaces.count(fid) == 0)
1203  {
1204  procFaces.push_back(fid);
1205  foundFaces.insert(fid);
1206  }
1207  }
1208  }
1209  else
1210  {
1211  ASSERTL0(false,
1212  "dynamic cast to a local 3D expansion failed");
1213  }
1214  }
1215 
1216  int unique_verts = foundVerts.size();
1217  int unique_edges = foundEdges.size();
1218  int unique_faces = foundFaces.size();
1219 
1220  // Now construct temporary GS objects. These will be used to
1221  // populate the arrays tmp3 and tmp4 with the multiplicity of
1222  // the vertices and edges respectively to identify those
1223  // vertices and edges which are located on partition boundary.
1224  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1225  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1226  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1227  Gs::gs_data *tmp1 = Gs::Init(vertArray, m_comm);
1228  Gs::gs_data *tmp2 = Gs::Init(edgeArray, m_comm);
1229  Gs::gs_data *tmp3 = Gs::Init(faceArray, m_comm);
1230  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1231  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1232  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1233  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1234  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1235  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1236 
1237  // Finally, fill the partVerts set with all non-Dirichlet
1238  // vertices which lie on a partition boundary.
1239  for (i = 0; i < unique_verts; ++i)
1240  {
1241  if (tmp4[i] > 1.0)
1242  {
1243  if (vertReorderedGraphVertId.count(procVerts[i]-1) == 0)
1244  {
1245  partVerts.insert(vertTempGraphVertId[procVerts[i]-1]);
1246  }
1247  }
1248  }
1249 
1250  for (i = 0; i < unique_edges; ++i)
1251  {
1252  if (tmp5[i] > 1.0)
1253  {
1254  if (edgeReorderedGraphVertId.count(procEdges[i]-1) == 0)
1255  {
1256  partVerts.insert(edgeTempGraphVertId[procEdges[i]-1]);
1257  }
1258  }
1259  }
1260 
1261  for (i = 0; i < unique_faces; ++i)
1262  {
1263  if (tmp6[i] > 1.0)
1264  {
1265  if (faceReorderedGraphVertId.count(procFaces[i]-1) == 0)
1266  {
1267  partVerts.insert(faceTempGraphVertId[procFaces[i]-1]);
1268  }
1269  }
1270  }
1271  }
1272 
1273  /**
1274  * STEP 3: Reorder graph for optimisation.
1275  */
1276  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
1277  int nGraphVerts = tempGraphVertId;
1278  Array<OneD, int> perm(nGraphVerts);
1279  Array<OneD, int> iperm(nGraphVerts);
1280  Array<OneD, int> vwgts(nGraphVerts);
1281  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1282  for(i = 0; i < nGraphVerts; ++i)
1283  {
1284  vwgts[i] = vwgts_map[i];
1285  }
1286 
1287  if(nGraphVerts)
1288  {
1289  switch(m_solnType)
1290  {
1291  case eDirectFullMatrix:
1292  case eIterativeFull:
1293  case eIterativeStaticCond:
1294  case ePETScStaticCond:
1295  case ePETScFullMatrix:
1296  case eXxtStaticCond:
1297  {
1298  NoReordering(boostGraphObj,perm,iperm);
1299  }
1300  break;
1301  case eDirectStaticCond:
1302  {
1303  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1304  }
1305  break;
1306  case ePETScMultiLevelStaticCond:
1307  case eDirectMultiLevelStaticCond:
1308  case eIterativeMultiLevelStaticCond:
1309  case eXxtMultiLevelStaticCond:
1310  {
1311  MultiLevelBisectionReordering(boostGraphObj,perm,iperm,bottomUpGraph,partVerts);
1312  }
1313  break;
1314  default:
1315  {
1316  ASSERTL0(false,"Unrecognised solution type: " + std::string(MultiRegions::GlobalSysSolnTypeMap[m_solnType]));
1317  }
1318  }
1319  }
1320 
1321  // For parallel multi-level static condensation determine the lowest
1322  // static condensation level amongst processors.
1323  if (m_solnType == eIterativeMultiLevelStaticCond)
1324  {
1325  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1326  m_comm->AllReduce(m_lowestStaticCondLevel,
1327  LibUtilities::ReduceMax);
1328  }
1329  else
1330  {
1331  m_lowestStaticCondLevel = 0;
1332  }
1333 
1334  /**
1335  * STEP 4: Fill the #vertReorderedGraphVertId and
1336  * #edgeReorderedGraphVertId with the optimal ordering from boost.
1337  */
1338  for(mapIt = vertTempGraphVertId.begin(); mapIt != vertTempGraphVertId.end(); mapIt++)
1339  {
1340  vertReorderedGraphVertId[mapIt->first] = iperm[mapIt->second] + graphVertId;
1341  }
1342  for(mapIt = edgeTempGraphVertId.begin(); mapIt != edgeTempGraphVertId.end(); mapIt++)
1343  {
1344  edgeReorderedGraphVertId[mapIt->first] = iperm[mapIt->second] + graphVertId;
1345  }
1346  for(mapIt = faceTempGraphVertId.begin(); mapIt != faceTempGraphVertId.end(); mapIt++)
1347  {
1348  faceReorderedGraphVertId[mapIt->first] = iperm[mapIt->second] + graphVertId;
1349  }
1350 
1351 
1352  /**
1353  * STEP 5: Set up an array which contains the offset information of
1354  * the different graph vertices.
1355  *
1356  * This basically means to identify to how many global degrees of
1357  * freedom the individual graph vertices correspond. Obviously,
1358  * the graph vertices corresponding to the mesh-vertices account
1359  * for a single global DOF. However, the graph vertices
1360  * corresponding to the element edges correspond to N-2 global DOF
1361  * where N is equal to the number of boundary modes on this edge.
1362  */
1363  Array<OneD, int> graphVertOffset(vertReorderedGraphVertId.size()+
1364  edgeReorderedGraphVertId.size()+
1365  faceReorderedGraphVertId.size()+1);
1366  graphVertOffset[0] = 0;
1367 
1368  for(i = 0; i < locExpVector.size(); ++i)
1369  {
1370  locExpansion = locExpVector[locExp.GetOffset_Elmt_Id(i)]
1371  ->as<LocalRegions::Expansion3D>();
1372 
1373  for(j = 0; j < locExpansion->GetNverts(); ++j)
1374  {
1375  meshVertId = locExpansion->GetGeom3D()->GetVid(j);
1376  graphVertOffset[vertReorderedGraphVertId[meshVertId]+1] = 1;
1377  }
1378 
1379  for(j = 0; j < locExpansion->GetNedges(); ++j)
1380  {
1381  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j) - 2;
1382  meshEdgeId = locExpansion->GetGeom3D()->GetEid(j);
1383  graphVertOffset[edgeReorderedGraphVertId[meshEdgeId]+1] = nEdgeInteriorCoeffs;
1384 
1385  bType = locExpansion->GetEdgeBasisType(j);
1386  // need a sign vector for modal expansions if nEdgeCoeffs >=4
1387  if( (nEdgeInteriorCoeffs+2 >= 4)&&
1388  ( (bType == LibUtilities::eModified_A)||
1389  (bType == LibUtilities::eModified_B) ) )
1390  {
1391  m_signChange = true;
1392  }
1393  }
1394 
1395  for(j = 0; j < locExpansion->GetNfaces(); ++j)
1396  {
1397  nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
1398  meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
1399  graphVertOffset[faceReorderedGraphVertId[meshFaceId]+1] = nFaceInteriorCoeffs;
1400  }
1401  }
1402 
1403  for(i = 1; i < graphVertOffset.num_elements(); i++)
1404  {
1405  graphVertOffset[i] += graphVertOffset[i-1];
1406  }
1407 
1408  // Allocate the proper amount of space for the class-data
1409  m_numLocalCoeffs = numLocalCoeffs;
1410  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1411  m_localToGlobalMap = Array<OneD, int>(m_numLocalCoeffs,-1);
1412  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
1413  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD,int>(nLocBndCondDofs,-1);
1414  // If required, set up the sign-vector
1415  if(m_signChange)
1416  {
1417  m_localToGlobalSign = Array<OneD, NekDouble>(m_numLocalCoeffs,1.0);
1418  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
1419  m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD,NekDouble>(nLocBndCondDofs,1.0);
1420  }
1421 
1422  m_staticCondLevel = 0;
1423  m_numPatches = locExpVector.size();
1424  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1425  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
1426  for(i = 0; i < m_numPatches; ++i)
1427  {
1428  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1429  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1430  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1431  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1432  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1433  }
1434 
1435  /**
1436  * STEP 6: Now, all ingredients are ready to set up the actual
1437  * local to global mapping.
1438  *
1439  * The remainder of the map consists of the element-interior
1440  * degrees of freedom. This leads to the block-diagonal submatrix
1441  * as each element-interior mode is globally orthogonal to modes
1442  * in all other elements.
1443  */
1444  cnt = 0;
1445 
1446  // Loop over all the elements in the domain
1447  for(i = 0; i < locExpVector.size(); ++i)
1448  {
1449  locExpansion = locExpVector[i]->as<LocalRegions::Expansion3D>();
1450  cnt = locExp.GetCoeff_Offset(i);
1451  for(j = 0; j < locExpansion->GetNverts(); ++j)
1452  {
1453  meshVertId = locExpansion->GetGeom3D()->GetVid(j);
1454 
1455  // Set the global DOF for vertex j of element i
1456  m_localToGlobalMap[cnt+locExpansion->GetVertexMap(j)] =
1457  graphVertOffset[vertReorderedGraphVertId[meshVertId]];
1458  }
1459 
1460  for(j = 0; j < locExpansion->GetNedges(); ++j)
1461  {
1462  nEdgeInteriorCoeffs = locExpansion->GetEdgeNcoeffs(j)-2;
1463  edgeOrient = locExpansion->GetGeom3D()->GetEorient(j);
1464  meshEdgeId = locExpansion->GetGeom3D()->GetEid(j);
1465 
1466  pIt = periodicEdges.find(meshEdgeId);
1467 
1468  // See if this edge is periodic. If it is, then we map all
1469  // edges to the one with lowest ID, and align all
1470  // coefficients to this edge orientation.
1471  if (pIt != periodicEdges.end())
1472  {
1473  int minId = pIt->second[0].id;
1474  int minIdK = 0;
1475  for (k = 1; k < pIt->second.size(); ++k)
1476  {
1477  if (pIt->second[k].id < minId)
1478  {
1479  minId = min(minId, pIt->second[k].id);
1480  minIdK = k;
1481  }
1482  }
1483 
1484  if( meshEdgeId != min(minId, meshEdgeId))
1485  {
1486  if (pIt->second[minIdK].orient == StdRegions::eBackwards)
1487  {
1488  // Swap edge orientation
1489  edgeOrient = (edgeOrient == StdRegions::eForwards) ?
1490  StdRegions::eBackwards : StdRegions::eForwards;
1491  }
1492  }
1493 
1494  }
1495 
1496  locExpansion->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1497 
1498  // Set the global DOF's for the interior modes of edge j
1499  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1500  {
1501  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1502  graphVertOffset[edgeReorderedGraphVertId[meshEdgeId]]+k;
1503  }
1504 
1505  // Fill the sign vector if required
1506  if(m_signChange)
1507  {
1508  for(k = 0; k < nEdgeInteriorCoeffs; ++k)
1509  {
1510  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1511  }
1512  }
1513  }
1514 
1515  for(j = 0; j < locExpansion->GetNfaces(); ++j)
1516  {
1517  nFaceInteriorCoeffs = locExpansion->GetFaceIntNcoeffs(j);
1518  faceOrient = locExpansion->GetGeom3D()->GetFaceOrient(j);
1519  meshFaceId = locExpansion->GetGeom3D()->GetFid(j);
1520 
1521  pIt = periodicFaces.find(meshFaceId);
1522 
1523  if (pIt != periodicFaces.end() &&
1524  meshFaceId == min(meshFaceId, pIt->second[0].id))
1525  {
1526  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1527  }
1528 
1529  locExpansion->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1530 
1531  // Set the global DOF's for the interior modes of face j
1532  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1533  {
1534  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1535  graphVertOffset[faceReorderedGraphVertId[meshFaceId]]+k;
1536  }
1537 
1538  if(m_signChange)
1539  {
1540  for(k = 0; k < nFaceInteriorCoeffs; ++k)
1541  {
1542  m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
1543  }
1544  }
1545  }
1546  }
1547 
1548  // Set up the mapping for the boundary conditions
1549  cnt = 0;
1550  int offset = 0;
1551  for(i = 0; i < bndCondExp.num_elements(); i++)
1552  {
1553  set<int> foundExtraVerts, foundExtraEdges;
1554  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1555  {
1556  bndCondFaceExp = bndCondExp[i]->GetExp(j)
1557  ->as<LocalRegions::Expansion2D>();
1558  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1559  for(k = 0; k < bndCondFaceExp->GetNverts(); k++)
1560  {
1561  meshVertId = bndCondFaceExp->GetGeom2D()->GetVid(k);
1562  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndCondFaceExp->GetVertexMap(k)] = graphVertOffset[vertReorderedGraphVertId[meshVertId]];
1563 
1564  if (bndConditions[i]->GetBoundaryConditionType() !=
1565  SpatialDomains::eDirichlet)
1566  {
1567  continue;
1568  }
1569 
1570  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1571  if (iter != extraDirVerts.end() &&
1572  foundExtraVerts.count(meshVertId) == 0)
1573  {
1574  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1575  bndCondFaceExp->GetVertexMap(k);
1576  int gid = graphVertOffset[
1577  vertReorderedGraphVertId[meshVertId]];
1578  ExtraDirDof t(loc, gid, 1.0);
1579  m_extraDirDofs[i].push_back(t);
1580  foundExtraVerts.insert(meshVertId);
1581  }
1582  }
1583 
1584  for(k = 0; k < bndCondFaceExp->GetNedges(); k++)
1585  {
1586  nEdgeInteriorCoeffs = bndCondFaceExp->GetEdgeNcoeffs(k)-2;
1587  edgeOrient = bndCondFaceExp->GetGeom2D()->GetEorient(k);
1588  meshEdgeId = bndCondFaceExp->GetGeom2D()->GetEid(k);
1589 
1590  pIt = periodicEdges.find(meshEdgeId);
1591 
1592  // See if this edge is periodic. If it is, then we map
1593  // all edges to the one with lowest ID, and align all
1594  // coefficients to this edge orientation.
1595  if (pIt != periodicEdges.end())
1596  {
1597  int minId = pIt->second[0].id;
1598  int minIdL = 0;
1599  for (l = 1; l < pIt->second.size(); ++l)
1600  {
1601  if (pIt->second[l].id < minId)
1602  {
1603  minId = min(minId, pIt->second[l].id);
1604  minIdL = l;
1605  }
1606  }
1607 
1608  if (pIt->second[minIdL].orient == StdRegions::eBackwards &&
1609  meshEdgeId != min(minId, meshEdgeId))
1610  {
1611  edgeOrient = (edgeOrient == StdRegions::eForwards) ? StdRegions::eBackwards : StdRegions::eForwards;
1612  }
1613  }
1614 
1615  bndCondFaceExp->GetEdgeInteriorMap(
1616  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1617 
1618  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1619  {
1620  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1621  graphVertOffset[edgeReorderedGraphVertId[meshEdgeId]]+l;
1622  }
1623 
1624  // Fill the sign vector if required
1625  if(m_signChange)
1626  {
1627  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1628  {
1629  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1630  }
1631  }
1632 
1633  if (bndConditions[i]->GetBoundaryConditionType() !=
1634  SpatialDomains::eDirichlet)
1635  {
1636  continue;
1637  }
1638 
1639  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1640  if (iter != extraDirEdges.end() &&
1641  foundExtraEdges.count(meshEdgeId) == 0 &&
1642  nEdgeInteriorCoeffs > 0)
1643  {
1644  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1645  {
1646  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1647  edgeInteriorMap[l];
1648  int gid = graphVertOffset[
1649  edgeReorderedGraphVertId[meshEdgeId]]+l;
1650  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1651  m_extraDirDofs[i].push_back(t);
1652  }
1653  foundExtraEdges.insert(meshEdgeId);
1654  }
1655  }
1656 
1657  meshFaceId = bndCondFaceExp->GetGeom2D()->GetFid();
1658  intDofCnt = 0;
1659  for(k = 0; k < bndCondFaceExp->GetNcoeffs(); k++)
1660  {
1661  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1662  {
1663  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
1664  graphVertOffset[faceReorderedGraphVertId[meshFaceId]]+intDofCnt;
1665  intDofCnt++;
1666  }
1667  }
1668  }
1669  offset += bndCondExp[i]->GetNcoeffs();
1670  }
1671 
1672  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1673  m_numGlobalBndCoeffs = globalId;
1674 
1675 
1676  /**
1677  * STEP 7: The boundary condition mapping is generated from the
1678  * same vertex renumbering.
1679  */
1680  cnt=0;
1681  for(i = 0; i < m_numLocalCoeffs; ++i)
1682  {
1683  if(m_localToGlobalMap[i] == -1)
1684  {
1685  m_localToGlobalMap[i] = globalId++;
1686  }
1687  else
1688  {
1689  if(m_signChange)
1690  {
1691  m_localToGlobalBndSign[cnt]=m_localToGlobalSign[i];
1692  }
1693  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
1694  }
1695  }
1696 
1697  m_numGlobalCoeffs = globalId;
1698 
1699  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1700 
1701  // Since we can now have multiple entries to m_extraDirDofs due to
1702  // periodic boundary conditions we make a call to work out the
1703  // multiplicity of all entries and invert (Need to be after
1704  // SetupUniversalC0ContMap)
1705  Array<OneD, NekDouble> valence(m_numGlobalBndCoeffs,0.0);
1706 
1707  // Fill in Dirichlet coefficients that are to be sent to other
1708  // processors with a value of 1
1709  map<int, vector<ExtraDirDof> >::iterator Tit;
1710 
1711  // Generate valence for extraDirDofs
1712  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1713  {
1714  for (i = 0; i < Tit->second.size(); ++i)
1715  {
1716  valence[Tit->second[i].get<1>()] = 1.0;
1717  }
1718  }
1719 
1720  UniversalAssembleBnd(valence);
1721 
1722  // Set third argument of tuple to inverse of valence.
1723  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1724  {
1725  for (i = 0; i < Tit->second.size(); ++i)
1726  {
1727  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1728  }
1729  }
1730 
1731 
1732  // Set up the local to global map for the next level when using
1733  // multi-level static condensation
1734  if ((m_solnType == eDirectMultiLevelStaticCond ||
1735  m_solnType == eIterativeMultiLevelStaticCond) && nGraphVerts)
1736  {
1737  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1738  {
1739  Array<OneD, int> vwgts_perm(nGraphVerts);
1740  for(i = 0; i < nGraphVerts; ++i)
1741  {
1742  vwgts_perm[i] = vwgts[perm[i]];
1743  }
1744 
1745  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1746  m_nextLevelLocalToGlobalMap = MemoryManager<
1747  AssemblyMap>::AllocateSharedPtr(this,bottomUpGraph);
1748  }
1749  }
1750 
1751  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
1752  m_localToGlobalMap.end());
1753 
1754  // Add up hash values if parallel
1755  int hash = m_hash;
1756  m_comm->GetRowComm()->AllReduce(hash,
1757  LibUtilities::ReduceSum);
1758  m_hash = hash;
1759  }
1760  } // namespace
1761 } // namespace
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2