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