Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
AssemblyMapCG.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMapCG.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 Local to Global mapping routines, base class
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ExpList.h>
38 #include <LocalRegions/Expansion.h>
42 
43 
44 #include <boost/config.hpp>
45 #include <boost/graph/adjacency_list.hpp>
46 #include <boost/graph/cuthill_mckee_ordering.hpp>
47 #include <boost/graph/properties.hpp>
48 #include <boost/graph/bandwidth.hpp>
49 
50 using namespace std;
51 
52 namespace Nektar
53 {
54  namespace MultiRegions
55  {
56  /**
57  * @class AssemblyMapCG
58  * Mappings are created for three possible global solution types:
59  * - Direct full matrix
60  * - Direct static condensation
61  * - Direct multi-level static condensation
62  * In the latter case, mappings are created recursively for the
63  * different levels of static condensation.
64  *
65  * These mappings are used by GlobalLinSys to generate the global
66  * system.
67  */
68 
69  /**
70  *
71  */
72  AssemblyMapCG::AssemblyMapCG(
74  const std::string variable):
75  AssemblyMap(pSession,variable)
76  {
77  pSession->LoadParameter(
78  "MaxStaticCondLevel", m_maxStaticCondLevel, 100);
79  }
80 
82  const ExpList &locExp,
83  const BndCondExp &bndCondExp,
84  const Array<OneD, const BndCond> &bndConditions,
85  const bool checkIfSystemSingular,
86  const PeriodicMap &periodicVerts,
87  const PeriodicMap &periodicEdges,
88  const PeriodicMap &periodicFaces,
89  DofGraph &graph,
91  set<int> &extraDirVerts,
92  set<int> &extraDirEdges,
93  int &firstNonDirGraphVertId,
94  int &nExtraDirichlet,
95  int mdswitch)
96  {
97  int graphVertId = 0;
98  int vMaxVertId = -1;
99  int i, j, k, l, cnt;
100  int meshVertId, meshEdgeId, meshFaceId;
101  int meshVertId2, meshEdgeId2;
102 
104  const LocalRegions::ExpansionVector &locExpVector =
105  *(locExp.GetExp());
106  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
107  PeriodicMap::const_iterator pIt;
108 
110  m_systemSingular = checkIfSystemSingular;
111 
112  for(i = 0; i < bndCondExp.num_elements(); i++)
113  {
114  // Check to see if any value on boundary has Dirichlet value.
115  cnt = 0;
116  for(k = 0; k < bndConditions.num_elements(); ++k)
117  {
118  if (bndConditions[k][i]->GetBoundaryConditionType() ==
120  {
121  cnt++;
122  }
123  if (bndConditions[k][i]->GetBoundaryConditionType() !=
125  {
126  m_systemSingular = false;
127  }
128  }
129 
130  // Find the maximum boundary vertex ID on this process. This is
131  // used later to pin a vertex if the system is singular.
132  for (j = 0; j < bndCondExp[i]->GetNumElmts(); ++j)
133  {
134  bndExp = bndCondExp[i]->GetExp(j)->as<LocalRegions::Expansion>();
135  for (k = 0; k < bndExp->GetNverts(); ++k)
136  {
137  if (vMaxVertId < bndExp->GetGeom()->GetVid(k))
138  {
139  vMaxVertId = bndExp->GetGeom()->GetVid(k);
140  }
141  }
142 
143  }
144 
145  // If all boundaries are Dirichlet take out of mask
146  if(cnt == bndConditions.num_elements())
147  {
148  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
149  {
150  bndExp = bndCondExp[i]->GetExp(j);
151 
152  for (k = 0; k < bndExp->GetNverts(); k++)
153  {
154  meshVertId = bndExp->GetGeom()->GetVid(k);
155  if (graph[0].count(meshVertId) == 0)
156  {
157  graph[0][meshVertId] = graphVertId++;
158  }
159  }
160 
161  for (k = 0; k < bndExp->GetNedges(); k++)
162  {
163  meshEdgeId = bndExp->GetGeom()->GetEid(k);
164  if (graph[1].count(meshEdgeId) == 0)
165  {
166  graph[1][meshEdgeId] = graphVertId++;
167  }
168  }
169 
170  // Possibly not a face.
171  meshFaceId = bndExp->GetGeom()->GetGlobalID();
172  const int bndDim = bndExp->GetNumBases();
173  if (graph[bndDim].count(meshFaceId) == 0)
174  {
175  graph[bndDim][meshFaceId] = graphVertId++;
176  }
177  m_numLocalDirBndCoeffs += bndExp->GetNcoeffs();
178  }
179  }
180 
181  m_numLocalBndCondCoeffs += bndCondExp[i]->GetNcoeffs();
182  }
183 
184  // Number of dirichlet edges and faces (not considering periodic
185  // BCs)
186  m_numDirEdges = graph[1].size();
187  m_numDirFaces = graph[2].size();
188 
189  /*
190  * The purpose of this routine is to deal with those degrees of
191  * freedom that are Dirichlet, but do not have a local Dirichlet
192  * boundary condition expansion set.
193  *
194  * For example, in 2D, consider a triangulation of a square into two
195  * triangles. Now imagine one edge of the square is Dirichlet and
196  * the problem is run on two processors. On one processor, one
197  * triangle vertex is Dirichlet, but doesn't know this since the
198  * Dirichlet composite lives on the other processor.
199  *
200  * When the global linear system is solved therefore, there is an
201  * inconsistency that at best leads to an inaccurate answer or a
202  * divergence of the system.
203  *
204  * This routine identifies such cases for 2D, and also for 3D where
205  * e.g. edges may have the same problem (consider an extrusion of
206  * the case above, for example).
207  */
208 
209  // Collate information on Dirichlet vertices from all processes
210  int n = vComm->GetSize();
211  int p = vComm->GetRank();
212 
213  // At this point, graph only contains information from Dirichlet
214  // boundaries. Therefore make a global list of the vert and edge
215  // information on all processors.
216  Array<OneD, int> vertcounts (n, 0);
217  Array<OneD, int> vertoffsets(n, 0);
218  Array<OneD, int> edgecounts (n, 0);
219  Array<OneD, int> edgeoffsets(n, 0);
220  vertcounts[p] = graph[0].size();
221  edgecounts[p] = graph[1].size();
222  vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
223  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
224 
225  for (i = 1; i < n; ++i)
226  {
227  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
228  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
229  }
230 
231  int nTotVerts = Vmath::Vsum(n,vertcounts,1);
232  int nTotEdges = Vmath::Vsum(n,edgecounts,1);
233 
234  Array<OneD, int> vertlist(nTotVerts, 0);
235  Array<OneD, int> edgelist(nTotEdges, 0);
237 
238  // construct list of global ids of global vertices
239  for (it = graph[0].begin(), i = 0;
240  it != graph[0].end();
241  ++it, ++i)
242  {
243  vertlist[vertoffsets[p] + i] = it->first;
244  }
245 
246  // construct list of global ids of global edges
247  for (it = graph[1].begin(), i = 0;
248  it != graph[1].end();
249  ++it, ++i)
250  {
251  edgelist[edgeoffsets[p] + i] = it->first;
252  }
253  vComm->AllReduce(vertlist, LibUtilities::ReduceSum);
254  vComm->AllReduce(edgelist, LibUtilities::ReduceSum);
255 
256  // Now we have a list of all Dirichlet vertices and edges on all
257  // processors.
258  nExtraDirichlet = 0;
259  map<int, int> extraDirVertIds, extraDirEdgeIds;
260 
261  // Ensure Dirchlet vertices are consistently recorded between
262  // processes (e.g. Dirichlet region meets Neumann region across a
263  // partition boundary requires vertex on partition to be Dirichlet).
264  //
265  // To do this we look over all elements and vertices in local
266  // partition and see if they match the values stored in the vertlist
267  // from othe processors and if so record the meshVertId/meshEdgeId
268  // and the processor it comes from.
269  for (i = 0; i < n; ++i)
270  {
271  if (i == p)
272  {
273  continue;
274  }
275 
276  for(j = 0; j < locExpVector.size(); j++)
277  {
278  exp = locExpVector[locExp.GetOffset_Elmt_Id(j)];
279 
280  for(k = 0; k < exp->GetNverts(); k++)
281  {
282  meshVertId = exp->GetGeom()->GetVid(k);
283  if(graph[0].count(meshVertId) == 0)
284  {
285  for (l = 0; l < vertcounts[i]; ++l)
286  {
287  if (vertlist[vertoffsets[i]+l] == meshVertId)
288  {
289  extraDirVertIds[meshVertId] = i;
290  graph[0][meshVertId] = graphVertId++;
291  nExtraDirichlet++;
292  }
293  }
294  }
295  }
296 
297  for(k = 0; k < exp->GetNedges(); k++)
298  {
299  meshEdgeId = exp->GetGeom()->GetEid(k);
300  if(graph[1].count(meshEdgeId) == 0)
301  {
302  for (l = 0; l < edgecounts[i]; ++l)
303  {
304  if (edgelist[edgeoffsets[i]+l] == meshEdgeId)
305  {
306  extraDirEdgeIds[meshEdgeId] = i;
307  graph[1][meshEdgeId] = graphVertId++;
308  nExtraDirichlet += exp->GetEdgeNcoeffs(k)-2;
309  }
310  }
311  }
312  }
313  }
314  }
315 
316  // Low Energy preconditioner needs to know how many extra Dirichlet
317  // edges are on this process so store map in array.
318  map<int, int>::const_iterator mapConstIt;
319  m_extraDirEdges = Array<OneD, int>(extraDirEdgeIds.size(), -1);
320  for (mapConstIt = extraDirEdgeIds.begin(), i = 0;
321  mapConstIt != extraDirEdgeIds.end(); mapConstIt++)
322  {
323  meshEdgeId = mapConstIt->first;
324  m_extraDirEdges[i++] = meshEdgeId;
325  }
326 
327  // Now we have a list of all vertices and edges that are Dirichlet
328  // and not defined on the local partition as well as which processor
329  // they are stored on.
330  //
331  // Make a full list of all such entities on all processors and which
332  // processor they belong to.
333  for (i = 0; i < n; ++i)
334  {
335  vertcounts [i] = 0;
336  vertoffsets[i] = 0;
337  edgecounts [i] = 0;
338  edgeoffsets[i] = 0;
339  }
340 
341  vertcounts[p] = extraDirVertIds.size();
342  edgecounts[p] = extraDirEdgeIds.size();
343  vComm->AllReduce(vertcounts, LibUtilities::ReduceSum);
344  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
345  nTotVerts = Vmath::Vsum(n, vertcounts, 1);
346  nTotEdges = Vmath::Vsum(n, edgecounts, 1);
347 
348  vertoffsets[0] = edgeoffsets[0] = 0;
349 
350  for (i = 1; i < n; ++i)
351  {
352  vertoffsets[i] = vertoffsets[i-1] + vertcounts[i-1];
353  edgeoffsets[i] = edgeoffsets[i-1] + edgecounts[i-1];
354  }
355 
356  Array<OneD, int> vertids (nTotVerts, 0);
357  Array<OneD, int> edgeids (nTotEdges, 0);
358  Array<OneD, int> vertprocs(nTotVerts, 0);
359  Array<OneD, int> edgeprocs(nTotEdges, 0);
360 
361  for (it = extraDirVertIds.begin(), i = 0;
362  it != extraDirVertIds.end(); ++it, ++i)
363  {
364  vertids [vertoffsets[p]+i] = it->first;
365  vertprocs[vertoffsets[p]+i] = it->second;
366  }
367 
368  for (it = extraDirEdgeIds.begin(), i = 0;
369  it != extraDirEdgeIds.end(); ++it, ++i)
370  {
371  edgeids [edgeoffsets[p]+i] = it->first;
372  edgeprocs[edgeoffsets[p]+i] = it->second;
373  }
374 
375  vComm->AllReduce(vertids, LibUtilities::ReduceSum);
376  vComm->AllReduce(vertprocs, LibUtilities::ReduceSum);
377  vComm->AllReduce(edgeids, LibUtilities::ReduceSum);
378  vComm->AllReduce(edgeprocs, LibUtilities::ReduceSum);
379 
380  // Set up list of vertices that need to be shared to other
381  // partitions
382  for (i = 0; i < nTotVerts; ++i)
383  {
384  if (vComm->GetRank() == vertprocs[i])
385  {
386  extraDirVerts.insert(vertids[i]);
387  }
388  }
389 
390  // Set up list of edges that need to be shared to other partitions
391  for (i = 0; i < nTotEdges; ++i)
392  {
393  if (vComm->GetRank() == edgeprocs[i])
394  {
395  extraDirEdges.insert(edgeids[i]);
396  }
397  }
398 
399  // Check between processes if the whole system is singular
400  int s = m_systemSingular ? 1 : 0;
401  vComm->AllReduce(s, LibUtilities::ReduceMin);
402  m_systemSingular = s == 1 ? true : false;
403 
404  // Find the minimum boundary vertex ID on each process
405  Array<OneD, int> bcminvertid(n, 0);
406  bcminvertid[p] = vMaxVertId;
407  vComm->AllReduce(bcminvertid, LibUtilities::ReduceMax);
408 
409  // Find the process rank with the minimum boundary vertex ID
410  int maxIdx = Vmath::Imax(n, bcminvertid, 1);
411 
412  // If the system is singular, the process with the maximum number of
413  // BCs will set a Dirichlet vertex to make system non-singular.
414  // Note: we find the process with maximum boundary regions to ensure
415  // we do not try to set a Dirichlet vertex on a partition with no
416  // intersection with the boundary.
417  meshVertId = 0;
418 
419  if (m_systemSingular && checkIfSystemSingular && maxIdx == p)
420  {
421  if (m_session->DefinesParameter("SingularVertex"))
422  {
423  m_session->LoadParameter("SingularVertex", meshVertId);
424  }
425  else if (vMaxVertId == -1)
426  {
427  // All boundaries are periodic.
428  meshVertId = locExpVector[0]->GetGeom()->GetVid(0);
429  }
430  else
431  {
432  // Set pinned vertex to that with minimum vertex ID to
433  // ensure consistency in parallel.
434  meshVertId = bcminvertid[p];
435  }
436 
437  if (graph[0].count(meshVertId) == 0)
438  {
439  graph[0][meshVertId] = graphVertId++;
440  }
441  }
442 
443  vComm->AllReduce(meshVertId, LibUtilities::ReduceSum);
444 
445  // When running in parallel, we need to ensure that the singular
446  // mesh vertex is communicated to any periodic vertices, otherwise
447  // the system may diverge.
448  if(m_systemSingular && checkIfSystemSingular)
449  {
450  // Firstly, we check that no other processors have this
451  // vertex. If they do, then we mark the vertex as also being
452  // Dirichlet.
453  if (maxIdx != p)
454  {
455  for (i = 0; i < locExpVector.size(); ++i)
456  {
457  for (j = 0; j < locExpVector[i]->GetNverts(); ++j)
458  {
459  if (locExpVector[i]->GetGeom()->GetVid(j) !=
460  meshVertId)
461  {
462  continue;
463  }
464 
465  if (graph[0].count(meshVertId) == 0)
466  {
467  graph[0][meshVertId] =
468  graphVertId++;
469  }
470  }
471  }
472  }
473 
474  // In the case that meshVertId is periodic with other vertices,
475  // this process and all other processes need to make sure that
476  // the periodic vertices are also marked as Dirichlet.
477  int gId;
478 
479  // At least one process (maxBCidx) will have already associated
480  // a graphVertId with meshVertId. Others won't even have any of
481  // the vertices. The logic below is designed to handle both
482  // cases.
483  if (graph[0].count(meshVertId) == 0)
484  {
485  gId = -1;
486  }
487  else
488  {
489  gId = graph[0][meshVertId];
490  }
491 
492  for (pIt = periodicVerts.begin();
493  pIt != periodicVerts.end(); ++pIt)
494  {
495  // Either the vertex is local to this processor (in which
496  // case it will be in the pIt->first position) or else
497  // meshVertId might be contained within another processor's
498  // vertex list. The if statement below covers both cases. If
499  // we find it, set as Dirichlet with the vertex id gId.
500  if (pIt->first == meshVertId)
501  {
502  gId = gId < 0 ? graphVertId++ : gId;
503  graph[0][meshVertId] = gId;
504 
505  for (i = 0; i < pIt->second.size(); ++i)
506  {
507  if (pIt->second[i].isLocal)
508  {
509  graph[0][pIt->second[i].id] = graph[0][meshVertId];
510  }
511  }
512  }
513  else
514  {
515  bool found = false;
516  for (i = 0; i < pIt->second.size(); ++i)
517  {
518  if (pIt->second[i].id == meshVertId)
519  {
520  found = true;
521  break;
522  }
523  }
524 
525  if (found)
526  {
527  gId = gId < 0 ? graphVertId++ : gId;
528  graph[0][pIt->first] = gId;
529 
530  for (i = 0; i < pIt->second.size(); ++i)
531  {
532  if (pIt->second[i].isLocal)
533  {
534  graph[0][pIt->second[i].id] = graph[0][pIt->first];
535  }
536  }
537  }
538  }
539  }
540  }
541 
542  // Add extra dirichlet boundary conditions to count.
543  m_numLocalDirBndCoeffs += nExtraDirichlet;
544  firstNonDirGraphVertId = graphVertId;
545 
546  typedef boost::adjacency_list<
547  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
548  BoostGraph boostGraphObj;
549 
550  vector<map<int,int> > tempGraph(3);
551  map<int, int> vwgts_map;
552  Array<OneD, int> localVerts;
553  Array<OneD, int> localEdges;
554  Array<OneD, int> localFaces;
555 
556  int tempGraphVertId = 0;
557  int localVertOffset = 0;
558  int localEdgeOffset = 0;
559  int localFaceOffset = 0;
560  int nTotalVerts = 0;
561  int nTotalEdges = 0;
562  int nTotalFaces = 0;
563  int nVerts;
564  int nEdges;
565  int nFaces;
566  int vertCnt;
567  int edgeCnt;
568  int faceCnt;
569 
571  m_numNonDirEdges = 0;
572  m_numNonDirFaces = 0;
576 
577  map<int,int> EdgeSize;
578  map<int,int> FaceSize;
579 
580  /// - Count verts, edges, face and add up edges and face sizes
581  for(i = 0; i < locExpVector.size(); ++i)
582  {
583  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
584  nTotalVerts += exp->GetNverts();
585  nTotalEdges += exp->GetNedges();
586  nTotalFaces += exp->GetNfaces();
587 
588  nEdges = exp->GetNedges();
589  for(j = 0; j < nEdges; ++j)
590  {
591  meshEdgeId = exp->GetGeom()->GetEid(j);
592  if (EdgeSize.count(meshEdgeId) > 0)
593  {
594  EdgeSize[meshEdgeId] =
595  min(EdgeSize[meshEdgeId],
596  exp->GetEdgeNcoeffs(j) - 2);
597  }
598  else
599  {
600  EdgeSize[meshEdgeId] = exp->GetEdgeNcoeffs(j) - 2;
601  }
602  }
603 
604  nFaces = exp->GetNfaces();
605  faceCnt = 0;
606  for(j = 0; j < nFaces; ++j)
607  {
608  meshFaceId = exp->GetGeom()->GetFid(j);
609  if (FaceSize.count(meshFaceId) > 0)
610  {
611  FaceSize[meshFaceId] =
612  min(FaceSize[meshFaceId],
613  exp->GetFaceIntNcoeffs(j));
614  }
615  else
616  {
617  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
618  }
619  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
620  }
621  }
622 
623  /// - Periodic vertices
624  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
625  {
626  meshVertId = pIt->first;
627 
628  // This periodic vertex is joined to a Dirichlet condition.
629  if (graph[0].count(pIt->first) != 0)
630  {
631  for (i = 0; i < pIt->second.size(); ++i)
632  {
633  meshVertId2 = pIt->second[i].id;
634  if (graph[0].count(meshVertId2) == 0 &&
635  pIt->second[i].isLocal)
636  {
637  graph[0][meshVertId2] =
638  graph[0][meshVertId];
639  }
640  }
641  continue;
642  }
643 
644  // One of the attached vertices is Dirichlet.
645  bool isDirichlet = false;
646  for (i = 0; i < pIt->second.size(); ++i)
647  {
648  if (!pIt->second[i].isLocal)
649  {
650  continue;
651  }
652 
653  meshVertId2 = pIt->second[i].id;
654  if (graph[0].count(meshVertId2) > 0)
655  {
656  isDirichlet = true;
657  break;
658  }
659  }
660 
661  if (isDirichlet)
662  {
663  graph[0][meshVertId] =
664  graph[0][pIt->second[i].id];
665 
666  for (j = 0; j < pIt->second.size(); ++j)
667  {
668  meshVertId2 = pIt->second[i].id;
669  if (j == i || !pIt->second[j].isLocal ||
670  graph[0].count(meshVertId2) > 0)
671  {
672  continue;
673  }
674 
675  graph[0][meshVertId2] =
676  graph[0][pIt->second[i].id];
677  }
678 
679  continue;
680  }
681 
682  // Otherwise, see if a vertex ID has already been set.
683  for (i = 0; i < pIt->second.size(); ++i)
684  {
685  if (!pIt->second[i].isLocal)
686  {
687  continue;
688  }
689 
690  if (tempGraph[0].count(pIt->second[i].id) > 0)
691  {
692  break;
693  }
694  }
695 
696  if (i == pIt->second.size())
697  {
698  boost::add_vertex(boostGraphObj);
699  tempGraph[0][meshVertId] = tempGraphVertId++;
701  }
702  else
703  {
704  tempGraph[0][meshVertId] = tempGraph[0][pIt->second[i].id];
705  }
706  }
707 
708  // Store the temporary graph vertex id's of all element edges and
709  // vertices in these 3 arrays below
710  localVerts = Array<OneD, int>(nTotalVerts,-1);
711  localEdges = Array<OneD, int>(nTotalEdges,-1);
712  localFaces = Array<OneD, int>(nTotalFaces,-1);
713 
714  // Set up vertex numbering
715  for(i = 0; i < locExpVector.size(); ++i)
716  {
717  exp = locExpVector[i];
718  vertCnt = 0;
719  nVerts = exp->GetNverts();
720  for(j = 0; j < nVerts; ++j)
721  {
722  meshVertId = exp->GetGeom()->GetVid(j);
723  if(graph[0].count(meshVertId) == 0)
724  {
725  if(tempGraph[0].count(meshVertId) == 0)
726  {
727  boost::add_vertex(boostGraphObj);
728  tempGraph[0][meshVertId] = tempGraphVertId++;
730  }
731  localVerts[localVertOffset+vertCnt++] = tempGraph[0][meshVertId];
732  vwgts_map[ tempGraph[0][meshVertId] ] = 1;
733  }
734  }
735 
736  localVertOffset+=nVerts;
737  }
738 
739  /// - Periodic edges
740  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
741  {
742  meshEdgeId = pIt->first;
743 
744  // This periodic edge is joined to a Dirichlet condition.
745  if (graph[1].count(pIt->first) != 0)
746  {
747  for (i = 0; i < pIt->second.size(); ++i)
748  {
749  meshEdgeId2 = pIt->second[i].id;
750  if (graph[1].count(meshEdgeId2) == 0 &&
751  pIt->second[i].isLocal)
752  {
753  graph[1][meshEdgeId2] =
754  graph[1][meshEdgeId];
755  }
756  }
757  continue;
758  }
759 
760  // One of the attached edges is Dirichlet.
761  bool isDirichlet = false;
762  for (i = 0; i < pIt->second.size(); ++i)
763  {
764  if (!pIt->second[i].isLocal)
765  {
766  continue;
767  }
768 
769  meshEdgeId2 = pIt->second[i].id;
770  if (graph[1].count(meshEdgeId2) > 0)
771  {
772  isDirichlet = true;
773  break;
774  }
775  }
776 
777  if (isDirichlet)
778  {
779  graph[1][meshEdgeId] =
780  graph[1][pIt->second[i].id];
781 
782  for (j = 0; j < pIt->second.size(); ++j)
783  {
784  meshEdgeId2 = pIt->second[i].id;
785  if (j == i || !pIt->second[j].isLocal ||
786  graph[1].count(meshEdgeId2) > 0)
787  {
788  continue;
789  }
790 
791  graph[1][meshEdgeId2] =
792  graph[1][pIt->second[i].id];
793  }
794 
795  continue;
796  }
797 
798  // Otherwise, see if a edge ID has already been set.
799  for (i = 0; i < pIt->second.size(); ++i)
800  {
801  if (!pIt->second[i].isLocal)
802  {
803  continue;
804  }
805 
806  if (tempGraph[1].count(pIt->second[i].id) > 0)
807  {
808  break;
809  }
810  }
811 
812  if (i == pIt->second.size())
813  {
814  boost::add_vertex(boostGraphObj);
815  tempGraph[1][meshEdgeId] = tempGraphVertId++;
816  m_numNonDirEdgeModes += EdgeSize[meshEdgeId];
818  }
819  else
820  {
821  tempGraph[1][meshEdgeId] = tempGraph[1][pIt->second[i].id];
822  }
823  }
824 
825  int nEdgeIntCoeffs, nFaceIntCoeffs;
826 
827  // Set up edge numbering
828  for(i = 0; i < locExpVector.size(); ++i)
829  {
830  exp = locExpVector[i];
831  edgeCnt = 0;
832  nEdges = exp->GetNedges();
833 
834  for(j = 0; j < nEdges; ++j)
835  {
836  meshEdgeId = exp->GetGeom()->GetEid(j);
837  nEdgeIntCoeffs = EdgeSize[meshEdgeId];
838  if(graph[1].count(meshEdgeId) == 0)
839  {
840  if(tempGraph[1].count(meshEdgeId) == 0)
841  {
842  boost::add_vertex(boostGraphObj);
843  tempGraph[1][meshEdgeId] = tempGraphVertId++;
844  m_numNonDirEdgeModes+=nEdgeIntCoeffs;
845 
847  }
848  localEdges[localEdgeOffset+edgeCnt++] = tempGraph[1][meshEdgeId];
849  vwgts_map[ tempGraph[1][meshEdgeId] ] = nEdgeIntCoeffs;
850  }
851  }
852 
853  localEdgeOffset+=nEdges;
854  }
855 
856  /// - Periodic faces
857  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
858  {
859  if (!pIt->second[0].isLocal)
860  {
861  // The face mapped to is on another process.
862  meshFaceId = pIt->first;
863  ASSERTL0(graph[2].count(meshFaceId) == 0,
864  "This periodic boundary edge has been specified before");
865  boost::add_vertex(boostGraphObj);
866  tempGraph[2][meshFaceId] = tempGraphVertId++;
867  nFaceIntCoeffs = FaceSize[meshFaceId];
868  m_numNonDirFaceModes+=nFaceIntCoeffs;
870  }
871  else if (pIt->first < pIt->second[0].id)
872  {
873  ASSERTL0(graph[2].count(pIt->first) == 0,
874  "This periodic boundary face has been specified before");
875  ASSERTL0(graph[2].count(pIt->second[0].id) == 0,
876  "This periodic boundary face has been specified before");
877 
878  boost::add_vertex(boostGraphObj);
879  tempGraph[2][pIt->first] = tempGraphVertId;
880  tempGraph[2][pIt->second[0].id] = tempGraphVertId++;
881  nFaceIntCoeffs = FaceSize[pIt->first];
882  m_numNonDirFaceModes+=nFaceIntCoeffs;
884  }
885  }
886 
887  // setup face numbering
888  for(i = 0; i < locExpVector.size(); ++i)
889  {
890  exp = locExpVector[i];
891  nFaces = exp->GetNfaces();
892  faceCnt = 0;
893  for(j = 0; j < nFaces; ++j)
894  {
895  nFaceIntCoeffs = exp->GetFaceIntNcoeffs(j);
896  meshFaceId = exp->GetGeom()->GetFid(j);
897  if(graph[2].count(meshFaceId) == 0)
898  {
899  if(tempGraph[2].count(meshFaceId) == 0)
900  {
901  boost::add_vertex(boostGraphObj);
902  tempGraph[2][meshFaceId] = tempGraphVertId++;
903  m_numNonDirFaceModes+=nFaceIntCoeffs;
904 
906  }
907  localFaces[localFaceOffset+faceCnt++] = tempGraph[2][meshFaceId];
908  vwgts_map[ tempGraph[2][meshFaceId] ] = nFaceIntCoeffs;
909  }
910  }
911  m_numLocalBndCoeffs += exp->NumBndryCoeffs();
912 
913  localFaceOffset+=nFaces;
914  }
915 
916  localVertOffset=0;
917  localEdgeOffset=0;
918  localFaceOffset=0;
919  for(i = 0; i < locExpVector.size(); ++i)
920  {
921  exp = locExpVector[i];
922  nVerts = exp->GetNverts();
923  nEdges = exp->GetNedges();
924  nFaces = exp->GetNfaces();
925 
926  // Now loop over all local faces, edges and vertices of this
927  // element and define that all other faces, edges and verices of
928  // this element are adjacent to them.
929 
930  // Vertices
931  for(j = 0; j < nVerts; j++)
932  {
933  if(localVerts[j+localVertOffset]==-1)
934  {
935  break;
936  }
937  // associate to other vertices
938  for(k = 0; k < nVerts; k++)
939  {
940  if(localVerts[k+localVertOffset]==-1)
941  {
942  break;
943  }
944  if(k!=j)
945  {
946  boost::add_edge( (size_t) localVerts[j+localVertOffset],
947  (size_t) localVerts[k+localVertOffset],boostGraphObj);
948  }
949  }
950  // associate to other edges
951  for(k = 0; k < nEdges; k++)
952  {
953  if(localEdges[k+localEdgeOffset]==-1)
954  {
955  break;
956  }
957  boost::add_edge( (size_t) localVerts[j+localVertOffset],
958  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
959  }
960  // associate to other faces
961  for(k = 0; k < nFaces; k++)
962  {
963  if(localFaces[k+localFaceOffset]==-1)
964  {
965  break;
966  }
967  boost::add_edge( (size_t) localVerts[j+localVertOffset],
968  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
969  }
970  }
971 
972  // Edges
973  for(j = 0; j < nEdges; j++)
974  {
975  if(localEdges[j+localEdgeOffset]==-1)
976  {
977  break;
978  }
979  // Associate to other edges
980  for(k = 0; k < nEdges; k++)
981  {
982  if(localEdges[k+localEdgeOffset]==-1)
983  {
984  break;
985  }
986  if(k!=j)
987  {
988  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
989  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
990  }
991  }
992  // Associate to vertices
993  for(k = 0; k < nVerts; k++)
994  {
995  if(localVerts[k+localVertOffset]==-1)
996  {
997  break;
998  }
999  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1000  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1001  }
1002  // Associate to faces
1003  for(k = 0; k < nFaces; k++)
1004  {
1005  if(localFaces[k+localFaceOffset]==-1)
1006  {
1007  break;
1008  }
1009  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1010  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1011  }
1012  }
1013 
1014  // Faces
1015  for(j = 0; j < nFaces; j++)
1016  {
1017  if(localFaces[j+localFaceOffset]==-1)
1018  {
1019  break;
1020  }
1021  // Associate to other faces
1022  for(k = 0; k < nFaces; k++)
1023  {
1024  if(localFaces[k+localFaceOffset]==-1)
1025  {
1026  break;
1027  }
1028  if(k!=j)
1029  {
1030  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1031  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1032  }
1033  }
1034  // Associate to vertices
1035  for(k = 0; k < nVerts; k++)
1036  {
1037  if(localVerts[k+localVertOffset]==-1)
1038  {
1039  break;
1040  }
1041  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1042  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1043  }
1044  // Associate to edges
1045  for(k = 0; k < nEdges; k++)
1046  {
1047  if(localEdges[k+localEdgeOffset]==-1)
1048  {
1049  break;
1050  }
1051  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1052  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1053  }
1054  }
1055 
1056  localVertOffset+=nVerts;
1057  localEdgeOffset+=nEdges;
1058  localFaceOffset+=nFaces;
1059  }
1060 
1061  // Container to store vertices of the graph which correspond to
1062  // degrees of freedom along the boundary and periodic BCs.
1063  set<int> partVerts;
1064 
1067  {
1068  vector<long> procVerts, procEdges, procFaces;
1069  set <int> foundVerts, foundEdges, foundFaces;
1070 
1071  // Loop over element and construct the procVerts and procEdges
1072  // vectors, which store the geometry IDs of mesh vertices and
1073  // edges respectively which are local to this process.
1074  for(i = cnt = 0; i < locExpVector.size(); ++i)
1075  {
1076  int elmtid = locExp.GetOffset_Elmt_Id(i);
1077  exp = locExpVector[elmtid];
1078  for (j = 0; j < exp->GetNverts(); ++j)
1079  {
1080  int vid = exp->GetGeom()->GetVid(j)+1;
1081  if (foundVerts.count(vid) == 0)
1082  {
1083  procVerts.push_back(vid);
1084  foundVerts.insert(vid);
1085  }
1086  }
1087 
1088  for (j = 0; j < exp->GetNedges(); ++j)
1089  {
1090  int eid = exp->GetGeom()->GetEid(j)+1;
1091 
1092  if (foundEdges.count(eid) == 0)
1093  {
1094  procEdges.push_back(eid);
1095  foundEdges.insert(eid);
1096  }
1097  }
1098 
1099  for (j = 0; j < exp->GetNfaces(); ++j)
1100  {
1101  int fid = exp->GetGeom()->GetFid(j)+1;
1102 
1103  if (foundFaces.count(fid) == 0)
1104  {
1105  procFaces.push_back(fid);
1106  foundFaces.insert(fid);
1107  }
1108  }
1109  }
1110 
1111  int unique_verts = foundVerts.size();
1112  int unique_edges = foundEdges.size();
1113  int unique_faces = foundFaces.size();
1114 
1115  bool verbose = m_session->DefinesCmdLineArgument("verbose");
1116 
1117  // Now construct temporary GS objects. These will be used to
1118  // populate the arrays tmp3 and tmp4 with the multiplicity of
1119  // the vertices and edges respectively to identify those
1120  // vertices and edges which are located on partition boundary.
1121  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1122  Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm, verbose);
1123  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1124  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1125  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1126  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1127  Gs::Finalise(tmp1);
1128 
1129  if (unique_edges > 0)
1130  {
1131  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1132  Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm, verbose);
1133  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1134  Gs::Finalise(tmp2);
1135  }
1136 
1137  if (unique_faces > 0)
1138  {
1139  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1140  Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm, verbose);
1141  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1142  Gs::Finalise(tmp3);
1143  }
1144 
1145  // Finally, fill the partVerts set with all non-Dirichlet
1146  // vertices which lie on a partition boundary.
1147  for (i = 0; i < unique_verts; ++i)
1148  {
1149  if (tmp4[i] > 1.0)
1150  {
1151  if (graph[0].count(procVerts[i]-1) == 0)
1152  {
1153  partVerts.insert(tempGraph[0][procVerts[i]-1]);
1154  }
1155  }
1156  }
1157 
1158  for (i = 0; i < unique_edges; ++i)
1159  {
1160  if (tmp5[i] > 1.0)
1161  {
1162  if (graph[1].count(procEdges[i]-1) == 0)
1163  {
1164  partVerts.insert(tempGraph[1][procEdges[i]-1]);
1165  }
1166  }
1167  }
1168 
1169  for (i = 0; i < unique_faces; ++i)
1170  {
1171  if (tmp6[i] > 1.0)
1172  {
1173  if (graph[2].count(procFaces[i]-1) == 0)
1174  {
1175  partVerts.insert(tempGraph[2][procFaces[i]-1]);
1176  }
1177  }
1178  }
1179 
1180  // Now fill with all vertices on periodic BCs
1181  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
1182  {
1183  if (graph[0].count(pIt->first) == 0)
1184  {
1185  partVerts.insert(tempGraph[0][pIt->first]);
1186  }
1187  }
1188  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1189  {
1190  if (graph[1].count(pIt->first) == 0)
1191  {
1192  partVerts.insert(tempGraph[1][pIt->first]);
1193  }
1194  }
1195  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1196  {
1197  if (graph[2].count(pIt->first) == 0)
1198  {
1199  partVerts.insert(tempGraph[2][pIt->first]);
1200  }
1201  }
1202  }
1203 
1204  int nGraphVerts = tempGraphVertId;
1205  Array<OneD, int> perm(nGraphVerts);
1206  Array<OneD, int> iperm(nGraphVerts);
1207  Array<OneD, int> vwgts(nGraphVerts);
1208  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1209  for(i = 0; i < nGraphVerts; ++i)
1210  {
1211  vwgts[i] = vwgts_map[i];
1212  }
1213 
1214  if(nGraphVerts)
1215  {
1216  switch(m_solnType)
1217  {
1218  case eDirectFullMatrix:
1219  case eIterativeFull:
1220  case eIterativeStaticCond:
1221  case ePETScStaticCond:
1222  case ePETScFullMatrix:
1223  case eXxtFullMatrix:
1224  case eXxtStaticCond:
1225  {
1226  NoReordering(boostGraphObj,perm,iperm);
1227  break;
1228  }
1229 
1230  case eDirectStaticCond:
1231  {
1232  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1233  break;
1234  }
1235 
1240  {
1242  boostGraphObj, perm, iperm, bottomUpGraph,
1243  partVerts, mdswitch);
1244  break;
1245  }
1246  default:
1247  {
1248  ASSERTL0(false,
1249  "Unrecognised solution type: " + std::string(
1251  }
1252  }
1253  }
1254 
1255  // For parallel multi-level static condensation determine the lowest
1256  // static condensation level amongst processors.
1260  m_solnType == eXxtMultiLevelStaticCond) && bottomUpGraph)
1261  {
1262  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1263  vComm->AllReduce(m_lowestStaticCondLevel,
1265  }
1266  else
1267  {
1269  }
1270 
1271  /**
1272  * STEP 4: Fill the #graph[0] and
1273  * #graph[1] with the optimal ordering from boost.
1274  */
1276  for(mapIt = tempGraph[0].begin(); mapIt != tempGraph[0].end(); mapIt++)
1277  {
1278  graph[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1279  }
1280  for(mapIt = tempGraph[1].begin(); mapIt != tempGraph[1].end(); mapIt++)
1281  {
1282  graph[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1283  }
1284  for(mapIt = tempGraph[2].begin(); mapIt != tempGraph[2].end(); mapIt++)
1285  {
1286  graph[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1287  }
1288 
1289  return nGraphVerts;
1290  }
1291 
1292  /**
1293  *
1294  */
1296  const LibUtilities::SessionReaderSharedPtr &pSession,
1297  const int numLocalCoeffs,
1298  const ExpList &locExp,
1299  const BndCondExp &bndCondExp,
1300  const BndCond &bndConditions,
1301  const bool checkIfSystemSingular,
1302  const std::string variable,
1303  const PeriodicMap &periodicVerts,
1304  const PeriodicMap &periodicEdges,
1305  const PeriodicMap &periodicFaces)
1306  : AssemblyMap(pSession, variable)
1307  {
1308  int i, j, k, l;
1309  int p, q, numModes0, numModes1;
1310  int cnt = 0;
1311  int intDofCnt;
1312  int meshVertId, meshEdgeId, meshEdgeId2, meshFaceId, meshFaceId2;
1313  int globalId;
1314  int nEdgeInteriorCoeffs;
1315  int firstNonDirGraphVertId;
1316  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1319  StdRegions::Orientation edgeOrient;
1320  StdRegions::Orientation faceOrient;
1321  Array<OneD, unsigned int> edgeInteriorMap;
1322  Array<OneD, int> edgeInteriorSign;
1323  Array<OneD, unsigned int> faceInteriorMap;
1324  Array<OneD, int> faceInteriorSign;
1325  PeriodicMap::const_iterator pIt;
1326 
1327  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1328 
1329  bool verbose = m_session->DefinesCmdLineArgument("verbose");
1330 
1331  m_signChange = false;
1332 
1333  // Stores vertex, edge and face reordered vertices.
1334  DofGraph graph(3);
1335  DofGraph dofs(3);
1336  vector<map<int, int> > faceModes(2);
1337  map<int, LibUtilities::ShapeType> faceType;
1338 
1339  set<int> extraDirVerts, extraDirEdges;
1341 
1342  // Construct list of number of degrees of freedom for each vertex,
1343  // edge and face.
1344  for (i = 0; i < locExpVector.size(); ++i)
1345  {
1346  exp = locExpVector[i];
1347 
1348  for(j = 0; j < exp->GetNverts(); ++j)
1349  {
1350  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1351  }
1352 
1353  for(j = 0; j < exp->GetNedges(); ++j)
1354  {
1355  if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1356  {
1357  if (dofs[1][exp->GetGeom()->GetEid(j)] !=
1358  exp->GetEdgeNcoeffs(j)-2)
1359  {
1360  ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
1361  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B) ||
1362  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_C),
1363  "CG with variable order only available with modal expansion");
1364  }
1365  dofs[1][exp->GetGeom()->GetEid(j)] =
1366  min(dofs[1][exp->GetGeom()->GetEid(j)],
1367  exp->GetEdgeNcoeffs(j)-2);
1368  }
1369  else
1370  {
1371  dofs[1][exp->GetGeom()->GetEid(j)] =
1372  exp->GetEdgeNcoeffs(j) - 2;
1373  }
1374  }
1375 
1376  for(j = 0; j < exp->GetNfaces(); ++j)
1377  {
1378  faceOrient = exp->GetGeom()->GetForient(j);
1379  meshFaceId = exp->GetGeom()->GetFid(j);
1380  exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
1381 
1382  if (faceModes[0].count(meshFaceId) > 0)
1383  {
1384  faceModes[0][meshFaceId] =
1385  min(faceModes[0][meshFaceId], numModes0);
1386 
1387  faceModes[1][meshFaceId] =
1388  min(faceModes[1][meshFaceId], numModes1);
1389  }
1390  else
1391  {
1392  faceModes[0][meshFaceId] = numModes0;
1393  faceModes[1][meshFaceId] = numModes1;
1394 
1395  // Get shape of this face
1397  geom = boost::dynamic_pointer_cast<SpatialDomains::
1398  Geometry3D> (exp->GetGeom());
1399  faceType[meshFaceId] =
1400  geom->GetFace(j)->GetShapeType();
1401  }
1402  }
1403  }
1404 
1405  // Add non-local periodic dofs to the map
1406  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1407  {
1408  for (i = 0; i < pIt->second.size(); ++i)
1409  {
1410  meshEdgeId2 = pIt->second[i].id;
1411  if (dofs[1].count(meshEdgeId2) == 0)
1412  {
1413  dofs[1][meshEdgeId2] = 1e6;
1414  }
1415  }
1416  }
1417  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1418  {
1419  for (i = 0; i < pIt->second.size(); ++i)
1420  {
1421  meshFaceId2 = pIt->second[i].id;
1422  if (faceModes[0].count(meshFaceId2) == 0)
1423  {
1424  faceModes[0][meshFaceId2] = 1e6;
1425  faceModes[1][meshFaceId2] = 1e6;
1426  }
1427  }
1428  }
1429 
1430  // Now use information from all partitions to determine
1431  // the correct size
1432  map<int, int>::iterator dofIt, dofIt2;
1433  // edges
1434  Array<OneD, long> edgeId (dofs[1].size());
1435  Array<OneD, NekDouble> edgeDof (dofs[1].size());
1436  for(dofIt = dofs[1].begin(), i=0; dofIt != dofs[1].end(); dofIt++, i++)
1437  {
1438  edgeId[i] = dofIt->first + 1;
1439  edgeDof[i] = (NekDouble) dofIt->second;
1440  }
1441  Gs::gs_data *tmp = Gs::Init(edgeId, vComm, verbose);
1442  Gs::Gather(edgeDof, Gs::gs_min, tmp);
1443  Gs::Finalise(tmp);
1444  for (i=0; i < dofs[1].size(); i++)
1445  {
1446  dofs[1][edgeId[i]-1] = (int) (edgeDof[i]+0.5);
1447  }
1448  // Periodic edges
1449  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1450  {
1451  meshEdgeId = pIt->first;
1452  for (i = 0; i < pIt->second.size(); ++i)
1453  {
1454  meshEdgeId2 = pIt->second[i].id;
1455  if (dofs[1][meshEdgeId2] < dofs[1][meshEdgeId])
1456  {
1457  dofs[1][meshEdgeId] = dofs[1][meshEdgeId2];
1458  }
1459  }
1460  }
1461  // faces
1462  Array<OneD, long> faceId (faceModes[0].size());
1463  Array<OneD, NekDouble> faceP (faceModes[0].size());
1464  Array<OneD, NekDouble> faceQ (faceModes[0].size());
1465  for(dofIt = faceModes[0].begin(), dofIt2 = faceModes[1].begin(),i=0;
1466  dofIt != faceModes[0].end(); dofIt++, dofIt2++, i++)
1467  {
1468  faceId[i] = dofIt->first+1;
1469  faceP[i] = (NekDouble) dofIt->second;
1470  faceQ[i] = (NekDouble) dofIt2->second;
1471  }
1472  Gs::gs_data *tmp2 = Gs::Init(faceId, vComm, verbose);
1473  Gs::Gather(faceP, Gs::gs_min, tmp2);
1474  Gs::Gather(faceQ, Gs::gs_min, tmp2);
1475  Gs::Finalise(tmp2);
1476  for (i=0; i < faceModes[0].size(); i++)
1477  {
1478  faceModes[0][faceId[i]-1] = (int) (faceP[i]+0.5);
1479  faceModes[1][faceId[i]-1] = (int) (faceQ[i]+0.5);
1480  }
1481  // Periodic faces
1482  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1483  {
1484  meshFaceId = pIt->first;
1485  for (i = 0; i < pIt->second.size(); ++i)
1486  {
1487  meshFaceId2 = pIt->second[i].id;
1488  if (faceModes[0][meshFaceId2] < faceModes[0][meshFaceId])
1489  {
1490  faceModes[0][meshFaceId] = faceModes[0][meshFaceId2];
1491  }
1492  if (faceModes[1][meshFaceId2] < faceModes[1][meshFaceId])
1493  {
1494  faceModes[1][meshFaceId] = faceModes[1][meshFaceId2];
1495  }
1496  }
1497  }
1498  // Calculate number of dof in each face
1499  int P, Q;
1500  for (i=0; i < faceModes[0].size(); i++)
1501  {
1502  P = faceModes[0][faceId[i]-1];
1503  Q = faceModes[1][faceId[i]-1];
1504  if (faceType[faceId[i]-1] == LibUtilities::eQuadrilateral)
1505  {
1506  // Quad face
1507  dofs[2][faceId[i]-1] =
1510  }
1511  else
1512  {
1513  // Tri face
1514  dofs[2][faceId[i]-1] =
1517  }
1518  }
1519 
1520  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1521 
1522  // Note that nExtraDirichlet is not used in the logic below; it just
1523  // needs to be set so that the coupled solver in
1524  // IncNavierStokesSolver can work.
1525  int nExtraDirichlet;
1526  int mdswitch;
1527  m_session->LoadParameter(
1528  "MDSwitch", mdswitch, 10);
1529  int nGraphVerts =
1530  CreateGraph(locExp, bndCondExp, bndCondVec,
1531  checkIfSystemSingular, periodicVerts, periodicEdges,
1532  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1533  extraDirEdges, firstNonDirGraphVertId,
1534  nExtraDirichlet, mdswitch);
1535 
1536  /*
1537  * Set up an array which contains the offset information of the
1538  * different graph vertices.
1539  *
1540  * This basically means to identify to how many global degrees of
1541  * freedom the individual graph vertices correspond. Obviously,
1542  * the graph vertices corresponding to the mesh-vertices account
1543  * for a single global DOF. However, the graph vertices
1544  * corresponding to the element edges correspond to N-2 global DOF
1545  * where N is equal to the number of boundary modes on this edge.
1546  */
1547  Array<OneD, int> graphVertOffset(
1548  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1549 
1550  graphVertOffset[0] = 0;
1551 
1552  for(i = 0; i < locExpVector.size(); ++i)
1553  {
1554  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1555 
1556  for(j = 0; j < exp->GetNverts(); ++j)
1557  {
1558  meshVertId = exp->GetGeom()->GetVid(j);
1559  graphVertOffset[graph[0][meshVertId]+1] = 1;
1560  }
1561 
1562  for(j = 0; j < exp->GetNedges(); ++j)
1563  {
1564  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1565  meshEdgeId = exp->GetGeom()->GetEid(j);
1566  graphVertOffset[graph[1][meshEdgeId]+1]
1567  = dofs[1][meshEdgeId];
1568 
1569  bType = exp->GetEdgeBasisType(j);
1570 
1571  // need a sign vector for modal expansions if nEdgeCoeffs
1572  // >=3 (not 4 because of variable order case)
1573  if(nEdgeInteriorCoeffs &&
1574  (bType == LibUtilities::eModified_A ||
1575  bType == LibUtilities::eModified_B))
1576  {
1577  m_signChange = true;
1578  }
1579  }
1580 
1581  for(j = 0; j < exp->GetNfaces(); ++j)
1582  {
1583  meshFaceId = exp->GetGeom()->GetFid(j);
1584  graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
1585  }
1586  }
1587 
1588  for(i = 1; i < graphVertOffset.num_elements(); i++)
1589  {
1590  graphVertOffset[i] += graphVertOffset[i-1];
1591  }
1592 
1593  // Allocate the proper amount of space for the class-data
1594  m_numLocalCoeffs = numLocalCoeffs;
1595  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1599  // If required, set up the sign-vector
1600  if(m_signChange)
1601  {
1605  }
1606 
1607  m_staticCondLevel = 0;
1608  m_numPatches = locExpVector.size();
1611  for(i = 0; i < m_numPatches; ++i)
1612  {
1613  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1614  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1615  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1616  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1617  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1618  }
1619 
1620  /**
1621  * STEP 6: Now, all ingredients are ready to set up the actual
1622  * local to global mapping.
1623  *
1624  * The remainder of the map consists of the element-interior
1625  * degrees of freedom. This leads to the block-diagonal submatrix
1626  * as each element-interior mode is globally orthogonal to modes
1627  * in all other elements.
1628  */
1629  cnt = 0;
1630 
1631  // Loop over all the elements in the domain
1632  for(i = 0; i < locExpVector.size(); ++i)
1633  {
1634  exp = locExpVector[i];
1635  cnt = locExp.GetCoeff_Offset(i);
1636  for(j = 0; j < exp->GetNverts(); ++j)
1637  {
1638  meshVertId = exp->GetGeom()->GetVid(j);
1639 
1640  // Set the global DOF for vertex j of element i
1641  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1642  graphVertOffset[graph[0][meshVertId]];
1643  }
1644 
1645  for(j = 0; j < exp->GetNedges(); ++j)
1646  {
1647  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1648  edgeOrient = exp->GetGeom()->GetEorient(j);
1649  meshEdgeId = exp->GetGeom()->GetEid(j);
1650 
1651  pIt = periodicEdges.find(meshEdgeId);
1652 
1653  // See if this edge is periodic. If it is, then we map all
1654  // edges to the one with lowest ID, and align all
1655  // coefficients to this edge orientation.
1656  if (pIt != periodicEdges.end())
1657  {
1658  pair<int, StdRegions::Orientation> idOrient =
1660  meshEdgeId, edgeOrient, pIt->second);
1661  edgeOrient = idOrient.second;
1662  }
1663 
1664  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1665 
1666  // Set the global DOF's for the interior modes of edge j
1667  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1668  {
1669  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1670  graphVertOffset[graph[1][meshEdgeId]]+k;
1671  }
1672  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1673  {
1674  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1675  graphVertOffset[graph[1][meshEdgeId]];
1676  }
1677 
1678  // Fill the sign vector if required
1679  if(m_signChange)
1680  {
1681  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1682  {
1683  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1684  }
1685  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1686  {
1687  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = 0.0;
1688  }
1689  }
1690  }
1691 
1692  for(j = 0; j < exp->GetNfaces(); ++j)
1693  {
1694  faceOrient = exp->GetGeom()->GetForient(j);
1695  meshFaceId = exp->GetGeom()->GetFid(j);
1696 
1697  pIt = periodicFaces.find(meshFaceId);
1698 
1699  if (pIt != periodicFaces.end() &&
1700  meshFaceId == min(meshFaceId, pIt->second[0].id))
1701  {
1702  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1703  }
1704 
1705  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1706 
1707  // Set the global DOF's for the interior modes of face j
1708  exp->GetFaceNumModes(j, faceOrient, numModes0, numModes1);
1709  switch(faceType[meshFaceId])
1710  {
1712  {
1713  int kLoc=0;
1714  k = 0;
1715  for( q = 2; q < numModes1; q++)
1716  {
1717  for( p = 2; p < numModes0; p++)
1718  {
1719  if( (p < faceModes[0][meshFaceId]) &&
1720  (q < faceModes[1][meshFaceId]))
1721  {
1722  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1723  graphVertOffset[graph[2][meshFaceId]]+k;
1724  if(m_signChange)
1725  {
1726  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1727  (NekDouble) faceInteriorSign[kLoc];
1728  }
1729  k++;
1730  }
1731  else
1732  {
1733  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1734  graphVertOffset[graph[2][meshFaceId]];
1735  if(m_signChange)
1736  {
1737  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1738  }
1739  }
1740  kLoc++;
1741  }
1742  }
1743  }
1744  break;
1746  {
1747  int kLoc=0;
1748  k = 0;
1749  for( p = 2; p < numModes0; p++)
1750  {
1751  for( q = 1; q < numModes1-p; q++)
1752  {
1753  if( (p < faceModes[0][meshFaceId]) &&
1754  (p+q < faceModes[1][meshFaceId]))
1755  {
1756  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1757  graphVertOffset[graph[2][meshFaceId]]+k;
1758  if(m_signChange)
1759  {
1760  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] =
1761  (NekDouble) faceInteriorSign[kLoc];
1762  }
1763  k++;
1764  }
1765  else
1766  {
1767  m_localToGlobalMap[cnt+faceInteriorMap[kLoc]] =
1768  graphVertOffset[graph[2][meshFaceId]];
1769  if(m_signChange)
1770  {
1771  m_localToGlobalSign[cnt+faceInteriorMap[kLoc]] = 0.0;
1772  }
1773  }
1774  kLoc++;
1775  }
1776  }
1777  }
1778  break;
1779  default:
1780  ASSERTL0(false,"Shape not recognised");
1781  break;
1782  }
1783  }
1784  }
1785 
1786  // Set up the mapping for the boundary conditions
1787  cnt = 0;
1788  int offset = 0;
1789  for(i = 0; i < bndCondExp.num_elements(); i++)
1790  {
1791  set<int> foundExtraVerts, foundExtraEdges;
1792  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1793  {
1794  bndExp = bndCondExp[i]->GetExp(j);
1795  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1796  for(k = 0; k < bndExp->GetNverts(); k++)
1797  {
1798  meshVertId = bndExp->GetGeom()->GetVid(k);
1799  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1800 
1801  if (bndConditions[i]->GetBoundaryConditionType() !=
1803  {
1804  continue;
1805  }
1806 
1807  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1808  if (iter != extraDirVerts.end() &&
1809  foundExtraVerts.count(meshVertId) == 0)
1810  {
1811  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1812  bndExp->GetVertexMap(k);
1813  int gid = graphVertOffset[
1814  graph[0][meshVertId]];
1815  ExtraDirDof t(loc, gid, 1.0);
1816  m_extraDirDofs[i].push_back(t);
1817  foundExtraVerts.insert(meshVertId);
1818  }
1819  }
1820 
1821  for(k = 0; k < bndExp->GetNedges(); k++)
1822  {
1823  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1824  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1825  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1826 
1827  pIt = periodicEdges.find(meshEdgeId);
1828 
1829  // See if this edge is periodic. If it is, then we map
1830  // all edges to the one with lowest ID, and align all
1831  // coefficients to this edge orientation.
1832  if (pIt != periodicEdges.end())
1833  {
1834  pair<int, StdRegions::Orientation> idOrient =
1836  meshEdgeId, edgeOrient, pIt->second);
1837  edgeOrient = idOrient.second;
1838  }
1839 
1840  bndExp->GetEdgeInteriorMap(
1841  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1842 
1843  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1844  {
1845  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1846  graphVertOffset[graph[1][meshEdgeId]]+l;
1847  }
1848  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1849  {
1850  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1851  graphVertOffset[graph[1][meshEdgeId]];
1852  }
1853 
1854  // Fill the sign vector if required
1855  if(m_signChange)
1856  {
1857  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1858  {
1859  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1860  }
1861  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1862  {
1863  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = 0.0;
1864  }
1865  }
1866 
1867  if (bndConditions[i]->GetBoundaryConditionType() !=
1869  {
1870  continue;
1871  }
1872 
1873  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1874  if (iter != extraDirEdges.end() &&
1875  foundExtraEdges.count(meshEdgeId) == 0 &&
1876  nEdgeInteriorCoeffs > 0)
1877  {
1878  for(l = 0; l < dofs[1][meshEdgeId]; ++l)
1879  {
1880  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1881  edgeInteriorMap[l];
1882  int gid = graphVertOffset[
1883  graph[1][meshEdgeId]]+l;
1884  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1885  m_extraDirDofs[i].push_back(t);
1886  }
1887  for(l = dofs[1][meshEdgeId]; l < nEdgeInteriorCoeffs; ++l)
1888  {
1889  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1890  edgeInteriorMap[l];
1891  int gid = graphVertOffset[
1892  graph[1][meshEdgeId]];
1893  ExtraDirDof t(loc, gid, 0.0);
1894  m_extraDirDofs[i].push_back(t);
1895  }
1896  foundExtraEdges.insert(meshEdgeId);
1897  }
1898  }
1899 
1900  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1901  intDofCnt = 0;
1902  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1903  {
1904  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1905  {
1906  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
1907  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1908  intDofCnt++;
1909  }
1910  }
1911  }
1912  offset += bndCondExp[i]->GetNcoeffs();
1913  }
1914 
1915  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1916  m_numGlobalBndCoeffs = globalId;
1917 
1918 
1919  /*
1920  * The boundary condition mapping is generated from the same vertex
1921  * renumbering.
1922  */
1923  cnt=0;
1924  for(i = 0; i < m_numLocalCoeffs; ++i)
1925  {
1926  if(m_localToGlobalMap[i] == -1)
1927  {
1928  m_localToGlobalMap[i] = globalId++;
1929  }
1930  else
1931  {
1932  if(m_signChange)
1933  {
1935  }
1936  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
1937  }
1938  }
1939 
1940  m_numGlobalCoeffs = globalId;
1941 
1942  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1943 
1944  // Since we can now have multiple entries to m_extraDirDofs due to
1945  // periodic boundary conditions we make a call to work out the
1946  // multiplicity of all entries and invert (Need to be after
1947  // SetupUniversalC0ContMap)
1949 
1950  // Fill in Dirichlet coefficients that are to be sent to other
1951  // processors with a value of 1
1952  map<int, vector<ExtraDirDof> >::iterator Tit;
1953 
1954  // Generate valence for extraDirDofs
1955  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1956  {
1957  for (i = 0; i < Tit->second.size(); ++i)
1958  {
1959  valence[Tit->second[i].get<1>()] = 1.0;
1960  }
1961  }
1962 
1963  UniversalAssembleBnd(valence);
1964 
1965  // Set third argument of tuple to inverse of valence.
1966  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1967  {
1968  for (i = 0; i < Tit->second.size(); ++i)
1969  {
1970  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1971  }
1972  }
1973 
1974  // Set up the local to global map for the next level when using
1975  // multi-level static condensation
1979  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1980  {
1981  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1982  {
1983  Array<OneD, int> vwgts_perm(
1984  graph[0].size() + graph[1].size() + graph[2].size()
1985  - firstNonDirGraphVertId);
1986 
1987  for (i = 0; i < locExpVector.size(); ++i)
1988  {
1989  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1990 
1991  for (j = 0; j < exp->GetNverts(); ++j)
1992  {
1993  meshVertId = exp->GetGeom()->GetVid(j);
1994 
1995  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1996  {
1997  vwgts_perm[graph[0][meshVertId] -
1998  firstNonDirGraphVertId] =
1999  dofs[0][meshVertId];
2000  }
2001  }
2002 
2003  for (j = 0; j < exp->GetNedges(); ++j)
2004  {
2005  meshEdgeId = exp->GetGeom()->GetEid(j);
2006 
2007  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
2008  {
2009  vwgts_perm[graph[1][meshEdgeId] -
2010  firstNonDirGraphVertId] =
2011  dofs[1][meshEdgeId];
2012  }
2013  }
2014 
2015  for (j = 0; j < exp->GetNfaces(); ++j)
2016  {
2017  meshFaceId = exp->GetGeom()->GetFid(j);
2018 
2019  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
2020  {
2021  vwgts_perm[graph[2][meshFaceId] -
2022  firstNonDirGraphVertId] =
2023  dofs[2][meshFaceId];
2024  }
2025  }
2026  }
2027 
2028  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
2030  AllocateSharedPtr(this, bottomUpGraph);
2031  }
2032  }
2033 
2034  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
2035  m_localToGlobalMap.end());
2036 
2037  // Add up hash values if parallel
2038  int hash = m_hash;
2039  vComm->AllReduce(hash, LibUtilities::ReduceSum);
2040  m_hash = hash;
2041 
2044  }
2045 
2046  /**
2047  *
2048  */
2050  {
2053  }
2054 
2055  /**
2056  * @brief Determine orientation of an edge to its periodic equivalents,
2057  * as well as the ID of the representative edge.
2058  *
2059  * Since an edge may be periodic with more than one other edge (e.g. a
2060  * periodic cube has sets of four periodic edges in each coordinate
2061  * direction), we have to define a 'representative' edge. In this
2062  * assembly map we define it to be the one with the minimum ID. This
2063  * routine is set up to calculate the orientation of a given edge with
2064  * ID @p meshEdgeId with respect to the edge ID.
2065  *
2066  * @param meshEdgeId ID of a periodic edge.
2067  * @param edgeOrient Edge orientation of meshEdgeId with respect to
2068  * its parent element.
2069  * @param periodicEdges The map of all periodic edges.
2070  *
2071  * @return Pair containing the ID of the periodic edge and the
2072  * orientation of @p meshEdgeID with respect to this edge.
2073  */
2074  pair<int, StdRegions::Orientation> DeterminePeriodicEdgeOrientId(
2075  int meshEdgeId,
2076  StdRegions::Orientation edgeOrient,
2077  const vector<PeriodicEntity> &periodicEdges)
2078  {
2079  int minId = periodicEdges[0].id;
2080  int minIdK = 0;
2081  int k;
2082 
2083  for (k = 1; k < periodicEdges.size(); ++k)
2084  {
2085  if (periodicEdges[k].id < minId)
2086  {
2087  minId = min(minId, periodicEdges[k].id);
2088  minIdK = k;
2089  }
2090  }
2091 
2092  minId = min(minId, meshEdgeId);
2093 
2094  if (meshEdgeId != minId)
2095  {
2096  if (periodicEdges[minIdK].orient == StdRegions::eBackwards)
2097  {
2098  // Swap edge orientation
2099  edgeOrient = (edgeOrient == StdRegions::eForwards) ?
2101  }
2102  }
2103 
2104  return make_pair(minId, edgeOrient);
2105  }
2106 
2107  /**
2108  * @brief Determine relative orientation between two faces.
2109  *
2110  * Given the orientation of a local element to its local face, defined
2111  * as @p faceOrient, and @p perFaceOrient which states the alignment of
2112  * one periodic face to the other global face, this routine determines
2113  * the orientation that takes this local element face to the
2114  * global/unique face.
2115  *
2116  * @param faceOrient Orientation of the face with respect to its
2117  * parent element.
2118  * @param perFaceOrient Orientation of the representative/global face.
2119  *
2120  * @return Orientation between the two faces.
2121  */
2123  StdRegions::Orientation faceOrient,
2124  StdRegions::Orientation perFaceOrient)
2125  {
2126  StdRegions::Orientation returnval = faceOrient;
2127 
2128  if(perFaceOrient != StdRegions::eDir1FwdDir1_Dir2FwdDir2)
2129  {
2130  int tmp1 = (int)faceOrient - 5;
2131  int tmp2 = (int)perFaceOrient - 5;
2132 
2133  int flipDir1Map [8] = {2,3,0,1,6,7,4,5};
2134  int flipDir2Map [8] = {1,0,3,2,5,4,7,6};
2135  int transposeMap[8] = {4,5,6,7,0,2,1,3};
2136 
2137  // Transpose orientation
2138  if (tmp2 > 3)
2139  {
2140  tmp1 = transposeMap[tmp1];
2141  }
2142 
2143  // Reverse orientation in direction 1.
2144  if (tmp2 == 2 || tmp2 == 3 || tmp2 == 6 || tmp2 == 7)
2145  {
2146  tmp1 = flipDir1Map[tmp1];
2147  }
2148 
2149  // Reverse orientation in direction 2
2150  if (tmp2 % 2 == 1)
2151  {
2152  tmp1 = flipDir2Map[tmp1];
2153  }
2154 
2155  returnval = (StdRegions::Orientation)(tmp1+5);
2156  }
2157  return returnval;
2158  }
2159 
2160 
2161  /**
2162  * Sets up the global to universal mapping of degrees of freedom across
2163  * processors.
2164  */
2166  const ExpList &locExp,
2167  const PeriodicMap &perVerts,
2168  const PeriodicMap &perEdges,
2169  const PeriodicMap &perFaces)
2170  {
2172  int nVert = 0;
2173  int nEdge = 0;
2174  int nFace = 0;
2175  int maxEdgeDof = 0;
2176  int maxFaceDof = 0;
2177  int maxIntDof = 0;
2178  int dof = 0;
2179  int cnt;
2180  int i,j,k,l;
2181  int meshVertId;
2182  int meshEdgeId;
2183  int meshFaceId;
2184  int elementId;
2185  int vGlobalId;
2186  int maxBndGlobalId = 0;
2187  StdRegions::Orientation edgeOrient;
2188  StdRegions::Orientation faceOrient;
2189  Array<OneD, unsigned int> edgeInteriorMap;
2190  Array<OneD, int> edgeInteriorSign;
2191  Array<OneD, unsigned int> faceInteriorMap;
2192  Array<OneD, int> faceInteriorSign;
2193  Array<OneD, unsigned int> interiorMap;
2194  PeriodicMap::const_iterator pIt;
2195 
2196  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
2197  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
2198  const bool verbose = locExp.GetSession()->DefinesCmdLineArgument("verbose");
2199 
2204 
2205  // Loop over all the elements in the domain to gather mesh data
2206  for(i = 0; i < locExpVector.size(); ++i)
2207  {
2208  exp = locExpVector[i];
2209  nVert += exp->GetNverts();
2210  nEdge += exp->GetNedges();
2211  nFace += exp->GetNfaces();
2212  // Loop over all edges (and vertices) of element i
2213  for(j = 0; j < exp->GetNedges(); ++j)
2214  {
2215  dof = exp->GetEdgeNcoeffs(j)-2;
2216  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2217  }
2218  for(j = 0; j < exp->GetNfaces(); ++j)
2219  {
2220  dof = exp->GetFaceIntNcoeffs(j);
2221  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2222  }
2223  exp->GetInteriorMap(interiorMap);
2224  dof = interiorMap.num_elements();
2225  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2226  }
2227 
2228  // Tell other processes about how many dof we have
2229  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2230  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2231  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2232  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2233  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2234  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2235 
2236  // Assemble global to universal mapping for this process
2237  for(i = 0; i < locExpVector.size(); ++i)
2238  {
2239  exp = locExpVector[i];
2240  cnt = locExp.GetCoeff_Offset(i);
2241 
2242  // Loop over all vertices of element i
2243  for(j = 0; j < exp->GetNverts(); ++j)
2244  {
2245  meshVertId = exp->GetGeom()->GetVid(j);
2246  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2247 
2248  pIt = perVerts.find(meshVertId);
2249  if (pIt != perVerts.end())
2250  {
2251  for (k = 0; k < pIt->second.size(); ++k)
2252  {
2253  meshVertId = min(meshVertId, pIt->second[k].id);
2254  }
2255  }
2256 
2257  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2258  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2259  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2260  }
2261 
2262  // Loop over all edges of element i
2263  for(j = 0; j < exp->GetNedges(); ++j)
2264  {
2265  meshEdgeId = exp->GetGeom()->GetEid(j);
2266  pIt = perEdges.find(meshEdgeId);
2267  edgeOrient = exp->GetGeom()->GetEorient(j);
2268 
2269  if (pIt != perEdges.end())
2270  {
2271  pair<int, StdRegions::Orientation> idOrient =
2273  meshEdgeId, edgeOrient, pIt->second);
2274  meshEdgeId = idOrient.first;
2275  edgeOrient = idOrient.second;
2276  }
2277 
2278  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
2279  dof = exp->GetEdgeNcoeffs(j)-2;
2280 
2281  // Set the global DOF's for the interior modes of edge j
2282  // for varP, ignore modes with sign == 0
2283  for(k = 0, l = 0; k < dof; ++k)
2284  {
2285  if (m_signChange)
2286  {
2287  if (m_localToGlobalSign[cnt+edgeInteriorMap[k]]==0)
2288  {
2289  continue;
2290  }
2291  }
2292  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2293  m_globalToUniversalMap[vGlobalId]
2294  = nVert + meshEdgeId * maxEdgeDof + l + 1;
2295  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2296  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2297  l++;
2298  }
2299  }
2300 
2301  // Loop over all faces of element i
2302  for(j = 0; j < exp->GetNfaces(); ++j)
2303  {
2304  faceOrient = exp->GetGeom()->GetForient(j);
2305 
2306  meshFaceId = exp->GetGeom()->GetFid(j);
2307 
2308  pIt = perFaces.find(meshFaceId);
2309  if (pIt != perFaces.end())
2310  {
2311  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2312  {
2313  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
2314  }
2315  meshFaceId = min(meshFaceId, pIt->second[0].id);
2316  }
2317 
2318 
2319  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2320  dof = exp->GetFaceIntNcoeffs(j);
2321 
2322  for(k = 0, l = 0; k < dof; ++k)
2323  {
2324  if (m_signChange)
2325  {
2326  if (m_localToGlobalSign[cnt+faceInteriorMap[k]]==0)
2327  {
2328  continue;
2329  }
2330  }
2331  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2332  m_globalToUniversalMap[vGlobalId]
2333  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2334  + l + 1;
2335  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2336 
2337  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2338  l++;
2339  }
2340  }
2341 
2342  // Add interior DOFs to complete universal numbering
2343  exp->GetInteriorMap(interiorMap);
2344  dof = interiorMap.num_elements();
2345  elementId = (exp->GetGeom())->GetGlobalID();
2346  for (k = 0; k < dof; ++k)
2347  {
2348  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2349  m_globalToUniversalMap[vGlobalId]
2350  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2351  }
2352  }
2353 
2354  // Set up the GSLib universal assemble mapping
2355  // Internal DOF do not participate in any data
2356  // exchange, so we keep these set to the special GSLib id=0 so
2357  // they are ignored.
2358  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2359  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2360  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2361  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2362  {
2363  tmp[i] = m_globalToUniversalMap[i];
2364  }
2365 
2366  m_gsh = Gs::Init(tmp, vCommRow, verbose);
2367  m_bndGsh = Gs::Init(tmp2, vCommRow, verbose);
2368  Gs::Unique(tmp, vCommRow);
2369  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2370  {
2371  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2372  }
2373  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2374  {
2375  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2376  }
2377  }
2378 
2379  /**
2380  * @brief Construct an AssemblyMapCG object which corresponds to the
2381  * linear space of the current object.
2382  *
2383  * This function is used to create a linear-space assembly map, which is
2384  * then used in the linear space preconditioner in the conjugate
2385  * gradient solve.
2386  */
2388  const ExpList &locexp, GlobalSysSolnType solnType)
2389  {
2390  AssemblyMapCGSharedPtr returnval;
2391 
2392  int i, j;
2393  int nverts = 0;
2394  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2395  = locexp.GetExp();
2396  int nelmts = exp->size();
2397  const bool verbose = locexp.GetSession()->DefinesCmdLineArgument("verbose");
2398 
2399  // Get Default Map and turn off any searched values.
2400  returnval = MemoryManager<AssemblyMapCG>
2402  returnval->m_solnType = solnType;
2403  returnval->m_preconType = eNull;
2404  returnval->m_maxStaticCondLevel = 0;
2405  returnval->m_signChange = false;
2406  returnval->m_comm = m_comm;
2407 
2408  // Count the number of vertices
2409  for (i = 0; i < nelmts; ++i)
2410  {
2411  nverts += (*exp)[i]->GetNverts();
2412  }
2413 
2414  returnval->m_numLocalCoeffs = nverts;
2415  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2416 
2417  // Store original global ids in this map
2418  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2419 
2420  int cnt = 0;
2421  int cnt1 = 0;
2422  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2423 
2424  // Set up local to global map;
2425  for (i = 0; i < nelmts; ++i)
2426  {
2427  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2428  {
2429  returnval->m_localToGlobalMap[cnt] =
2430  returnval->m_localToGlobalBndMap[cnt] =
2431  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2432  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2433 
2434  // Set up numLocalDirBndCoeffs
2435  if ((returnval->m_localToGlobalMap[cnt]) <
2437  {
2438  returnval->m_numLocalDirBndCoeffs++;
2439  }
2440  cnt++;
2441  }
2442  cnt1 += (*exp)[i]->GetNcoeffs();
2443  }
2444 
2445  cnt = 0;
2446  // Reset global numbering and count number of dofs
2447  for (i = 0; i < m_numGlobalCoeffs; ++i)
2448  {
2449  if (GlobCoeffs[i] != -1)
2450  {
2451  GlobCoeffs[i] = cnt++;
2452  }
2453  }
2454 
2455  // Set up number of globalCoeffs;
2456  returnval->m_numGlobalCoeffs = cnt;
2457 
2458  // Set up number of global Dirichlet boundary coefficients
2459  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2460  {
2461  if (GlobCoeffs[i] != -1)
2462  {
2463  returnval->m_numGlobalDirBndCoeffs++;
2464  }
2465  }
2466 
2467  // Set up global to universal map
2468  if (m_globalToUniversalMap.num_elements())
2469  {
2471  = m_session->GetComm()->GetRowComm();
2472  int nglocoeffs = returnval->m_numGlobalCoeffs;
2473  returnval->m_globalToUniversalMap
2474  = Array<OneD, int> (nglocoeffs);
2475  returnval->m_globalToUniversalMapUnique
2476  = Array<OneD, int> (nglocoeffs);
2477 
2478  // Reset local to global map and setup universal map
2479  for (i = 0; i < nverts; ++i)
2480  {
2481  cnt = returnval->m_localToGlobalMap[i];
2482  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2483 
2484  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2486  }
2487 
2488  Nektar::Array<OneD, long> tmp(nglocoeffs);
2489  Vmath::Zero(nglocoeffs, tmp, 1);
2490  for (unsigned int i = 0; i < nglocoeffs; ++i)
2491  {
2492  tmp[i] = returnval->m_globalToUniversalMap[i];
2493  }
2494  returnval->m_gsh = Gs::Init(tmp, vCommRow, verbose);
2495  Gs::Unique(tmp, vCommRow);
2496  for (unsigned int i = 0; i < nglocoeffs; ++i)
2497  {
2498  returnval->m_globalToUniversalMapUnique[i]
2499  = (tmp[i] >= 0 ? 1 : 0);
2500  }
2501  }
2502  else // not sure this option is ever needed.
2503  {
2504  for (i = 0; i < nverts; ++i)
2505  {
2506  cnt = returnval->m_localToGlobalMap[i];
2507  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2508  }
2509  }
2510 
2511  return returnval;
2512  }
2513 
2514  /**
2515  * The bandwidth calculated here corresponds to what is referred to as
2516  * half-bandwidth. If the elements of the matrix are designated as
2517  * a_ij, it corresponds to the maximum value of |i-j| for non-zero
2518  * a_ij. As a result, the value also corresponds to the number of
2519  * sub- or super-diagonals.
2520  *
2521  * The bandwith can be calculated elementally as it corresponds to the
2522  * maximal elemental bandwith (i.e. the maximal difference in global
2523  * DOF index for every element).
2524  *
2525  * We caluclate here the bandwith of the full global system.
2526  */
2528  {
2529  int i,j;
2530  int cnt = 0;
2531  int locSize;
2532  int maxId;
2533  int minId;
2534  int bwidth = -1;
2535  for(i = 0; i < m_numPatches; ++i)
2536  {
2538  maxId = -1;
2539  minId = m_numLocalCoeffs+1;
2540  for(j = 0; j < locSize; j++)
2541  {
2543  {
2544  if(m_localToGlobalMap[cnt+j] > maxId)
2545  {
2546  maxId = m_localToGlobalMap[cnt+j];
2547  }
2548 
2549  if(m_localToGlobalMap[cnt+j] < minId)
2550  {
2551  minId = m_localToGlobalMap[cnt+j];
2552  }
2553  }
2554  }
2555  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2556 
2557  cnt+=locSize;
2558  }
2559 
2560  m_fullSystemBandWidth = bwidth;
2561  }
2562 
2563 
2565  {
2566  return m_localToGlobalMap[i];
2567  }
2568 
2570  {
2571  return m_globalToUniversalMap[i];
2572  }
2573 
2575  {
2576  return m_globalToUniversalMapUnique[i];
2577  }
2578 
2579  const Array<OneD,const int>&
2581  {
2582  return m_localToGlobalMap;
2583  }
2584 
2585  const Array<OneD,const int>&
2587  {
2588  return m_globalToUniversalMap;
2589  }
2590 
2591  const Array<OneD,const int>&
2593  {
2595  }
2596 
2598  const int i) const
2599  {
2600  if(m_signChange)
2601  {
2602  return m_localToGlobalSign[i];
2603  }
2604  else
2605  {
2606  return 1.0;
2607  }
2608  }
2609 
2611  {
2612  return m_localToGlobalSign;
2613  }
2614 
2616  const Array<OneD, const NekDouble>& loc,
2617  Array<OneD, NekDouble>& global,
2618  bool useComm) const
2619  {
2621  if(global.data() == loc.data())
2622  {
2623  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2624  }
2625  else
2626  {
2627  local = loc; // create reference
2628  }
2629 
2630 
2631  if(m_signChange)
2632  {
2633  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2634  }
2635  else
2636  {
2637  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2638  }
2639 
2640  // ensure all values are unique by calling a max
2641  if(useComm)
2642  {
2643  Gs::Gather(global, Gs::gs_max, m_gsh);
2644  }
2645  }
2646 
2648  const NekVector<NekDouble>& loc,
2649  NekVector< NekDouble>& global,
2650  bool useComm) const
2651  {
2652  LocalToGlobal(loc.GetPtr(),global.GetPtr(),useComm);
2653  }
2654 
2656  const Array<OneD, const NekDouble>& global,
2657  Array<OneD, NekDouble>& loc) const
2658  {
2660  if(global.data() == loc.data())
2661  {
2662  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2663  }
2664  else
2665  {
2666  glo = global; // create reference
2667  }
2668 
2669 
2670  if(m_signChange)
2671  {
2672  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2673  }
2674  else
2675  {
2676  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2677  }
2678  }
2679 
2681  const NekVector<NekDouble>& global,
2682  NekVector< NekDouble>& loc) const
2683  {
2684  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2685  }
2686 
2688  const Array<OneD, const NekDouble> &loc,
2689  Array<OneD, NekDouble> &global) const
2690  {
2692  if(global.data() == loc.data())
2693  {
2694  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2695  }
2696  else
2697  {
2698  local = loc; // create reference
2699  }
2700  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2701 
2702  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2703 
2704  if(m_signChange)
2705  {
2706  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2707  }
2708  else
2709  {
2710  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2711  }
2712  UniversalAssemble(global);
2713  }
2714 
2716  const NekVector<NekDouble>& loc,
2717  NekVector< NekDouble>& global) const
2718  {
2719  Assemble(loc.GetPtr(),global.GetPtr());
2720  }
2721 
2723  Array<OneD, NekDouble>& pGlobal) const
2724  {
2725  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2726  }
2727 
2729  NekVector< NekDouble>& pGlobal) const
2730  {
2731  UniversalAssemble(pGlobal.GetPtr());
2732  }
2733 
2735  Array<OneD, NekDouble>& pGlobal,
2736  int offset) const
2737  {
2738  Array<OneD, NekDouble> tmp(offset);
2739  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2740  UniversalAssemble(pGlobal);
2741  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2742  }
2743 
2745  {
2746  return m_fullSystemBandWidth;
2747  }
2748 
2750  {
2751  return m_numNonDirVertexModes;
2752  }
2753 
2755  {
2756  return m_numNonDirEdgeModes;
2757  }
2758 
2760  {
2761  return m_numNonDirFaceModes;
2762  }
2763 
2765  {
2766  return m_numDirEdges;
2767  }
2768 
2770  {
2771  return m_numDirFaces;
2772  }
2773 
2775  {
2776  return m_numNonDirEdges;
2777  }
2778 
2780  {
2781  return m_numNonDirFaces;
2782  }
2783 
2785  {
2786  return m_extraDirEdges;
2787  }
2788  } // namespace
2789 } // namespace
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool m_systemSingular
Flag indicating if the system is singular or not.
Definition: AssemblyMap.h:322
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
AssemblyMapCG(const LibUtilities::SessionReaderSharedPtr &pSession, const std::string variable="DefaultVar")
Default constructor.
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:120
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:316
Principle Modified Functions .
Definition: BasisType.h:51
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:2084
void MultiLevelBisectionReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm, BottomUpSubStructuredGraphSharedPtr &substructgraph, std::set< int > partVerts, int mdswitch)
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: AssemblyMap.h:308
static void Gather(Nektar::Array< OneD, NekDouble > pU, gs_op pOp, gs_data *pGsh, Nektar::Array< OneD, NekDouble > pBuffer=NullNekDouble1DArray)
Performs a gather-scatter operation of the provided values.
Definition: GsLib.hpp:239
int m_maxStaticCondLevel
Maximum static condensation level.
void Gathr(int n, const T *x, const int *y, T *z)
Gather vector z[i] = x[y[i]].
Definition: Vmath.cpp:644
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
static void Finalise(gs_data *pGsh)
Deallocates the GSLib mapping data.
Definition: GsLib.hpp:222
boost::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:53
virtual int v_GetNumNonDirEdgeModes() const
Array< OneD, int > m_globalToUniversalMapUnique
Integer map of unique process coeffs to universal space (signed)
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:779
Array< OneD, int > m_extraDirEdges
Extra dirichlet edges in parallel.
void SetUpUniversalC0ContMap(const ExpList &locExp, const PeriodicMap &perVerts=NullPeriodicMap, const PeriodicMap &perEdges=NullPeriodicMap, const PeriodicMap &perFaces=NullPeriodicMap)
int CreateGraph(const ExpList &locExp, const BndCondExp &bndCondExp, const Array< OneD, const BndCond > &bndConditions, const bool checkIfSystemSingular, const PeriodicMap &periodicVerts, const PeriodicMap &periodicEdges, const PeriodicMap &periodicFaces, DofGraph &graph, BottomUpSubStructuredGraphSharedPtr &bottomUpGraph, std::set< int > &extraDirVerts, std::set< int > &extraDirEdges, int &firstNonDirGraphVertId, int &nExtraDirichlet, int mdswitch=1)
virtual int v_GetNumNonDirFaceModes() const
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2075
int m_numNonDirVertexModes
Number of non Dirichlet vertex modes.
Principle Modified Functions .
Definition: BasisType.h:49
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
STL namespace.
Array< OneD, int > m_globalToUniversalMap
Integer map of process coeffs to universal space.
int m_numNonDirEdges
Number of Dirichlet edges.
int m_numLocalCoeffs
Total number of local coefficients.
Definition: AssemblyMap.h:333
void CuthillMckeeReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
int m_numLocalBndCondCoeffs
Number of local boundary condition coefficients.
virtual void v_Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
virtual void v_GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
std::map< int, std::vector< ExtraDirDof > > m_extraDirDofs
Map indicating degrees of freedom which are Dirichlet but whose value is stored on another processor...
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
int Imax(int n, const T *x, const int incx)
Return the index of the maximum element in x.
Definition: Vmath.cpp:755
void CalculateFullSystemBandWidth()
Calculate the bandwith of the full matrix system.
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Array< OneD, int > m_localToGlobalMap
Integer map of local coeffs to global space.
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap()
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:397
Base class for constructing local to global mapping of degrees of freedom.
Definition: AssemblyMap.h:59
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
const char *const GlobalSysSolnTypeMap[]
int m_numNonDirFaceModes
Number of non Dirichlet face modes.
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:311
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
void LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm=true) const
Array< OneD, unsigned int > m_numLocalBndCoeffsPerPatch
The number of bnd dofs per patch.
Definition: AssemblyMap.h:390
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:365
Array< OneD, NekDouble > m_bndCondCoeffsToGlobalCoeffsSign
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:356
boost::tuple< int, int, NekDouble > ExtraDirDof
Definition: AssemblyMapCG.h:54
virtual AssemblyMapSharedPtr v_LinearSpaceMap(const ExpList &locexp, GlobalSysSolnType solnType)
Construct an AssemblyMapCG object which corresponds to the linear space of the current object...
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
Principle Modified Functions .
Definition: BasisType.h:50
StdRegions::Orientation DeterminePeriodicFaceOrient(StdRegions::Orientation faceOrient, StdRegions::Orientation perFaceOrient)
Determine relative orientation between two faces.
void GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
int m_numDirFaces
Number of Dirichlet faces.
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:673
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:111
static void Unique(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm)
Updates pId to negate all-but-one references to each universal ID.
Definition: GsLib.hpp:200
boost::shared_ptr< LibUtilities::SessionReader > GetSession() const
Returns the session object.
Definition: ExpList.h:924
virtual int v_GetNumNonDirFaces() const
double NekDouble
void Assmb(int n, const T *x, const int *y, T *z)
Assemble z[y[i]] += x[i]; z should be zero'd first.
Definition: Vmath.cpp:709
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
virtual int v_GetFullSystemBandWidth() const
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:392
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:399
void CalculateBndSystemBandWidth()
Calculates the bandwidth of the boundary system.
Array< OneD, int > m_localToGlobalBndMap
Integer map of local boundary coeffs to global space.
Definition: AssemblyMap.h:350
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
No Solution type specified.
virtual const Array< OneD, const int > & v_GetExtraDirEdges()
int getNumberOfBndCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:139
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:354
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
std::vector< std::map< int, int > > DofGraph
Definition: AssemblyMapCG.h:56
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:314
int m_numDirEdges
Number of Dirichlet edges.
3D geometry information
Definition: Geometry3D.h:70
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:386
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:352
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:161
Array< OneD, NekDouble > m_localToGlobalSign
Integer sign of local coeffs to global space.
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:132
int GetOffset_Elmt_Id(int n) const
Get the element id associated with the n th consecutive block of data in m_phys and m_coeffs...
Definition: ExpList.h:2100
int m_numNonDirEdgeModes
Number of non Dirichlet edge modes.
virtual int v_GetNumNonDirVertexModes() const
LibUtilities::SessionReaderSharedPtr m_session
Session object.
Definition: AssemblyMap.h:305
Array< OneD, int > m_globalToUniversalBndMap
Integer map of process coeffs to universal space.
Definition: AssemblyMap.h:360
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:362
void UniversalAssembleBnd(Array< OneD, NekDouble > &pGlobal) const
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
pair< int, StdRegions::Orientation > DeterminePeriodicEdgeOrientId(int meshEdgeId, StdRegions::Orientation edgeOrient, const vector< PeriodicEntity > &periodicEdges)
Determine orientation of an edge to its periodic equivalents, as well as the ID of the representative...
int m_numNonDirFaces
Number of Dirichlet faces.
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign() const
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique()
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:737
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
virtual int v_GetNumNonDirEdges() const
boost::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:52
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap()
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
int m_fullSystemBandWidth
Bandwith of the full matrix system (no static condensation).
void Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:388
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
boost::shared_ptr< Geometry3D > Geometry3DSharedPtr
Definition: Geometry3D.h:52