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