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  gId = gId < 0 ? graphVertId++ : gId;
502  graph[0][meshVertId] = gId;
503 
504  for (i = 0; i < pIt->second.size(); ++i)
505  {
506  if (pIt->second[i].isLocal)
507  {
508  graph[0][pIt->second[i].id] = graph[0][meshVertId];
509  }
510  }
511  }
512  else
513  {
514  bool found = false;
515  for (i = 0; i < pIt->second.size(); ++i)
516  {
517  if (pIt->second[i].id == meshVertId)
518  {
519  found = true;
520  break;
521  }
522  }
523 
524  if (found)
525  {
526  gId = gId < 0 ? graphVertId++ : gId;
527  graph[0][pIt->first] = gId;
528 
529  for (i = 0; i < pIt->second.size(); ++i)
530  {
531  if (pIt->second[i].isLocal)
532  {
533  graph[0][pIt->second[i].id] = graph[0][pIt->first];
534  }
535  }
536  }
537  }
538  }
539  }
540 
541  // Add extra dirichlet boundary conditions to count.
542  m_numLocalDirBndCoeffs += nExtraDirichlet;
543  firstNonDirGraphVertId = graphVertId;
544 
545  typedef boost::adjacency_list<
546  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
547  BoostGraph boostGraphObj;
548 
549  vector<map<int,int> > tempGraph(3);
550  map<int, int> vwgts_map;
551  Array<OneD, int> localVerts;
552  Array<OneD, int> localEdges;
553  Array<OneD, int> localFaces;
554 
555  int tempGraphVertId = 0;
556  int localVertOffset = 0;
557  int localEdgeOffset = 0;
558  int localFaceOffset = 0;
559  int nTotalVerts = 0;
560  int nTotalEdges = 0;
561  int nTotalFaces = 0;
562  int nVerts;
563  int nEdges;
564  int nFaces;
565  int vertCnt;
566  int edgeCnt;
567  int faceCnt;
568 
570  m_numNonDirEdges = 0;
571  m_numNonDirFaces = 0;
575 
576  map<int,int> EdgeSize;
577  map<int,int> FaceSize;
578 
579  /// - Count verts, edges, face and add up edges and face sizes
580  for(i = 0; i < locExpVector.size(); ++i)
581  {
582  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
583  nTotalVerts += exp->GetNverts();
584  nTotalEdges += exp->GetNedges();
585  nTotalFaces += exp->GetNfaces();
586 
587  nEdges = exp->GetNedges();
588  for(j = 0; j < nEdges; ++j)
589  {
590  meshEdgeId = exp->GetGeom()->GetEid(j);
591  if (EdgeSize.count(meshEdgeId) > 0)
592  {
593  EdgeSize[meshEdgeId] =
594  min(EdgeSize[meshEdgeId],
595  exp->GetEdgeNcoeffs(j) - 2);
596  }
597  else
598  {
599  EdgeSize[meshEdgeId] = exp->GetEdgeNcoeffs(j) - 2;
600  }
601  }
602 
603  nFaces = exp->GetNfaces();
604  faceCnt = 0;
605  for(j = 0; j < nFaces; ++j)
606  {
607  meshFaceId = exp->GetGeom()->GetFid(j);
608  if (FaceSize.count(meshFaceId) > 0)
609  {
610  FaceSize[meshFaceId] =
611  min(FaceSize[meshFaceId],
612  exp->GetFaceIntNcoeffs(j));
613  }
614  else
615  {
616  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
617  }
618  FaceSize[meshFaceId] = exp->GetFaceIntNcoeffs(j);
619  }
620  }
621 
622  /// - Periodic vertices
623  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
624  {
625  meshVertId = pIt->first;
626 
627  // This periodic vertex is joined to a Dirichlet condition.
628  if (graph[0].count(pIt->first) != 0)
629  {
630  for (i = 0; i < pIt->second.size(); ++i)
631  {
632  meshVertId2 = pIt->second[i].id;
633  if (graph[0].count(meshVertId2) == 0 &&
634  pIt->second[i].isLocal)
635  {
636  graph[0][meshVertId2] =
637  graph[0][meshVertId];
638  }
639  }
640  continue;
641  }
642 
643  // One of the attached vertices is Dirichlet.
644  bool isDirichlet = false;
645  for (i = 0; i < pIt->second.size(); ++i)
646  {
647  if (!pIt->second[i].isLocal)
648  {
649  continue;
650  }
651 
652  meshVertId2 = pIt->second[i].id;
653  if (graph[0].count(meshVertId2) > 0)
654  {
655  isDirichlet = true;
656  break;
657  }
658  }
659 
660  if (isDirichlet)
661  {
662  graph[0][meshVertId] =
663  graph[0][pIt->second[i].id];
664 
665  for (j = 0; j < pIt->second.size(); ++j)
666  {
667  meshVertId2 = pIt->second[i].id;
668  if (j == i || !pIt->second[j].isLocal ||
669  graph[0].count(meshVertId2) > 0)
670  {
671  continue;
672  }
673 
674  graph[0][meshVertId2] =
675  graph[0][pIt->second[i].id];
676  }
677 
678  continue;
679  }
680 
681  // Otherwise, see if a vertex ID has already been set.
682  for (i = 0; i < pIt->second.size(); ++i)
683  {
684  if (!pIt->second[i].isLocal)
685  {
686  continue;
687  }
688 
689  if (tempGraph[0].count(pIt->second[i].id) > 0)
690  {
691  break;
692  }
693  }
694 
695  if (i == pIt->second.size())
696  {
697  boost::add_vertex(boostGraphObj);
698  tempGraph[0][meshVertId] = tempGraphVertId++;
700  }
701  else
702  {
703  tempGraph[0][meshVertId] = tempGraph[0][pIt->second[i].id];
704  }
705  }
706 
707  // Store the temporary graph vertex id's of all element edges and
708  // vertices in these 3 arrays below
709  localVerts = Array<OneD, int>(nTotalVerts,-1);
710  localEdges = Array<OneD, int>(nTotalEdges,-1);
711  localFaces = Array<OneD, int>(nTotalFaces,-1);
712 
713  // Set up vertex numbering
714  for(i = 0; i < locExpVector.size(); ++i)
715  {
716  exp = locExpVector[i];
717  vertCnt = 0;
718  nVerts = exp->GetNverts();
719  for(j = 0; j < nVerts; ++j)
720  {
721  meshVertId = exp->GetGeom()->GetVid(j);
722  if(graph[0].count(meshVertId) == 0)
723  {
724  if(tempGraph[0].count(meshVertId) == 0)
725  {
726  boost::add_vertex(boostGraphObj);
727  tempGraph[0][meshVertId] = tempGraphVertId++;
729  }
730  localVerts[localVertOffset+vertCnt++] = tempGraph[0][meshVertId];
731  vwgts_map[ tempGraph[0][meshVertId] ] = 1;
732  }
733  }
734 
735  localVertOffset+=nVerts;
736  }
737 
738  /// - Periodic edges
739  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
740  {
741  meshEdgeId = pIt->first;
742 
743  // This periodic edge is joined to a Dirichlet condition.
744  if (graph[1].count(pIt->first) != 0)
745  {
746  for (i = 0; i < pIt->second.size(); ++i)
747  {
748  meshEdgeId2 = pIt->second[i].id;
749  if (graph[1].count(meshEdgeId2) == 0 &&
750  pIt->second[i].isLocal)
751  {
752  graph[1][meshEdgeId2] =
753  graph[1][meshEdgeId];
754  }
755  }
756  continue;
757  }
758 
759  // One of the attached edges is Dirichlet.
760  bool isDirichlet = false;
761  for (i = 0; i < pIt->second.size(); ++i)
762  {
763  if (!pIt->second[i].isLocal)
764  {
765  continue;
766  }
767 
768  meshEdgeId2 = pIt->second[i].id;
769  if (graph[1].count(meshEdgeId2) > 0)
770  {
771  isDirichlet = true;
772  break;
773  }
774  }
775 
776  if (isDirichlet)
777  {
778  graph[1][meshEdgeId] =
779  graph[1][pIt->second[i].id];
780 
781  for (j = 0; j < pIt->second.size(); ++j)
782  {
783  meshEdgeId2 = pIt->second[i].id;
784  if (j == i || !pIt->second[j].isLocal ||
785  graph[1].count(meshEdgeId2) > 0)
786  {
787  continue;
788  }
789 
790  graph[1][meshEdgeId2] =
791  graph[1][pIt->second[i].id];
792  }
793 
794  continue;
795  }
796 
797  // Otherwise, see if a edge ID has already been set.
798  for (i = 0; i < pIt->second.size(); ++i)
799  {
800  if (!pIt->second[i].isLocal)
801  {
802  continue;
803  }
804 
805  if (tempGraph[1].count(pIt->second[i].id) > 0)
806  {
807  break;
808  }
809  }
810 
811  if (i == pIt->second.size())
812  {
813  boost::add_vertex(boostGraphObj);
814  tempGraph[1][meshEdgeId] = tempGraphVertId++;
815  m_numNonDirEdgeModes += EdgeSize[meshEdgeId];
817  }
818  else
819  {
820  tempGraph[1][meshEdgeId] = tempGraph[1][pIt->second[i].id];
821  }
822  }
823 
824  int nEdgeIntCoeffs, nFaceIntCoeffs;
825 
826  // Set up edge numbering
827  for(i = 0; i < locExpVector.size(); ++i)
828  {
829  exp = locExpVector[i];
830  edgeCnt = 0;
831  nEdges = exp->GetNedges();
832 
833  for(j = 0; j < nEdges; ++j)
834  {
835  nEdgeIntCoeffs = exp->GetEdgeNcoeffs(j) - 2;
836  meshEdgeId = exp->GetGeom()->GetEid(j);
837  if(graph[1].count(meshEdgeId) == 0)
838  {
839  if(tempGraph[1].count(meshEdgeId) == 0)
840  {
841  boost::add_vertex(boostGraphObj);
842  tempGraph[1][meshEdgeId] = tempGraphVertId++;
843  m_numNonDirEdgeModes+=nEdgeIntCoeffs;
844 
846  }
847  localEdges[localEdgeOffset+edgeCnt++] = tempGraph[1][meshEdgeId];
848  vwgts_map[ tempGraph[1][meshEdgeId] ] = nEdgeIntCoeffs;
849  }
850  }
851 
852  localEdgeOffset+=nEdges;
853  }
854 
855  /// - Periodic faces
856  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
857  {
858  if (!pIt->second[0].isLocal)
859  {
860  // The face mapped to is on another process.
861  meshFaceId = pIt->first;
862  ASSERTL0(graph[2].count(meshFaceId) == 0,
863  "This periodic boundary edge has been specified before");
864  boost::add_vertex(boostGraphObj);
865  tempGraph[2][meshFaceId] = tempGraphVertId++;
866  nFaceIntCoeffs = FaceSize[meshFaceId];
867  m_numNonDirFaceModes+=nFaceIntCoeffs;
869  }
870  else if (pIt->first < pIt->second[0].id)
871  {
872  ASSERTL0(graph[2].count(pIt->first) == 0,
873  "This periodic boundary face has been specified before");
874  ASSERTL0(graph[2].count(pIt->second[0].id) == 0,
875  "This periodic boundary face has been specified before");
876 
877  boost::add_vertex(boostGraphObj);
878  tempGraph[2][pIt->first] = tempGraphVertId;
879  tempGraph[2][pIt->second[0].id] = tempGraphVertId++;
880  nFaceIntCoeffs = FaceSize[pIt->first];
881  m_numNonDirFaceModes+=nFaceIntCoeffs;
883  }
884  }
885 
886  // setup face numbering
887  for(i = 0; i < locExpVector.size(); ++i)
888  {
889  exp = locExpVector[i];
890  nFaces = exp->GetNfaces();
891  faceCnt = 0;
892  for(j = 0; j < nFaces; ++j)
893  {
894  nFaceIntCoeffs = exp->GetFaceIntNcoeffs(j);
895  meshFaceId = exp->GetGeom()->GetFid(j);
896  if(graph[2].count(meshFaceId) == 0)
897  {
898  if(tempGraph[2].count(meshFaceId) == 0)
899  {
900  boost::add_vertex(boostGraphObj);
901  tempGraph[2][meshFaceId] = tempGraphVertId++;
902  m_numNonDirFaceModes+=nFaceIntCoeffs;
903 
905  }
906  localFaces[localFaceOffset+faceCnt++] = tempGraph[2][meshFaceId];
907  vwgts_map[ tempGraph[2][meshFaceId] ] = nFaceIntCoeffs;
908  }
909  }
910  m_numLocalBndCoeffs += exp->NumBndryCoeffs();
911 
912  localFaceOffset+=nFaces;
913  }
914 
915  localVertOffset=0;
916  localEdgeOffset=0;
917  localFaceOffset=0;
918  for(i = 0; i < locExpVector.size(); ++i)
919  {
920  exp = locExpVector[i];
921  nVerts = exp->GetNverts();
922  nEdges = exp->GetNedges();
923  nFaces = exp->GetNfaces();
924 
925  // Now loop over all local faces, edges and vertices of this
926  // element and define that all other faces, edges and verices of
927  // this element are adjacent to them.
928 
929  // Vertices
930  for(j = 0; j < nVerts; j++)
931  {
932  if(localVerts[j+localVertOffset]==-1)
933  {
934  break;
935  }
936  // associate to other vertices
937  for(k = 0; k < nVerts; k++)
938  {
939  if(localVerts[k+localVertOffset]==-1)
940  {
941  break;
942  }
943  if(k!=j)
944  {
945  boost::add_edge( (size_t) localVerts[j+localVertOffset],
946  (size_t) localVerts[k+localVertOffset],boostGraphObj);
947  }
948  }
949  // associate to other edges
950  for(k = 0; k < nEdges; k++)
951  {
952  if(localEdges[k+localEdgeOffset]==-1)
953  {
954  break;
955  }
956  boost::add_edge( (size_t) localVerts[j+localVertOffset],
957  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
958  }
959  // associate to other faces
960  for(k = 0; k < nFaces; k++)
961  {
962  if(localFaces[k+localFaceOffset]==-1)
963  {
964  break;
965  }
966  boost::add_edge( (size_t) localVerts[j+localVertOffset],
967  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
968  }
969  }
970 
971  // Edges
972  for(j = 0; j < nEdges; j++)
973  {
974  if(localEdges[j+localEdgeOffset]==-1)
975  {
976  break;
977  }
978  // Associate to other edges
979  for(k = 0; k < nEdges; k++)
980  {
981  if(localEdges[k+localEdgeOffset]==-1)
982  {
983  break;
984  }
985  if(k!=j)
986  {
987  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
988  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
989  }
990  }
991  // Associate to vertices
992  for(k = 0; k < nVerts; k++)
993  {
994  if(localVerts[k+localVertOffset]==-1)
995  {
996  break;
997  }
998  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
999  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1000  }
1001  // Associate to faces
1002  for(k = 0; k < nFaces; k++)
1003  {
1004  if(localFaces[k+localFaceOffset]==-1)
1005  {
1006  break;
1007  }
1008  boost::add_edge( (size_t) localEdges[j+localEdgeOffset],
1009  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1010  }
1011  }
1012 
1013  // Faces
1014  for(j = 0; j < nFaces; j++)
1015  {
1016  if(localFaces[j+localFaceOffset]==-1)
1017  {
1018  break;
1019  }
1020  // Associate to other faces
1021  for(k = 0; k < nFaces; k++)
1022  {
1023  if(localFaces[k+localFaceOffset]==-1)
1024  {
1025  break;
1026  }
1027  if(k!=j)
1028  {
1029  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1030  (size_t) localFaces[k+localFaceOffset],boostGraphObj);
1031  }
1032  }
1033  // Associate to vertices
1034  for(k = 0; k < nVerts; k++)
1035  {
1036  if(localVerts[k+localVertOffset]==-1)
1037  {
1038  break;
1039  }
1040  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1041  (size_t) localVerts[k+localVertOffset],boostGraphObj);
1042  }
1043  // Associate to edges
1044  for(k = 0; k < nEdges; k++)
1045  {
1046  if(localEdges[k+localEdgeOffset]==-1)
1047  {
1048  break;
1049  }
1050  boost::add_edge( (size_t) localFaces[j+localFaceOffset],
1051  (size_t) localEdges[k+localEdgeOffset],boostGraphObj);
1052  }
1053  }
1054 
1055  localVertOffset+=nVerts;
1056  localEdgeOffset+=nEdges;
1057  localFaceOffset+=nFaces;
1058  }
1059 
1060  // Container to store vertices of the graph which correspond to
1061  // degrees of freedom along the boundary and periodic BCs.
1062  set<int> partVerts;
1063 
1066  {
1067  vector<long> procVerts, procEdges, procFaces;
1068  set <int> foundVerts, foundEdges, foundFaces;
1069 
1070  // Loop over element and construct the procVerts and procEdges
1071  // vectors, which store the geometry IDs of mesh vertices and
1072  // edges respectively which are local to this process.
1073  for(i = cnt = 0; i < locExpVector.size(); ++i)
1074  {
1075  int elmtid = locExp.GetOffset_Elmt_Id(i);
1076  exp = locExpVector[elmtid];
1077  for (j = 0; j < exp->GetNverts(); ++j)
1078  {
1079  int vid = exp->GetGeom()->GetVid(j)+1;
1080  if (foundVerts.count(vid) == 0)
1081  {
1082  procVerts.push_back(vid);
1083  foundVerts.insert(vid);
1084  }
1085  }
1086 
1087  for (j = 0; j < exp->GetNedges(); ++j)
1088  {
1089  int eid = exp->GetGeom()->GetEid(j)+1;
1090 
1091  if (foundEdges.count(eid) == 0)
1092  {
1093  procEdges.push_back(eid);
1094  foundEdges.insert(eid);
1095  }
1096  }
1097 
1098  for (j = 0; j < exp->GetNfaces(); ++j)
1099  {
1100  int fid = exp->GetGeom()->GetFid(j)+1;
1101 
1102  if (foundFaces.count(fid) == 0)
1103  {
1104  procFaces.push_back(fid);
1105  foundFaces.insert(fid);
1106  }
1107  }
1108  }
1109 
1110  int unique_verts = foundVerts.size();
1111  int unique_edges = foundEdges.size();
1112  int unique_faces = foundFaces.size();
1113 
1114  // Now construct temporary GS objects. These will be used to
1115  // populate the arrays tmp3 and tmp4 with the multiplicity of
1116  // the vertices and edges respectively to identify those
1117  // vertices and edges which are located on partition boundary.
1118  Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
1119  Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm);
1120  Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
1121  Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
1122  Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
1123  Gs::Gather(tmp4, Gs::gs_add, tmp1);
1124  Gs::Finalise(tmp1);
1125 
1126  if (unique_edges > 0)
1127  {
1128  Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
1129  Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
1130  Gs::Gather(tmp5, Gs::gs_add, tmp2);
1131  Gs::Finalise(tmp2);
1132  }
1133 
1134  if (unique_faces > 0)
1135  {
1136  Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
1137  Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
1138  Gs::Gather(tmp6, Gs::gs_add, tmp3);
1139  Gs::Finalise(tmp3);
1140  }
1141 
1142  // Finally, fill the partVerts set with all non-Dirichlet
1143  // vertices which lie on a partition boundary.
1144  for (i = 0; i < unique_verts; ++i)
1145  {
1146  if (tmp4[i] > 1.0)
1147  {
1148  if (graph[0].count(procVerts[i]-1) == 0)
1149  {
1150  partVerts.insert(tempGraph[0][procVerts[i]-1]);
1151  }
1152  }
1153  }
1154 
1155  for (i = 0; i < unique_edges; ++i)
1156  {
1157  if (tmp5[i] > 1.0)
1158  {
1159  if (graph[1].count(procEdges[i]-1) == 0)
1160  {
1161  partVerts.insert(tempGraph[1][procEdges[i]-1]);
1162  }
1163  }
1164  }
1165 
1166  for (i = 0; i < unique_faces; ++i)
1167  {
1168  if (tmp6[i] > 1.0)
1169  {
1170  if (graph[2].count(procFaces[i]-1) == 0)
1171  {
1172  partVerts.insert(tempGraph[2][procFaces[i]-1]);
1173  }
1174  }
1175  }
1176 
1177  // Now fill with all vertices on periodic BCs
1178  for (pIt = periodicVerts.begin(); pIt != periodicVerts.end(); ++pIt)
1179  {
1180  if (graph[0].count(pIt->first) == 0)
1181  {
1182  partVerts.insert(tempGraph[0][pIt->first]);
1183  }
1184  }
1185  for (pIt = periodicEdges.begin(); pIt != periodicEdges.end(); ++pIt)
1186  {
1187  if (graph[1].count(pIt->first) == 0)
1188  {
1189  partVerts.insert(tempGraph[1][pIt->first]);
1190  }
1191  }
1192  for (pIt = periodicFaces.begin(); pIt != periodicFaces.end(); ++pIt)
1193  {
1194  if (graph[2].count(pIt->first) == 0)
1195  {
1196  partVerts.insert(tempGraph[2][pIt->first]);
1197  }
1198  }
1199  }
1200 
1201  int nGraphVerts = tempGraphVertId;
1202  Array<OneD, int> perm(nGraphVerts);
1203  Array<OneD, int> iperm(nGraphVerts);
1204  Array<OneD, int> vwgts(nGraphVerts);
1205  ASSERTL1(vwgts_map.size()==nGraphVerts,"Non matching dimensions");
1206  for(i = 0; i < nGraphVerts; ++i)
1207  {
1208  vwgts[i] = vwgts_map[i];
1209  }
1210 
1211  if(nGraphVerts)
1212  {
1213  switch(m_solnType)
1214  {
1215  case eDirectFullMatrix:
1216  case eIterativeFull:
1217  case eIterativeStaticCond:
1218  case ePETScStaticCond:
1219  case ePETScFullMatrix:
1220  case eXxtFullMatrix:
1221  case eXxtStaticCond:
1222  {
1223  NoReordering(boostGraphObj,perm,iperm);
1224  break;
1225  }
1226 
1227  case eDirectStaticCond:
1228  {
1229  CuthillMckeeReordering(boostGraphObj,perm,iperm);
1230  break;
1231  }
1232 
1237  {
1239  boostGraphObj, perm, iperm, bottomUpGraph,
1240  partVerts, mdswitch);
1241  break;
1242  }
1243  default:
1244  {
1245  ASSERTL0(false,
1246  "Unrecognised solution type: " + std::string(
1248  }
1249  }
1250  }
1251 
1252  // For parallel multi-level static condensation determine the lowest
1253  // static condensation level amongst processors.
1256  {
1257  m_lowestStaticCondLevel = bottomUpGraph->GetNlevels()-1;
1258  vComm->AllReduce(m_lowestStaticCondLevel,
1260  }
1261  else
1262  {
1264  }
1265 
1266  /**
1267  * STEP 4: Fill the #graph[0] and
1268  * #graph[1] with the optimal ordering from boost.
1269  */
1271  for(mapIt = tempGraph[0].begin(); mapIt != tempGraph[0].end(); mapIt++)
1272  {
1273  graph[0][mapIt->first] = iperm[mapIt->second] + graphVertId;
1274  }
1275  for(mapIt = tempGraph[1].begin(); mapIt != tempGraph[1].end(); mapIt++)
1276  {
1277  graph[1][mapIt->first] = iperm[mapIt->second] + graphVertId;
1278  }
1279  for(mapIt = tempGraph[2].begin(); mapIt != tempGraph[2].end(); mapIt++)
1280  {
1281  graph[2][mapIt->first] = iperm[mapIt->second] + graphVertId;
1282  }
1283 
1284  return nGraphVerts;
1285  }
1286 
1287  /**
1288  *
1289  */
1291  const LibUtilities::SessionReaderSharedPtr &pSession,
1292  const int numLocalCoeffs,
1293  const ExpList &locExp,
1294  const BndCondExp &bndCondExp,
1295  const BndCond &bndConditions,
1296  const bool checkIfSystemSingular,
1297  const std::string variable,
1298  const PeriodicMap &periodicVerts,
1299  const PeriodicMap &periodicEdges,
1300  const PeriodicMap &periodicFaces)
1301  : AssemblyMap(pSession, variable)
1302  {
1303  int i, j, k, l;
1304  int cnt = 0;
1305  int intDofCnt;
1306  int meshVertId, meshEdgeId, meshFaceId;
1307  int globalId;
1308  int nEdgeInteriorCoeffs;
1309  int nFaceInteriorCoeffs;
1310  int firstNonDirGraphVertId;
1311  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
1314  StdRegions::Orientation edgeOrient;
1315  StdRegions::Orientation faceOrient;
1316  Array<OneD, unsigned int> edgeInteriorMap;
1317  Array<OneD, int> edgeInteriorSign;
1318  Array<OneD, unsigned int> faceInteriorMap;
1319  Array<OneD, int> faceInteriorSign;
1320  PeriodicMap::const_iterator pIt;
1321 
1322  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1323 
1324  m_signChange = false;
1325 
1326  // Stores vertex, edge and face reordered vertices.
1327  DofGraph graph(3);
1328  DofGraph dofs(3);
1329 
1330  set<int> extraDirVerts, extraDirEdges;
1332 
1333  // Construct list of number of degrees of freedom for each vertex,
1334  // edge and face.
1335  for (i = 0; i < locExpVector.size(); ++i)
1336  {
1337  exp = locExpVector[i];
1338 
1339  for(j = 0; j < locExpVector[i]->GetNverts(); ++j)
1340  {
1341  dofs[0][exp->GetGeom()->GetVid(j)] = 1;
1342  }
1343 
1344  for(j = 0; j < locExpVector[i]->GetNedges(); ++j)
1345  {
1346  if (dofs[1].count(exp->GetGeom()->GetEid(j)) > 0)
1347  {
1348  if (dofs[1][exp->GetGeom()->GetEid(j)] !=
1349  locExpVector[i]->GetEdgeNcoeffs(j)-2)
1350  {
1351  ASSERTL0( (exp->GetEdgeBasisType(j) == LibUtilities::eModified_A) ||
1352  (exp->GetEdgeBasisType(j) == LibUtilities::eModified_B),
1353  "CG with variable order only available with modal expansion");
1354  }
1355  dofs[1][exp->GetGeom()->GetEid(j)] =
1356  min(dofs[1][exp->GetGeom()->GetEid(j)],
1357  locExpVector[i]->GetEdgeNcoeffs(j)-2);
1358  }
1359  else
1360  {
1361  dofs[1][exp->GetGeom()->GetEid(j)] =
1362  exp->GetEdgeNcoeffs(j) - 2;
1363  }
1364  }
1365 
1366  for(j = 0; j < locExpVector[i]->GetNfaces(); ++j)
1367  {
1368  if (dofs[2].count(exp->GetGeom()->GetFid(j)) > 0)
1369  {
1370  if (dofs[2][exp->GetGeom()->GetFid(j)] !=
1371  exp->GetFaceIntNcoeffs(j))
1372  {
1373  ASSERTL0( false,
1374  "CG with variable order not available in 3D");
1375  }
1376  dofs[2][exp->GetGeom()->GetFid(j)] =
1377  min(dofs[2][exp->GetGeom()->GetFid(j)],
1378  exp->GetFaceIntNcoeffs(j));
1379  }
1380  else
1381  {
1382  dofs[2][exp->GetGeom()->GetFid(j)] =
1383  exp->GetFaceIntNcoeffs(j);
1384  }
1385  }
1386  }
1387  // Now use information from all partitions to determine
1388  // the correct size
1390  // edges
1391  Array<OneD, long> edgeId (dofs[1].size());
1392  Array<OneD, NekDouble> edgeDof (dofs[1].size());
1393  for(dofIt = dofs[1].begin(), i=0; dofIt != dofs[1].end(); dofIt++, i++)
1394  {
1395  edgeId[i] = dofIt->first;
1396  edgeDof[i] = (NekDouble) dofIt->second;
1397  }
1398  Gs::gs_data *tmp = Gs::Init(edgeId, vComm);
1399  Gs::Gather(edgeDof, Gs::gs_min, tmp);
1400  Gs::Finalise(tmp);
1401  for (i=0; i < dofs[1].size(); i++)
1402  {
1403  dofs[1][edgeId[i]] = (int) (edgeDof[i]+0.5);
1404  }
1405  // faces
1406  Array<OneD, long> faceId (dofs[2].size());
1407  Array<OneD, NekDouble> faceDof (dofs[2].size());
1408  for(dofIt = dofs[2].begin(), i=0; dofIt != dofs[2].end(); dofIt++, i++)
1409  {
1410  faceId[i] = dofIt->first;
1411  faceDof[i] = (NekDouble) dofIt->second;
1412  }
1413  Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
1414  Gs::Gather(faceDof, Gs::gs_min, tmp2);
1415  Gs::Finalise(tmp2);
1416  for (i=0; i < dofs[2].size(); i++)
1417  {
1418  dofs[2][faceId[i]] = (int) (faceDof[i]+0.5);
1419  }
1420 
1421  Array<OneD, const BndCond> bndCondVec(1, bndConditions);
1422 
1423  // Note that nExtraDirichlet is not used in the logic below; it just
1424  // needs to be set so that the coupled solver in
1425  // IncNavierStokesSolver can work.
1426  int nExtraDirichlet;
1427  int nGraphVerts =
1428  CreateGraph(locExp, bndCondExp, bndCondVec,
1429  checkIfSystemSingular, periodicVerts, periodicEdges,
1430  periodicFaces, graph, bottomUpGraph, extraDirVerts,
1431  extraDirEdges, firstNonDirGraphVertId,
1432  nExtraDirichlet);
1433 
1434  /*
1435  * Set up an array which contains the offset information of the
1436  * different graph vertices.
1437  *
1438  * This basically means to identify to how many global degrees of
1439  * freedom the individual graph vertices correspond. Obviously,
1440  * the graph vertices corresponding to the mesh-vertices account
1441  * for a single global DOF. However, the graph vertices
1442  * corresponding to the element edges correspond to N-2 global DOF
1443  * where N is equal to the number of boundary modes on this edge.
1444  */
1445  Array<OneD, int> graphVertOffset(
1446  graph[0].size() + graph[1].size() + graph[2].size() + 1);
1447 
1448  graphVertOffset[0] = 0;
1449 
1450  for(i = 0; i < locExpVector.size(); ++i)
1451  {
1452  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1453 
1454  for(j = 0; j < exp->GetNverts(); ++j)
1455  {
1456  meshVertId = exp->GetGeom()->GetVid(j);
1457  graphVertOffset[graph[0][meshVertId]+1] = 1;
1458  }
1459 
1460  for(j = 0; j < exp->GetNedges(); ++j)
1461  {
1462  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j) - 2;
1463  meshEdgeId = exp->GetGeom()->GetEid(j);
1464  graphVertOffset[graph[1][meshEdgeId]+1]
1465  = dofs[1][meshEdgeId];
1466 
1467  bType = exp->GetEdgeBasisType(j);
1468 
1469  // need a sign vector for modal expansions if nEdgeCoeffs
1470  // >=3 (not 4 because of variable order case)
1471  if(nEdgeInteriorCoeffs &&
1472  (bType == LibUtilities::eModified_A ||
1473  bType == LibUtilities::eModified_B))
1474  {
1475  m_signChange = true;
1476  }
1477  }
1478 
1479  for(j = 0; j < exp->GetNfaces(); ++j)
1480  {
1481  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1482  meshFaceId = exp->GetGeom()->GetFid(j);
1483  graphVertOffset[graph[2][meshFaceId]+1] = dofs[2][meshFaceId];
1484  }
1485  }
1486 
1487  for(i = 1; i < graphVertOffset.num_elements(); i++)
1488  {
1489  graphVertOffset[i] += graphVertOffset[i-1];
1490  }
1491 
1492  // Allocate the proper amount of space for the class-data
1493  m_numLocalCoeffs = numLocalCoeffs;
1494  m_numGlobalDirBndCoeffs = graphVertOffset[firstNonDirGraphVertId];
1498  // If required, set up the sign-vector
1499  if(m_signChange)
1500  {
1504  }
1505 
1506  m_staticCondLevel = 0;
1507  m_numPatches = locExpVector.size();
1510  for(i = 0; i < m_numPatches; ++i)
1511  {
1512  m_numLocalBndCoeffsPerPatch[i] = (unsigned int)
1513  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1514  m_numLocalIntCoeffsPerPatch[i] = (unsigned int)
1515  locExpVector[locExp.GetOffset_Elmt_Id(i)]->GetNcoeffs() -
1516  locExpVector[locExp.GetOffset_Elmt_Id(i)]->NumBndryCoeffs();
1517  }
1518 
1519  /**
1520  * STEP 6: Now, all ingredients are ready to set up the actual
1521  * local to global mapping.
1522  *
1523  * The remainder of the map consists of the element-interior
1524  * degrees of freedom. This leads to the block-diagonal submatrix
1525  * as each element-interior mode is globally orthogonal to modes
1526  * in all other elements.
1527  */
1528  cnt = 0;
1529 
1530  // Loop over all the elements in the domain
1531  for(i = 0; i < locExpVector.size(); ++i)
1532  {
1533  exp = locExpVector[i];
1534  cnt = locExp.GetCoeff_Offset(i);
1535  for(j = 0; j < exp->GetNverts(); ++j)
1536  {
1537  meshVertId = exp->GetGeom()->GetVid(j);
1538 
1539  // Set the global DOF for vertex j of element i
1540  m_localToGlobalMap[cnt+exp->GetVertexMap(j)] =
1541  graphVertOffset[graph[0][meshVertId]];
1542  }
1543 
1544  for(j = 0; j < exp->GetNedges(); ++j)
1545  {
1546  nEdgeInteriorCoeffs = exp->GetEdgeNcoeffs(j)-2;
1547  edgeOrient = exp->GetGeom()->GetEorient(j);
1548  meshEdgeId = exp->GetGeom()->GetEid(j);
1549 
1550  pIt = periodicEdges.find(meshEdgeId);
1551 
1552  // See if this edge is periodic. If it is, then we map all
1553  // edges to the one with lowest ID, and align all
1554  // coefficients to this edge orientation.
1555  if (pIt != periodicEdges.end())
1556  {
1557  pair<int, StdRegions::Orientation> idOrient =
1559  meshEdgeId, edgeOrient, pIt->second);
1560  edgeOrient = idOrient.second;
1561  }
1562 
1563  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1564 
1565  // Set the global DOF's for the interior modes of edge j
1566  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1567  {
1568  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1569  graphVertOffset[graph[1][meshEdgeId]]+k;
1570  }
1571  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1572  {
1573  m_localToGlobalMap[cnt+edgeInteriorMap[k]] =
1574  graphVertOffset[graph[1][meshEdgeId]];
1575  }
1576 
1577  // Fill the sign vector if required
1578  if(m_signChange)
1579  {
1580  for(k = 0; k < dofs[1][exp->GetGeom()->GetEid(j)]; ++k)
1581  {
1582  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = (NekDouble) edgeInteriorSign[k];
1583  }
1584  for(k = dofs[1][exp->GetGeom()->GetEid(j)]; k < nEdgeInteriorCoeffs; ++k)
1585  {
1586  m_localToGlobalSign[cnt+edgeInteriorMap[k]] = 0.0;
1587  }
1588  }
1589  }
1590 
1591  for(j = 0; j < exp->GetNfaces(); ++j)
1592  {
1593  nFaceInteriorCoeffs = exp->GetFaceIntNcoeffs(j);
1594  faceOrient = exp->GetGeom()->GetForient(j);
1595  meshFaceId = exp->GetGeom()->GetFid(j);
1596 
1597  pIt = periodicFaces.find(meshFaceId);
1598 
1599  if (pIt != periodicFaces.end() &&
1600  meshFaceId == min(meshFaceId, pIt->second[0].id))
1601  {
1602  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
1603  }
1604 
1605  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
1606 
1607  // Set the global DOF's for the interior modes of face j
1608  for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
1609  {
1610  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1611  graphVertOffset[graph[2][meshFaceId]]+k;
1612  }
1613  for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
1614  {
1615  m_localToGlobalMap[cnt+faceInteriorMap[k]] =
1616  graphVertOffset[graph[2][meshFaceId]];
1617  }
1618 
1619  if(m_signChange)
1620  {
1621  for(k = 0; k < dofs[2][exp->GetGeom()->GetFid(j)]; ++k)
1622  {
1623  m_localToGlobalSign[cnt+faceInteriorMap[k]] = (NekDouble) faceInteriorSign[k];
1624  }
1625  for(k = dofs[2][exp->GetGeom()->GetFid(j)]; k < nFaceInteriorCoeffs; ++k)
1626  {
1627  m_localToGlobalSign[cnt+faceInteriorMap[k]] = 0;
1628  }
1629  }
1630 
1631  }
1632  }
1633 
1634  // Set up the mapping for the boundary conditions
1635  cnt = 0;
1636  int offset = 0;
1637  for(i = 0; i < bndCondExp.num_elements(); i++)
1638  {
1639  set<int> foundExtraVerts, foundExtraEdges;
1640  for(j = 0; j < bndCondExp[i]->GetNumElmts(); j++)
1641  {
1642  bndExp = bndCondExp[i]->GetExp(j);
1643  cnt = offset + bndCondExp[i]->GetCoeff_Offset(j);
1644  for(k = 0; k < bndExp->GetNverts(); k++)
1645  {
1646  meshVertId = bndExp->GetGeom()->GetVid(k);
1647  m_bndCondCoeffsToGlobalCoeffsMap[cnt+bndExp->GetVertexMap(k)] = graphVertOffset[graph[0][meshVertId]];
1648 
1649  if (bndConditions[i]->GetBoundaryConditionType() !=
1651  {
1652  continue;
1653  }
1654 
1655  set<int>::iterator iter = extraDirVerts.find(meshVertId);
1656  if (iter != extraDirVerts.end() &&
1657  foundExtraVerts.count(meshVertId) == 0)
1658  {
1659  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1660  bndExp->GetVertexMap(k);
1661  int gid = graphVertOffset[
1662  graph[0][meshVertId]];
1663  ExtraDirDof t(loc, gid, 1.0);
1664  m_extraDirDofs[i].push_back(t);
1665  foundExtraVerts.insert(meshVertId);
1666  }
1667  }
1668 
1669  for(k = 0; k < bndExp->GetNedges(); k++)
1670  {
1671  nEdgeInteriorCoeffs = bndExp->GetEdgeNcoeffs(k)-2;
1672  edgeOrient = bndExp->GetGeom()->GetEorient(k);
1673  meshEdgeId = bndExp->GetGeom()->GetEid(k);
1674 
1675  pIt = periodicEdges.find(meshEdgeId);
1676 
1677  // See if this edge is periodic. If it is, then we map
1678  // all edges to the one with lowest ID, and align all
1679  // coefficients to this edge orientation.
1680  if (pIt != periodicEdges.end())
1681  {
1682  pair<int, StdRegions::Orientation> idOrient =
1684  meshEdgeId, edgeOrient, pIt->second);
1685  edgeOrient = idOrient.second;
1686  }
1687 
1688  bndExp->GetEdgeInteriorMap(
1689  k,edgeOrient,edgeInteriorMap,edgeInteriorSign);
1690 
1691  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1692  {
1693  m_bndCondCoeffsToGlobalCoeffsMap[cnt+edgeInteriorMap[l]] =
1694  graphVertOffset[graph[1][meshEdgeId]]+l;
1695  }
1696 
1697  // Fill the sign vector if required
1698  if(m_signChange)
1699  {
1700  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1701  {
1702  m_bndCondCoeffsToGlobalCoeffsSign[cnt+edgeInteriorMap[l]] = (NekDouble) edgeInteriorSign[l];
1703  }
1704  }
1705 
1706  if (bndConditions[i]->GetBoundaryConditionType() !=
1708  {
1709  continue;
1710  }
1711 
1712  set<int>::iterator iter = extraDirEdges.find(meshEdgeId);
1713  if (iter != extraDirEdges.end() &&
1714  foundExtraEdges.count(meshEdgeId) == 0 &&
1715  nEdgeInteriorCoeffs > 0)
1716  {
1717  for(l = 0; l < nEdgeInteriorCoeffs; ++l)
1718  {
1719  int loc = bndCondExp[i]->GetCoeff_Offset(j) +
1720  edgeInteriorMap[l];
1721  int gid = graphVertOffset[
1722  graph[1][meshEdgeId]]+l;
1723  ExtraDirDof t(loc, gid, edgeInteriorSign[l]);
1724  m_extraDirDofs[i].push_back(t);
1725  }
1726  foundExtraEdges.insert(meshEdgeId);
1727  }
1728  }
1729 
1730  meshFaceId = bndExp->GetGeom()->GetGlobalID();
1731  intDofCnt = 0;
1732  for(k = 0; k < bndExp->GetNcoeffs(); k++)
1733  {
1734  if(m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] == -1)
1735  {
1736  m_bndCondCoeffsToGlobalCoeffsMap[cnt+k] =
1737  graphVertOffset[graph[bndExp->GetNumBases()][meshFaceId]]+intDofCnt;
1738  intDofCnt++;
1739  }
1740  }
1741  }
1742  offset += bndCondExp[i]->GetNcoeffs();
1743  }
1744 
1745  globalId = Vmath::Vmax(m_numLocalCoeffs,&m_localToGlobalMap[0],1)+1;
1746  m_numGlobalBndCoeffs = globalId;
1747 
1748 
1749  /*
1750  * The boundary condition mapping is generated from the same vertex
1751  * renumbering.
1752  */
1753  cnt=0;
1754  for(i = 0; i < m_numLocalCoeffs; ++i)
1755  {
1756  if(m_localToGlobalMap[i] == -1)
1757  {
1758  m_localToGlobalMap[i] = globalId++;
1759  }
1760  else
1761  {
1762  if(m_signChange)
1763  {
1765  }
1766  m_localToGlobalBndMap[cnt++]=m_localToGlobalMap[i];
1767  }
1768  }
1769 
1770  m_numGlobalCoeffs = globalId;
1771 
1772  SetUpUniversalC0ContMap(locExp, periodicVerts, periodicEdges, periodicFaces);
1773 
1774  // Since we can now have multiple entries to m_extraDirDofs due to
1775  // periodic boundary conditions we make a call to work out the
1776  // multiplicity of all entries and invert (Need to be after
1777  // SetupUniversalC0ContMap)
1779 
1780  // Fill in Dirichlet coefficients that are to be sent to other
1781  // processors with a value of 1
1782  map<int, vector<ExtraDirDof> >::iterator Tit;
1783 
1784  // Generate valence for extraDirDofs
1785  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1786  {
1787  for (i = 0; i < Tit->second.size(); ++i)
1788  {
1789  valence[Tit->second[i].get<1>()] = 1.0;
1790  }
1791  }
1792 
1793  UniversalAssembleBnd(valence);
1794 
1795  // Set third argument of tuple to inverse of valence.
1796  for (Tit = m_extraDirDofs.begin(); Tit != m_extraDirDofs.end(); ++Tit)
1797  {
1798  for (i = 0; i < Tit->second.size(); ++i)
1799  {
1800  boost::get<2>(Tit->second.at(i)) /= valence[Tit->second.at(i).get<1>()];
1801  }
1802  }
1803 
1804  // Set up the local to global map for the next level when using
1805  // multi-level static condensation
1809  m_solnType == ePETScMultiLevelStaticCond) && nGraphVerts)
1810  {
1811  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1812  {
1813  Array<OneD, int> vwgts_perm(
1814  dofs[0].size() + dofs[1].size() + dofs[2].size()
1815  - firstNonDirGraphVertId);
1816 
1817  for (i = 0; i < locExpVector.size(); ++i)
1818  {
1819  exp = locExpVector[locExp.GetOffset_Elmt_Id(i)];
1820 
1821  for (j = 0; j < exp->GetNverts(); ++j)
1822  {
1823  meshVertId = exp->GetGeom()->GetVid(j);
1824 
1825  if (graph[0][meshVertId] >= firstNonDirGraphVertId)
1826  {
1827  vwgts_perm[graph[0][meshVertId] -
1828  firstNonDirGraphVertId] =
1829  dofs[0][meshVertId];
1830  }
1831  }
1832 
1833  for (j = 0; j < exp->GetNedges(); ++j)
1834  {
1835  meshEdgeId = exp->GetGeom()->GetEid(j);
1836 
1837  if (graph[1][meshEdgeId] >= firstNonDirGraphVertId)
1838  {
1839  vwgts_perm[graph[1][meshEdgeId] -
1840  firstNonDirGraphVertId] =
1841  dofs[1][meshEdgeId];
1842  }
1843  }
1844 
1845  for (j = 0; j < exp->GetNfaces(); ++j)
1846  {
1847  meshFaceId = exp->GetGeom()->GetFid(j);
1848 
1849  if (graph[2][meshFaceId] >= firstNonDirGraphVertId)
1850  {
1851  vwgts_perm[graph[2][meshFaceId] -
1852  firstNonDirGraphVertId] =
1853  dofs[2][meshFaceId];
1854  }
1855  }
1856  }
1857 
1858  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1860  AllocateSharedPtr(this, bottomUpGraph);
1861  }
1862  }
1863 
1864  m_hash = boost::hash_range(m_localToGlobalMap.begin(),
1865  m_localToGlobalMap.end());
1866 
1867  // Add up hash values if parallel
1868  int hash = m_hash;
1869  vComm->AllReduce(hash, LibUtilities::ReduceSum);
1870  m_hash = hash;
1871 
1874  }
1875 
1876  /**
1877  *
1878  */
1880  {
1883  }
1884 
1885  /**
1886  * @brief Determine orientation of an edge to its periodic equivalents,
1887  * as well as the ID of the representative edge.
1888  *
1889  * Since an edge may be periodic with more than one other edge (e.g. a
1890  * periodic cube has sets of four periodic edges in each coordinate
1891  * direction), we have to define a 'representative' edge. In this
1892  * assembly map we define it to be the one with the minimum ID. This
1893  * routine is set up to calculate the orientation of a given edge with
1894  * ID @p meshEdgeId with respect to the edge ID.
1895  *
1896  * @param meshEdgeId ID of a periodic edge.
1897  * @param edgeOrient Edge orientation of meshEdgeId with respect to
1898  * its parent element.
1899  * @param periodicEdges The map of all periodic edges.
1900  *
1901  * @return Pair containing the ID of the periodic edge and the
1902  * orientation of @p meshEdgeID with respect to this edge.
1903  */
1904  pair<int, StdRegions::Orientation> DeterminePeriodicEdgeOrientId(
1905  int meshEdgeId,
1906  StdRegions::Orientation edgeOrient,
1907  const vector<PeriodicEntity> &periodicEdges)
1908  {
1909  int minId = periodicEdges[0].id;
1910  int minIdK = 0;
1911  int k;
1912 
1913  for (k = 1; k < periodicEdges.size(); ++k)
1914  {
1915  if (periodicEdges[k].id < minId)
1916  {
1917  minId = min(minId, periodicEdges[k].id);
1918  minIdK = k;
1919  }
1920  }
1921 
1922  minId = min(minId, meshEdgeId);
1923 
1924  if (meshEdgeId != minId)
1925  {
1926  if (periodicEdges[minIdK].orient == StdRegions::eBackwards)
1927  {
1928  // Swap edge orientation
1929  edgeOrient = (edgeOrient == StdRegions::eForwards) ?
1931  }
1932  }
1933 
1934  return make_pair(minId, edgeOrient);
1935  }
1936 
1937  /**
1938  * @brief Determine relative orientation between two faces.
1939  *
1940  * Given the orientation of a local element to its local face, defined
1941  * as @p faceOrient, and @p perFaceOrient which states the alignment of
1942  * one periodic face to the other global face, this routine determines
1943  * the orientation that takes this local element face to the
1944  * global/unique face.
1945  *
1946  * @param faceOrient Orientation of the face with respect to its
1947  * parent element.
1948  * @param perFaceOrient Orientation of the representative/global face.
1949  *
1950  * @return Orientation between the two faces.
1951  */
1953  StdRegions::Orientation faceOrient,
1954  StdRegions::Orientation perFaceOrient)
1955  {
1956  StdRegions::Orientation returnval = faceOrient;
1957 
1958  if(perFaceOrient != StdRegions::eDir1FwdDir1_Dir2FwdDir2)
1959  {
1960  int tmp1 = (int)faceOrient - 5;
1961  int tmp2 = (int)perFaceOrient - 5;
1962 
1963  int flipDir1Map [8] = {2,3,0,1,6,7,4,5};
1964  int flipDir2Map [8] = {1,0,3,2,5,4,7,6};
1965  int transposeMap[8] = {4,5,6,7,0,2,1,3};
1966 
1967  // Transpose orientation
1968  if (tmp2 > 3)
1969  {
1970  tmp1 = transposeMap[tmp1];
1971  }
1972 
1973  // Reverse orientation in direction 1.
1974  if (tmp2 == 2 || tmp2 == 3 || tmp2 == 6 || tmp2 == 7)
1975  {
1976  tmp1 = flipDir1Map[tmp1];
1977  }
1978 
1979  // Reverse orientation in direction 2
1980  if (tmp2 % 2 == 1)
1981  {
1982  tmp1 = flipDir2Map[tmp1];
1983  }
1984 
1985  returnval = (StdRegions::Orientation)(tmp1+5);
1986  }
1987  return returnval;
1988  }
1989 
1990 
1991  /**
1992  * Sets up the global to universal mapping of degrees of freedom across
1993  * processors.
1994  */
1996  const ExpList &locExp,
1997  const PeriodicMap &perVerts,
1998  const PeriodicMap &perEdges,
1999  const PeriodicMap &perFaces)
2000  {
2002  int nVert = 0;
2003  int nEdge = 0;
2004  int nFace = 0;
2005  int maxEdgeDof = 0;
2006  int maxFaceDof = 0;
2007  int maxIntDof = 0;
2008  int dof = 0;
2009  int cnt;
2010  int i,j,k;
2011  int meshVertId;
2012  int meshEdgeId;
2013  int meshFaceId;
2014  int elementId;
2015  int vGlobalId;
2016  int maxBndGlobalId = 0;
2017  StdRegions::Orientation edgeOrient;
2018  StdRegions::Orientation faceOrient;
2019  Array<OneD, unsigned int> edgeInteriorMap;
2020  Array<OneD, int> edgeInteriorSign;
2021  Array<OneD, unsigned int> faceInteriorMap;
2022  Array<OneD, int> faceInteriorSign;
2023  Array<OneD, unsigned int> interiorMap;
2024  PeriodicMap::const_iterator pIt;
2025 
2026  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
2027  LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
2028 
2033 
2034  // Loop over all the elements in the domain to gather mesh data
2035  for(i = 0; i < locExpVector.size(); ++i)
2036  {
2037  exp = locExpVector[i];
2038  nVert += exp->GetNverts();
2039  nEdge += exp->GetNedges();
2040  nFace += exp->GetNfaces();
2041  // Loop over all edges (and vertices) of element i
2042  for(j = 0; j < exp->GetNedges(); ++j)
2043  {
2044  dof = exp->GetEdgeNcoeffs(j)-2;
2045  maxEdgeDof = (dof > maxEdgeDof ? dof : maxEdgeDof);
2046  }
2047  for(j = 0; j < exp->GetNfaces(); ++j)
2048  {
2049  dof = exp->GetFaceIntNcoeffs(j);
2050  maxFaceDof = (dof > maxFaceDof ? dof : maxFaceDof);
2051  }
2052  exp->GetInteriorMap(interiorMap);
2053  dof = interiorMap.num_elements();
2054  maxIntDof = (dof > maxIntDof ? dof : maxIntDof);
2055  }
2056 
2057  // Tell other processes about how many dof we have
2058  vCommRow->AllReduce(nVert, LibUtilities::ReduceSum);
2059  vCommRow->AllReduce(nEdge, LibUtilities::ReduceSum);
2060  vCommRow->AllReduce(nFace, LibUtilities::ReduceSum);
2061  vCommRow->AllReduce(maxEdgeDof, LibUtilities::ReduceMax);
2062  vCommRow->AllReduce(maxFaceDof, LibUtilities::ReduceMax);
2063  vCommRow->AllReduce(maxIntDof, LibUtilities::ReduceMax);
2064 
2065  // Assemble global to universal mapping for this process
2066  for(i = 0; i < locExpVector.size(); ++i)
2067  {
2068  exp = locExpVector[i];
2069  cnt = locExp.GetCoeff_Offset(i);
2070 
2071  // Loop over all vertices of element i
2072  for(j = 0; j < exp->GetNverts(); ++j)
2073  {
2074  meshVertId = exp->GetGeom()->GetVid(j);
2075  vGlobalId = m_localToGlobalMap[cnt+exp->GetVertexMap(j)];
2076 
2077  pIt = perVerts.find(meshVertId);
2078  if (pIt != perVerts.end())
2079  {
2080  for (k = 0; k < pIt->second.size(); ++k)
2081  {
2082  meshVertId = min(meshVertId, pIt->second[k].id);
2083  }
2084  }
2085 
2086  m_globalToUniversalMap[vGlobalId] = meshVertId + 1;
2087  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2088  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2089  }
2090 
2091  // Loop over all edges of element i
2092  for(j = 0; j < exp->GetNedges(); ++j)
2093  {
2094  meshEdgeId = exp->GetGeom()->GetEid(j);
2095  pIt = perEdges.find(meshEdgeId);
2096  edgeOrient = exp->GetGeom()->GetEorient(j);
2097 
2098  if (pIt != perEdges.end())
2099  {
2100  pair<int, StdRegions::Orientation> idOrient =
2102  meshEdgeId, edgeOrient, pIt->second);
2103  meshEdgeId = idOrient.first;
2104  edgeOrient = idOrient.second;
2105  }
2106 
2107  exp->GetEdgeInteriorMap(j,edgeOrient,edgeInteriorMap,edgeInteriorSign);
2108  dof = exp->GetEdgeNcoeffs(j)-2;
2109 
2110  // Set the global DOF's for the interior modes of edge j
2111  // run backwards because of variable P case "ghost" modes
2112  for(k = dof-1; k >= 0; --k)
2113  {
2114  vGlobalId = m_localToGlobalMap[cnt+edgeInteriorMap[k]];
2115  m_globalToUniversalMap[vGlobalId]
2116  = nVert + meshEdgeId * maxEdgeDof + k + 1;
2117  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2118  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2119  }
2120  }
2121 
2122  // Loop over all faces of element i
2123  for(j = 0; j < exp->GetNfaces(); ++j)
2124  {
2125  faceOrient = exp->GetGeom()->GetForient(j);
2126 
2127  meshFaceId = exp->GetGeom()->GetFid(j);
2128 
2129  pIt = perFaces.find(meshFaceId);
2130  if (pIt != perFaces.end())
2131  {
2132  if(meshFaceId == min(meshFaceId, pIt->second[0].id))
2133  {
2134  faceOrient = DeterminePeriodicFaceOrient(faceOrient,pIt->second[0].orient);
2135  }
2136  meshFaceId = min(meshFaceId, pIt->second[0].id);
2137  }
2138 
2139 
2140  exp->GetFaceInteriorMap(j,faceOrient,faceInteriorMap,faceInteriorSign);
2141  dof = exp->GetFaceIntNcoeffs(j);
2142 
2143 
2144  for(k = dof-1; k >= 0; --k)
2145  {
2146  vGlobalId = m_localToGlobalMap[cnt+faceInteriorMap[k]];
2147  m_globalToUniversalMap[vGlobalId]
2148  = nVert + nEdge*maxEdgeDof + meshFaceId * maxFaceDof
2149  + k + 1;
2150  m_globalToUniversalBndMap[vGlobalId]=m_globalToUniversalMap[vGlobalId];
2151 
2152  maxBndGlobalId = (vGlobalId > maxBndGlobalId ? vGlobalId : maxBndGlobalId);
2153  }
2154  }
2155 
2156  // Add interior DOFs to complete universal numbering
2157  exp->GetInteriorMap(interiorMap);
2158  dof = interiorMap.num_elements();
2159  elementId = (exp->GetGeom())->GetGlobalID();
2160  for (k = 0; k < dof; ++k)
2161  {
2162  vGlobalId = m_localToGlobalMap[cnt+interiorMap[k]];
2163  m_globalToUniversalMap[vGlobalId]
2164  = nVert + nEdge*maxEdgeDof + nFace*maxFaceDof + elementId*maxIntDof + k + 1;
2165  }
2166  }
2167 
2168  // Set up the GSLib universal assemble mapping
2169  // Internal DOF do not participate in any data
2170  // exchange, so we keep these set to the special GSLib id=0 so
2171  // they are ignored.
2172  Nektar::Array<OneD, long> tmp(m_numGlobalCoeffs);
2173  Vmath::Zero(m_numGlobalCoeffs, tmp, 1);
2174  Nektar::Array<OneD, long> tmp2(m_numGlobalBndCoeffs, tmp);
2175  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2176  {
2177  tmp[i] = m_globalToUniversalMap[i];
2178  }
2179 
2180  m_gsh = Gs::Init(tmp, vCommRow);
2181  m_bndGsh = Gs::Init(tmp2, vCommRow);
2182  Gs::Unique(tmp, vCommRow);
2183  for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
2184  {
2185  m_globalToUniversalMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
2186  }
2187  for (unsigned int i = 0; i < m_numGlobalBndCoeffs; ++i)
2188  {
2189  m_globalToUniversalBndMapUnique[i] = (tmp2[i] >= 0 ? 1 : 0);
2190  }
2191  }
2192 
2193  /**
2194  * @brief Construct an AssemblyMapCG object which corresponds to the
2195  * linear space of the current object.
2196  *
2197  * This function is used to create a linear-space assembly map, which is
2198  * then used in the linear space preconditioner in the conjugate
2199  * gradient solve.
2200  */
2202  const ExpList &locexp, GlobalSysSolnType solnType)
2203  {
2204  AssemblyMapCGSharedPtr returnval;
2205 
2206  int i, j;
2207  int nverts = 0;
2208  const boost::shared_ptr<LocalRegions::ExpansionVector> exp
2209  = locexp.GetExp();
2210  int nelmts = exp->size();
2211 
2212  // Get Default Map and turn off any searched values.
2213  returnval = MemoryManager<AssemblyMapCG>
2215  returnval->m_solnType = solnType;
2216  returnval->m_preconType = eNull;
2217  returnval->m_maxStaticCondLevel = 0;
2218  returnval->m_signChange = false;
2219  returnval->m_comm = m_comm;
2220 
2221  // Count the number of vertices
2222  for (i = 0; i < nelmts; ++i)
2223  {
2224  nverts += (*exp)[i]->GetNverts();
2225  }
2226 
2227  returnval->m_numLocalCoeffs = nverts;
2228  returnval->m_localToGlobalMap = Array<OneD, int>(nverts, -1);
2229 
2230  // Store original global ids in this map
2231  returnval->m_localToGlobalBndMap = Array<OneD, int>(nverts, -1);
2232 
2233  int cnt = 0;
2234  int cnt1 = 0;
2235  Array<OneD, int> GlobCoeffs(m_numGlobalCoeffs, -1);
2236 
2237  // Set up local to global map;
2238  for (i = 0; i < nelmts; ++i)
2239  {
2240  for (j = 0; j < (*exp)[i]->GetNverts(); ++j)
2241  {
2242  returnval->m_localToGlobalMap[cnt] =
2243  returnval->m_localToGlobalBndMap[cnt] =
2244  m_localToGlobalMap[cnt1 + (*exp)[i]->GetVertexMap(j,true)];
2245  GlobCoeffs[returnval->m_localToGlobalMap[cnt]] = 1;
2246 
2247  // Set up numLocalDirBndCoeffs
2248  if ((returnval->m_localToGlobalMap[cnt]) <
2250  {
2251  returnval->m_numLocalDirBndCoeffs++;
2252  }
2253  cnt++;
2254  }
2255  cnt1 += (*exp)[i]->GetNcoeffs();
2256  }
2257 
2258  cnt = 0;
2259  // Reset global numbering and count number of dofs
2260  for (i = 0; i < m_numGlobalCoeffs; ++i)
2261  {
2262  if (GlobCoeffs[i] != -1)
2263  {
2264  GlobCoeffs[i] = cnt++;
2265  }
2266  }
2267 
2268  // Set up number of globalCoeffs;
2269  returnval->m_numGlobalCoeffs = cnt;
2270 
2271  // Set up number of global Dirichlet boundary coefficients
2272  for (i = 0; i < m_numGlobalDirBndCoeffs; ++i)
2273  {
2274  if (GlobCoeffs[i] != -1)
2275  {
2276  returnval->m_numGlobalDirBndCoeffs++;
2277  }
2278  }
2279 
2280  // Set up global to universal map
2281  if (m_globalToUniversalMap.num_elements())
2282  {
2284  = m_session->GetComm()->GetRowComm();
2285  int nglocoeffs = returnval->m_numGlobalCoeffs;
2286  returnval->m_globalToUniversalMap
2287  = Array<OneD, int> (nglocoeffs);
2288  returnval->m_globalToUniversalMapUnique
2289  = Array<OneD, int> (nglocoeffs);
2290 
2291  // Reset local to global map and setup universal map
2292  for (i = 0; i < nverts; ++i)
2293  {
2294  cnt = returnval->m_localToGlobalMap[i];
2295  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2296 
2297  returnval->m_globalToUniversalMap[GlobCoeffs[cnt]] =
2299  }
2300 
2301  Nektar::Array<OneD, long> tmp(nglocoeffs);
2302  Vmath::Zero(nglocoeffs, tmp, 1);
2303  for (unsigned int i = 0; i < nglocoeffs; ++i)
2304  {
2305  tmp[i] = returnval->m_globalToUniversalMap[i];
2306  }
2307  returnval->m_gsh = Gs::Init(tmp, vCommRow);
2308  Gs::Unique(tmp, vCommRow);
2309  for (unsigned int i = 0; i < nglocoeffs; ++i)
2310  {
2311  returnval->m_globalToUniversalMapUnique[i]
2312  = (tmp[i] >= 0 ? 1 : 0);
2313  }
2314  }
2315  else // not sure this option is ever needed.
2316  {
2317  for (i = 0; i < nverts; ++i)
2318  {
2319  cnt = returnval->m_localToGlobalMap[i];
2320  returnval->m_localToGlobalMap[i] = GlobCoeffs[cnt];
2321  }
2322  }
2323 
2324  return returnval;
2325  }
2326 
2327  /**
2328  * The bandwidth calculated here corresponds to what is referred to as
2329  * half-bandwidth. If the elements of the matrix are designated as
2330  * a_ij, it corresponds to the maximum value of |i-j| for non-zero
2331  * a_ij. As a result, the value also corresponds to the number of
2332  * sub- or super-diagonals.
2333  *
2334  * The bandwith can be calculated elementally as it corresponds to the
2335  * maximal elemental bandwith (i.e. the maximal difference in global
2336  * DOF index for every element).
2337  *
2338  * We caluclate here the bandwith of the full global system.
2339  */
2341  {
2342  int i,j;
2343  int cnt = 0;
2344  int locSize;
2345  int maxId;
2346  int minId;
2347  int bwidth = -1;
2348  for(i = 0; i < m_numPatches; ++i)
2349  {
2351  maxId = -1;
2352  minId = m_numLocalCoeffs+1;
2353  for(j = 0; j < locSize; j++)
2354  {
2356  {
2357  if(m_localToGlobalMap[cnt+j] > maxId)
2358  {
2359  maxId = m_localToGlobalMap[cnt+j];
2360  }
2361 
2362  if(m_localToGlobalMap[cnt+j] < minId)
2363  {
2364  minId = m_localToGlobalMap[cnt+j];
2365  }
2366  }
2367  }
2368  bwidth = (bwidth>(maxId-minId))?bwidth:(maxId-minId);
2369 
2370  cnt+=locSize;
2371  }
2372 
2373  m_fullSystemBandWidth = bwidth;
2374  }
2375 
2376 
2378  {
2379  return m_localToGlobalMap[i];
2380  }
2381 
2383  {
2384  return m_globalToUniversalMap[i];
2385  }
2386 
2388  {
2389  return m_globalToUniversalMapUnique[i];
2390  }
2391 
2392  const Array<OneD,const int>&
2394  {
2395  return m_localToGlobalMap;
2396  }
2397 
2398  const Array<OneD,const int>&
2400  {
2401  return m_globalToUniversalMap;
2402  }
2403 
2404  const Array<OneD,const int>&
2406  {
2408  }
2409 
2411  const int i) const
2412  {
2413  if(m_signChange)
2414  {
2415  return m_localToGlobalSign[i];
2416  }
2417  else
2418  {
2419  return 1.0;
2420  }
2421  }
2422 
2424  {
2425  return m_localToGlobalSign;
2426  }
2427 
2429  const Array<OneD, const NekDouble>& loc,
2430  Array<OneD, NekDouble>& global) const
2431  {
2433  if(global.data() == loc.data())
2434  {
2435  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2436  }
2437  else
2438  {
2439  local = loc; // create reference
2440  }
2441 
2442 
2443  if(m_signChange)
2444  {
2445  Vmath::Scatr(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2446  }
2447  else
2448  {
2449  Vmath::Scatr(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2450  }
2451 
2452  // ensure all values are unique by calling a max
2453  Gs::Gather(global, Gs::gs_max, m_gsh);
2454  }
2455 
2457  const NekVector<NekDouble>& loc,
2458  NekVector< NekDouble>& global) const
2459  {
2460  LocalToGlobal(loc.GetPtr(),global.GetPtr());
2461  }
2462 
2464  const Array<OneD, const NekDouble>& global,
2465  Array<OneD, NekDouble>& loc) const
2466  {
2468  if(global.data() == loc.data())
2469  {
2470  glo = Array<OneD, NekDouble>(global.num_elements(),global.data());
2471  }
2472  else
2473  {
2474  glo = global; // create reference
2475  }
2476 
2477 
2478  if(m_signChange)
2479  {
2480  Vmath::Gathr(m_numLocalCoeffs, m_localToGlobalSign.get(), glo.get(), m_localToGlobalMap.get(), loc.get());
2481  }
2482  else
2483  {
2484  Vmath::Gathr(m_numLocalCoeffs, glo.get(), m_localToGlobalMap.get(), loc.get());
2485  }
2486  }
2487 
2489  const NekVector<NekDouble>& global,
2490  NekVector< NekDouble>& loc) const
2491  {
2492  GlobalToLocal(global.GetPtr(),loc.GetPtr());
2493  }
2494 
2496  const Array<OneD, const NekDouble> &loc,
2497  Array<OneD, NekDouble> &global) const
2498  {
2500  if(global.data() == loc.data())
2501  {
2502  local = Array<OneD, NekDouble>(loc.num_elements(),loc.data());
2503  }
2504  else
2505  {
2506  local = loc; // create reference
2507  }
2508  //ASSERTL1(loc.get() != global.get(),"Local and Global Arrays cannot be the same");
2509 
2510  Vmath::Zero(m_numGlobalCoeffs, global.get(), 1);
2511 
2512  if(m_signChange)
2513  {
2514  Vmath::Assmb(m_numLocalCoeffs, m_localToGlobalSign.get(), local.get(), m_localToGlobalMap.get(), global.get());
2515  }
2516  else
2517  {
2518  Vmath::Assmb(m_numLocalCoeffs, local.get(), m_localToGlobalMap.get(), global.get());
2519  }
2520  UniversalAssemble(global);
2521  }
2522 
2524  const NekVector<NekDouble>& loc,
2525  NekVector< NekDouble>& global) const
2526  {
2527  Assemble(loc.GetPtr(),global.GetPtr());
2528  }
2529 
2531  Array<OneD, NekDouble>& pGlobal) const
2532  {
2533  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
2534  }
2535 
2537  NekVector< NekDouble>& pGlobal) const
2538  {
2539  UniversalAssemble(pGlobal.GetPtr());
2540  }
2541 
2543  Array<OneD, NekDouble>& pGlobal,
2544  int offset) const
2545  {
2546  Array<OneD, NekDouble> tmp(offset);
2547  Vmath::Vcopy(offset, pGlobal, 1, tmp, 1);
2548  UniversalAssemble(pGlobal);
2549  Vmath::Vcopy(offset, tmp, 1, pGlobal, 1);
2550  }
2551 
2553  {
2554  return m_fullSystemBandWidth;
2555  }
2556 
2558  {
2559  return m_numNonDirVertexModes;
2560  }
2561 
2563  {
2564  return m_numNonDirEdgeModes;
2565  }
2566 
2568  {
2569  return m_numNonDirFaceModes;
2570  }
2571 
2573  {
2574  return m_numDirEdges;
2575  }
2576 
2578  {
2579  return m_numDirFaces;
2580  }
2581 
2583  {
2584  return m_numNonDirEdges;
2585  }
2586 
2588  {
2589  return m_numNonDirFaces;
2590  }
2591 
2593  {
2594  return m_extraDirEdges;
2595  }
2596  } // namespace
2597 } // namespace
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:188
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:218
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)