Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AssemblyMapDG.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AssemblyMapDG.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: Local to Global Base Class mapping routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <MultiRegions/ExpList.h>
38 #include <LocalRegions/SegExp.h>
39 #include <LocalRegions/TriExp.h>
40 #include <LocalRegions/QuadExp.h>
41 #include <LocalRegions/HexExp.h>
42 #include <LocalRegions/TetExp.h>
43 #include <LocalRegions/PrismExp.h>
44 #include <LocalRegions/PyrExp.h>
45 #include <LocalRegions/PointExp.h>
48 
49 #include <boost/config.hpp>
50 #include <boost/graph/adjacency_list.hpp>
51 #include <boost/graph/cuthill_mckee_ordering.hpp>
52 #include <boost/graph/properties.hpp>
53 #include <boost/graph/bandwidth.hpp>
54 
55 using namespace std;
56 
57 namespace Nektar
58 {
59  namespace MultiRegions
60  {
61  AssemblyMapDG::AssemblyMapDG():
62  m_numDirichletBndPhys(0)
63  {
64  }
65 
67  {
68  }
69 
73  const ExpListSharedPtr &trace,
74  const ExpList &locExp,
77  const PeriodicMap &periodicTrace,
78  const std::string variable):
79  AssemblyMap(pSession,variable)
80  {
81  int i, j, k, cnt, eid, id, id1, gid;
82  int order_e = 0;
83  int nTraceExp = trace->GetExpSize();
84  int nbnd = bndCondExp.num_elements();
85 
89 
90  const LocalRegions::ExpansionVector expList = *(locExp.GetExp());
91  int nel = expList.size();
92 
93  map<int, int> meshTraceId;
94 
95  m_signChange = true;
96 
97  // determine mapping from geometry edges to trace
98  for(i = 0; i < nTraceExp; ++i)
99  {
100  meshTraceId[trace->GetExp(i)->GetGeom()->GetGlobalID()] = i;
101  }
102 
103  // Count total number of trace elements
104  cnt = 0;
105  for(i = 0; i < nel; ++i)
106  {
107  cnt += expList[i]->GetNtrace();
108  }
109 
113 
114  // set up trace expansions links;
115  for(cnt = i = 0; i < nel; ++i)
116  {
117  m_elmtToTrace[i] = tracemap + cnt;
118 
119  for(j = 0; j < expList[i]->GetNtrace(); ++j)
120  {
121  id = expList[i]->GetGeom()->GetTid(j);
122 
123  if(meshTraceId.count(id) > 0)
124  {
125  m_elmtToTrace[i][j] =
126  trace->GetExp(meshTraceId.find(id)->second);
127  }
128  else
129  {
130  ASSERTL0(false, "Failed to find trace map");
131  }
132  }
133 
134  cnt += expList[i]->GetNtrace();
135  }
136 
137  // Set up boundary mapping
138  cnt = 0;
139  for(i = 0; i < nbnd; ++i)
140  {
141  cnt += bndCondExp[i]->GetExpSize();
142  }
143 
144  set<int> dirTrace;
145 
148 
149  cnt = 0;
150  for(i = 0; i < bndCondExp.num_elements(); ++i)
151  {
152  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
153  {
154  bndExp = bndCondExp[i]->GetExp(j);
155  traceGeom = bndExp->GetGeom();
156  id = traceGeom->GetGlobalID();
157 
158  if(bndCond[i]->GetBoundaryConditionType() ==
160  {
161  m_numLocalDirBndCoeffs += bndExp->GetNcoeffs();
162  m_numDirichletBndPhys += bndExp->GetTotPoints();
163  dirTrace.insert(id);
164  }
165  }
166 
167  cnt += j;
168  }
169 
170  // Set up integer mapping array and sign change for each degree of
171  // freedom + initialise some more data members.
172  m_staticCondLevel = 0;
174  m_numPatches = nel;
177 
178  int nbndry = 0;
179  for(i = 0; i < nel; ++i) // count number of elements in array
180  {
181  eid = locExp.GetOffset_Elmt_Id(i);
182  nbndry += expList[eid]->NumDGBndryCoeffs();
183  m_numLocalIntCoeffsPerPatch[i] = 0;
185  (unsigned int) expList[eid]->NumDGBndryCoeffs();
186  }
187 
189  m_numLocalBndCoeffs = nbndry;
190  m_numLocalCoeffs = nbndry;
193 
194  // Set up array for potential mesh optimsation
195  Array<OneD,int> traceElmtGid(nTraceExp, -1);
196  int nDir = 0;
197  cnt = 0;
198 
199  // We are now going to construct a graph of the mesh which can be
200  // reordered depending on the type of solver we would like to use.
201  typedef boost::adjacency_list<
202  boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
203 
204  BoostGraph boostGraphObj;
205  int trace_id, trace_id1;
206  int dirOffset = 0;
207 
208  // make trace trace renumbering map where first solved trace starts
209  // at 0 so we can set up graph.
210  for(i = 0; i < nTraceExp; ++i)
211  {
212  id = trace->GetExp(i)->GetGeom()->GetGlobalID();
213 
214  if (dirTrace.count(id) == 0)
215  {
216  // Initial put in element ordering (starting from zero) into
217  // traceElmtGid
218  boost::add_vertex(boostGraphObj);
219  traceElmtGid[i] = cnt++;
220  }
221  else
222  {
223  // Use existing offset for Dirichlet edges
224  traceElmtGid[i] = dirOffset;
225  dirOffset += trace->GetExp(i)->GetNcoeffs();
226  nDir++;
227  }
228  }
229 
230  // Set up boost Graph
231  for(i = 0; i < nel; ++i)
232  {
233  eid = locExp.GetOffset_Elmt_Id(i);
234 
235  for(j = 0; j < expList[eid]->GetNtrace(); ++j)
236  {
237  // Add trace to boost graph for non-Dirichlet Boundary
238  traceGeom = m_elmtToTrace[eid][j]->GetGeom();
239  id = traceGeom->GetGlobalID();
240  trace_id = meshTraceId.find(id)->second;
241 
242  if(dirTrace.count(id) == 0)
243  {
244  for(k = j+1; k < expList[eid]->GetNtrace(); ++k)
245  {
246  traceGeom = m_elmtToTrace[eid][k]->GetGeom();
247  id1 = traceGeom->GetGlobalID();
248  trace_id1 = meshTraceId.find(id1)->second;
249 
250  if(dirTrace.count(id1) == 0)
251  {
252  boost::add_edge((size_t)traceElmtGid[trace_id],
253  (size_t)traceElmtGid[trace_id1],
254  boostGraphObj);
255  }
256  }
257  }
258  }
259  }
260 
261  int nGraphVerts = nTraceExp - nDir;
262  Array<OneD, int> perm (nGraphVerts);
263  Array<OneD, int> iperm(nGraphVerts);
264  Array<OneD, int> vwgts(nGraphVerts);
266 
267  for(i = 0; i < nGraphVerts; ++i)
268  {
269  vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
270  }
271 
272  if(nGraphVerts)
273  {
274  switch(m_solnType)
275  {
276  case eDirectFullMatrix:
277  case eIterativeFull:
279  case eXxtFullMatrix:
280  case eXxtStaticCond:
281  case ePETScFullMatrix:
282  case ePETScStaticCond:
283  {
284  NoReordering(boostGraphObj,perm,iperm);
285  break;
286  }
287  case eDirectStaticCond:
288  {
289  CuthillMckeeReordering(boostGraphObj,perm,iperm);
290  break;
291  }
296  {
297  MultiLevelBisectionReordering(boostGraphObj,perm,iperm,
298  bottomUpGraph);
299  break;
300  }
301  default:
302  {
303  ASSERTL0(false,"Unrecognised solution type");
304  }
305  }
306  }
307 
308  // Recast the permutation so that it can be used as a map from old
309  // trace ID to new trace ID
311  for(i = 0; i < nTraceExp - nDir; ++i)
312  {
313  traceElmtGid[perm[i]+nDir] = cnt;
314  cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
315  }
316 
317  // Now have trace edges Gid position
318 
319  cnt = 0;
320  for(i = 0; i < nel; ++i)
321  {
322  // order list according to m_offset_elmt_id details in expList
323  eid = locExp.GetOffset_Elmt_Id(i);
324  exp = expList[eid];
325 
326  for(j = 0; j < exp->GetNtrace(); ++j)
327  {
328  traceGeom = m_elmtToTrace[eid][j]->GetGeom();
329  id = traceGeom->GetGlobalID();
330  gid = traceElmtGid[meshTraceId.find(id)->second];
331 
332  const int nDim = expList[eid]->GetNumBases();
333 
334  if (nDim == 1)
335  {
336  order_e = 1;
337  m_localToGlobalBndMap[cnt] = gid;
338  }
339  else if (nDim == 2)
340  {
341  order_e = expList[eid]->GetEdgeNcoeffs(j);
342 
343  if(expList[eid]->GetEorient(j) == StdRegions::eForwards)
344  {
345  for(k = 0; k < order_e; ++k)
346  {
347  m_localToGlobalBndMap[k+cnt] = gid + k;
348  }
349  }
350  else
351  {
352  switch(m_elmtToTrace[eid][j]->GetBasisType(0))
353  {
355  {
356  // reverse vertex order
357  m_localToGlobalBndMap[cnt] = gid + 1;
358  m_localToGlobalBndMap[cnt+1] = gid;
359  for (k = 2; k < order_e; ++k)
360  {
361  m_localToGlobalBndMap[k+cnt] = gid + k;
362  }
363 
364  // negate odd modes
365  for(k = 3; k < order_e; k+=2)
366  {
367  m_localToGlobalBndSign[cnt+k] = -1.0;
368  }
369  break;
370  }
372  {
373  // reverse order
374  for(k = 0; k < order_e; ++k)
375  {
376  m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
377  }
378  break;
379  }
381  {
382  // reverse order
383  for(k = 0; k < order_e; ++k)
384  {
385  m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
386  }
387  break;
388  }
389  default:
390  {
391  ASSERTL0(false,"Boundary type not permitted");
392  }
393  }
394  }
395  }
396  else if (nDim == 3)
397  {
398  order_e = expList[eid]->GetFaceNcoeffs(j);
399 
400  std::map<int, int> orientMap;
401 
402  Array<OneD, unsigned int> elmMap1 (order_e);
403  Array<OneD, int> elmSign1(order_e);
404  Array<OneD, unsigned int> elmMap2 (order_e);
405  Array<OneD, int> elmSign2(order_e);
406 
407  StdRegions::Orientation fo = expList[eid]->GetForient(j);
408 
409  // Construct mapping which will permute global IDs
410  // according to face orientations.
411  expList[eid]->GetFaceToElementMap(j,fo,elmMap1,elmSign1);
412  expList[eid]->GetFaceToElementMap(
413  j,StdRegions::eDir1FwdDir1_Dir2FwdDir2,elmMap2,elmSign2);
414 
415  for (k = 0; k < elmMap1.num_elements(); ++k)
416  {
417  // Find the elemental co-efficient in the original
418  // mapping.
419  int idx = -1;
420  for (int l = 0; l < elmMap2.num_elements(); ++l)
421  {
422  if (elmMap1[k] == elmMap2[l])
423  {
424  idx = l;
425  break;
426  }
427  }
428 
429  ASSERTL2(idx != -1, "Problem with face to element map!");
430  orientMap[k] = idx;
431  }
432 
433  for(k = 0; k < order_e; ++k)
434  {
435  m_localToGlobalBndMap [k+cnt] = gid + orientMap[k];
436  m_localToGlobalBndSign[k+cnt] = elmSign2[orientMap[k]];
437  }
438  }
439 
440  cnt += order_e;
441  }
442  }
443 
444  // set up m_bndCondCoeffsToGlobalCoeffsMap to align with map
445  cnt = 0;
446  for(i = 0; i < nbnd; ++i)
447  {
448  cnt += bndCondExp[i]->GetNcoeffs();
449  }
450 
452 
453  // Number of boundary expansions
454  int nbndexp = 0, bndOffset, bndTotal = 0;
455  for(cnt = i = 0; i < nbnd; ++i)
456  {
457  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
458  {
459  bndExp = bndCondExp[i]->GetExp(j);
460  id = bndExp->GetGeom()->GetGlobalID();
461  gid = traceElmtGid[meshTraceId.find(id)->second];
462  bndOffset = bndCondExp[i]->GetCoeff_Offset(j) + bndTotal;
463 
464  // Since boundary information is defined to be aligned with
465  // the geometry just use forward/forward (both coordinate
466  // directions) defintiion for gid's
467  for(k = 0; k < bndExp->GetNcoeffs(); ++k)
468  {
469  m_bndCondCoeffsToGlobalCoeffsMap[bndOffset+k] = gid + k;
470  }
471  }
472 
473  nbndexp += bndCondExp[i]->GetExpSize();
474  bndTotal += bndCondExp[i]->GetNcoeffs();
475  }
476 
477  m_numGlobalBndCoeffs = trace->GetNcoeffs();
479 
481 
486  && nGraphVerts)
487  {
488  if (m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
489  {
490  Array<OneD, int> vwgts_perm(nGraphVerts);
491 
492  for(int i = 0; i < nGraphVerts; i++)
493  {
494  vwgts_perm[i] = vwgts[perm[i]];
495  }
496 
497  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
499  AllocateSharedPtr(this, bottomUpGraph);
500  }
501  }
502 
503  cnt = 0;
505  for(i = 0; i < bndCondExp.num_elements(); ++i)
506  {
507  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
508  {
509  bndExp = bndCondExp[i]->GetExp(j);
510  traceGeom = bndExp->GetGeom();
511  id = traceGeom->GetGlobalID();
513  meshTraceId.find(id)->second;
514  }
515  }
516 
517  // Now set up mapping from global coefficients to universal.
518  ExpListSharedPtr tr = boost::dynamic_pointer_cast<ExpList>(trace);
519  SetUpUniversalDGMap (locExp);
520  SetUpUniversalTraceMap(locExp, tr, periodicTrace);
521 
522  m_hash = boost::hash_range(m_localToGlobalBndMap.begin(),
523  m_localToGlobalBndMap.end());
524  }
525 
526  /**
527  * Constructs a mapping between the process-local global numbering and
528  * a universal numbering of the trace space expansion. The universal
529  * numbering is defined by the mesh edge IDs to enforce consistency
530  * across processes.
531  *
532  * @param locExp List of local elemental expansions.
533  */
535  {
537  int eid = 0;
538  int cnt = 0;
539  int id = 0;
540  int order_e = 0;
541  int vGlobalId = 0;
542  int maxDof = 0;
543  int dof = 0;
544  int nDim = 0;
545  int i,j,k;
546 
547  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
548 
549  // Initialise the global to universal maps.
552 
553  // Loop over all the elements in the domain and compute max
554  // DOF. Reduce across all processes to get universal maximum.
555  for(i = 0; i < locExpVector.size(); ++i)
556  {
557  locExpansion = locExpVector[i];
558  nDim = locExpansion->GetShapeDimension();
559 
560  // Loop over all edges of element i
561  if (nDim == 1)
562  {
563  maxDof = (1 > maxDof ? 1 : maxDof);
564  }
565  else if (nDim == 2)
566  {
567  for (j = 0; j < locExpansion->GetNedges(); ++j)
568  {
569  dof = locExpansion->GetEdgeNcoeffs(j);
570  maxDof = (dof > maxDof ? dof : maxDof);
571  }
572  }
573  else if (nDim == 3)
574  {
575  for (j = 0; j < locExpansion->GetNfaces(); ++j)
576  {
577  dof = locExpansion->GetFaceNcoeffs(j);
578  maxDof = (dof > maxDof ? dof : maxDof);
579  }
580  }
581  }
582  m_comm->AllReduce(maxDof, LibUtilities::ReduceMax);
583 
584  // Now have trace edges Gid position
585  cnt = 0;
586  for(i = 0; i < locExpVector.size(); ++i)
587  {
588  locExpansion = locExpVector[i];
589  nDim = locExpansion->GetShapeDimension();
590 
591  // Order list according to m_offset_elmt_id details in Exp
592  // so that triangules are listed first and then quads
593  eid = locExp.GetOffset_Elmt_Id(i);
594 
595  // Populate mapping for each edge of the element.
596  if (nDim == 1)
597  {
598  int nverts = locExpansion->GetNverts();
599  for(j = 0; j < nverts; ++j)
600  {
601  LocalRegions::PointExpSharedPtr locPointExp =
602  m_elmtToTrace[eid][j]->as<LocalRegions::PointExp>();
603  id = locPointExp->GetGeom()->GetGlobalID();
604  vGlobalId = m_localToGlobalBndMap[cnt+j];
605  m_globalToUniversalBndMap[vGlobalId]
606  = id * maxDof + j + 1;
607  }
608  cnt += nverts;
609  }
610  else if (nDim == 2)
611  {
612  for(j = 0; j < locExpansion->GetNedges(); ++j)
613  {
615  m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
616 
617  id = locSegExp->GetGeom1D()->GetEid();
618  order_e = locExpVector[eid]->GetEdgeNcoeffs(j);
619 
620  map<int,int> orientMap;
621  Array<OneD, unsigned int> map1(order_e), map2(order_e);
622  Array<OneD, int> sign1(order_e), sign2(order_e);
623 
624  locExpVector[eid]->GetEdgeToElementMap(j, StdRegions::eForwards, map1, sign1);
625  locExpVector[eid]->GetEdgeToElementMap(j, locExpVector[eid]->GetEorient(j), map2, sign2);
626 
627  for (k = 0; k < map1.num_elements(); ++k)
628  {
629  // Find the elemental co-efficient in the original
630  // mapping.
631  int idx = -1;
632  for (int l = 0; l < map2.num_elements(); ++l)
633  {
634  if (map1[k] == map2[l])
635  {
636  idx = l;
637  break;
638  }
639  }
640 
641  ASSERTL2(idx != -1, "Problem with face to element map!");
642  orientMap[k] = idx;
643  }
644 
645  for(k = 0; k < order_e; ++k)
646  {
647  vGlobalId = m_localToGlobalBndMap[k+cnt];
648  m_globalToUniversalBndMap[vGlobalId]
649  = id * maxDof + orientMap[k] + 1;
650  }
651  cnt += order_e;
652  }
653  }
654  else if (nDim == 3)
655  {
656  for(j = 0; j < locExpansion->GetNfaces(); ++j)
657  {
659  m_elmtToTrace[eid][j]
661 
662  id = locFaceExp->GetGeom2D()->GetFid();
663  order_e = locExpVector[eid]->GetFaceNcoeffs(j);
664 
665  map<int,int> orientMap;
666  Array<OneD, unsigned int> map1(order_e), map2(order_e);
667  Array<OneD, int> sign1(order_e), sign2(order_e);
668 
669  locExpVector[eid]->GetFaceToElementMap(j, StdRegions::eDir1FwdDir1_Dir2FwdDir2, map1, sign1);
670  locExpVector[eid]->GetFaceToElementMap(j, locExpVector[eid]->GetForient(j), map2, sign2);
671 
672  for (k = 0; k < map1.num_elements(); ++k)
673  {
674  // Find the elemental co-efficient in the original
675  // mapping.
676  int idx = -1;
677  for (int l = 0; l < map2.num_elements(); ++l)
678  {
679  if (map1[k] == map2[l])
680  {
681  idx = l;
682  break;
683  }
684  }
685 
686  ASSERTL2(idx != -1, "Problem with face to element map!");
687  orientMap[k] = idx;
688  }
689 
690  for(k = 0; k < order_e; ++k)
691  {
692  vGlobalId = m_localToGlobalBndMap[k+cnt];
693  m_globalToUniversalBndMap[vGlobalId]
694  = id * maxDof + orientMap[k] + 1;
695  }
696  cnt += order_e;
697  }
698  }
699  }
700 
701  // Initialise GSlib and populate the unique map.
702  Array<OneD, long> tmp(m_globalToUniversalBndMap.num_elements());
703  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
704  {
705  tmp[i] = m_globalToUniversalBndMap[i];
706  }
707  m_bndGsh = m_gsh = Gs::Init(tmp, m_comm);
708  Gs::Unique(tmp, m_comm);
709  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
710  {
711  m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
712  }
713  }
714 
716  const ExpList &locExp,
717  const ExpListSharedPtr trace,
718  const PeriodicMap &perMap)
719  {
720  Array<OneD, int> tmp;
722  int i;
723  int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
724 
725  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
726 
727  int nTracePhys = trace->GetTotPoints();
728 
729  // Initialise the trace to universal maps.
731  Nektar::Array<OneD, int>(nTracePhys, -1);
733  Nektar::Array<OneD, int>(nTracePhys, -1);
734 
735  // Assume that each element of the expansion is of the same
736  // dimension.
737  nDim = locExpVector[0]->GetShapeDimension();
738 
739  if (nDim == 1)
740  {
741  maxQuad = (1 > maxQuad ? 1 : maxQuad);
742  }
743  else
744  {
745  for (i = 0; i < trace->GetExpSize(); ++i)
746  {
747  quad = trace->GetExp(i)->GetTotPoints();
748  if (quad > maxQuad)
749  {
750  maxQuad = quad;
751  }
752  }
753  }
754  m_comm->AllReduce(maxQuad, LibUtilities::ReduceMax);
755 
756  if (nDim == 1)
757  {
758  for (int i = 0; i < trace->GetExpSize(); ++i)
759  {
760  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
761  offset = trace->GetPhys_Offset(i);
762 
763  // Check to see if this vert is periodic. If it is, then we
764  // need use the unique eid of the two points
765  PeriodicMap::const_iterator it = perMap.find(eid);
766  if (perMap.count(eid) > 0)
767  {
768  PeriodicEntity ent = it->second[0];
769  if (ent.isLocal == false) // Not sure if true in 1D
770  {
771  eid = min(eid, ent.id);
772  }
773  }
774 
775  m_traceToUniversalMap[offset] = eid*maxQuad+1;
776  }
777  }
778  else
779  {
780  for (int i = 0; i < trace->GetExpSize(); ++i)
781  {
782  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
783  offset = trace->GetPhys_Offset(i);
784  quad = trace->GetExp(i)->GetTotPoints();
785 
786  // Check to see if this edge is periodic. If it is, then we
787  // need to reverse the trace order of one edge only in the
788  // universal map so that the data are reversed w.r.t each
789  // other. We do this by using the minimum of the two IDs.
790  PeriodicMap::const_iterator it = perMap.find(eid);
791  bool realign = false;
792  if (perMap.count(eid) > 0)
793  {
794  PeriodicEntity ent = it->second[0];
795  if (ent.isLocal == false)
796  {
797  realign = eid == min(eid, ent.id);
798  eid = min(eid, ent.id);
799  }
800  }
801 
802  for (int j = 0; j < quad; ++j)
803  {
804  m_traceToUniversalMap[j+offset] = eid*maxQuad+j+1;
805  }
806 
807  if (realign)
808  {
809  if (nDim == 2)
810  {
812  tmp = m_traceToUniversalMap+offset,
813  it->second[0].orient, quad);
814  }
815  else
816  {
818  tmp = m_traceToUniversalMap+offset,
819  it->second[0].orient,
820  trace->GetExp(i)->GetNumPoints(0),
821  trace->GetExp(i)->GetNumPoints(1));
822  }
823  }
824  }
825  }
826 
827  Array<OneD, long> tmp2(nTracePhys);
828  for (int i = 0; i < nTracePhys; ++i)
829  {
830  tmp2[i] = m_traceToUniversalMap[i];
831  }
832  m_traceGsh = Gs::Init(tmp2, m_comm);
833  Gs::Unique(tmp2, m_comm);
834  for (int i = 0; i < nTracePhys; ++i)
835  {
836  m_traceToUniversalMapUnique[i] = tmp2[i];
837  }
838  }
839 
841  Array<OneD, int> &toAlign,
843  int nquad1,
844  int nquad2)
845  {
846  if (orient == StdRegions::eBackwards)
847  {
848  ASSERTL1(nquad2 == 0, "nquad2 != 0 for reorienation");
849  for (int i = 0; i < nquad1/2; ++i)
850  {
851  swap(toAlign[i], toAlign[nquad1-1-i]);
852  }
853  }
854  else if (orient != StdRegions::eForwards)
855  {
856  ASSERTL1(nquad2 != 0, "nquad2 == 0 for reorienation");
857 
858  Array<OneD, int> tmp(nquad1*nquad2);
859 
860  // Copy transpose.
861  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
865  {
866  for (int i = 0; i < nquad2; ++i)
867  {
868  for (int j = 0; j < nquad1; ++j)
869  {
870  tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
871  }
872  }
873  }
874  else
875  {
876  for (int i = 0; i < nquad2; ++i)
877  {
878  for (int j = 0; j < nquad1; ++j)
879  {
880  tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
881  }
882  }
883  }
884 
885  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
889  {
890  // Reverse x direction
891  for (int i = 0; i < nquad2; ++i)
892  {
893  for (int j = 0; j < nquad1/2; ++j)
894  {
895  swap(tmp[i*nquad1 + j],
896  tmp[i*nquad1 + nquad1-j-1]);
897  }
898  }
899  }
900 
901  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
905  {
906  // Reverse y direction
907  for (int j = 0; j < nquad1; ++j)
908  {
909  for (int i = 0; i < nquad2/2; ++i)
910  {
911  swap(tmp[i*nquad1 + j],
912  tmp[(nquad2-i-1)*nquad1 + j]);
913  }
914  }
915  }
916  Vmath::Vcopy(nquad1*nquad2, tmp, 1, toAlign, 1);
917  }
918  }
919 
921  Array<OneD, NekDouble> &pGlobal) const
922  {
923  Gs::Gather(pGlobal, Gs::gs_add, m_traceGsh);
924  }
925 
927  {
928  return m_localToGlobalBndMap[i];
929  }
930 
932  {
933  return m_globalToUniversalBndMap[i];
934  }
935 
937  {
939  }
940 
942  {
943  return m_localToGlobalBndMap;
944  }
945 
947  {
949  }
950 
952  {
954  }
955 
957  const int i) const
958  {
959  return GetLocalToGlobalBndSign(i);
960  }
961 
963  const Array<OneD, const NekDouble>& loc,
964  Array<OneD, NekDouble>& global) const
965  {
966  AssembleBnd(loc,global);
967  }
968 
970  const NekVector<NekDouble>& loc,
971  NekVector< NekDouble>& global) const
972  {
973  AssembleBnd(loc,global);
974  }
975 
977  const Array<OneD, const NekDouble>& global,
978  Array<OneD, NekDouble>& loc) const
979  {
980  GlobalToLocalBnd(global,loc);
981  }
982 
984  const NekVector<NekDouble>& global,
985  NekVector< NekDouble>& loc) const
986  {
987  GlobalToLocalBnd(global,loc);
988  }
989 
991  const Array<OneD, const NekDouble> &loc,
992  Array<OneD, NekDouble> &global) const
993  {
994  AssembleBnd(loc,global);
995  }
996 
998  const NekVector<NekDouble>& loc,
999  NekVector< NekDouble>& global) const
1000  {
1001  AssembleBnd(loc,global);
1002  }
1003 
1005  Array<OneD, NekDouble>& pGlobal) const
1006  {
1007  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
1008  }
1009 
1011  NekVector< NekDouble>& pGlobal) const
1012  {
1013  UniversalAssemble(pGlobal.GetPtr());
1014  }
1015 
1017  {
1018  return GetBndSystemBandWidth();
1019  }
1020 
1022  {
1023  return m_traceToUniversalMap[i];
1024  }
1025 
1027  {
1028  return m_traceToUniversalMapUnique[i];
1029  }
1030 
1032  {
1033  return m_numDirichletBndPhys;
1034  }
1035 
1038  {
1039  ASSERTL1(i >= 0 && i < m_elmtToTrace.num_elements(),
1040  "i is out of range");
1041  return m_elmtToTrace[i];
1042  }
1043 
1046  {
1047  return m_elmtToTrace;
1048  }
1049  } //namespace
1050 } // namespace
AssemblyMapDG()
Default constructor.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:345
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:314
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
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_Assemble(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
Array< OneD, int > m_bndCondTraceToGlobalTraceMap
Integer map of bnd cond trace number to global trace number.
Definition: AssemblyMap.h:356
bool isLocal
Flag specifying if this entity is local to this partition.
int GetBndSystemBandWidth() const
Returns the bandwidth of the boundary system.
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
Principle Modified Functions .
Definition: BasisType.h:49
boost::shared_ptr< BottomUpSubStructuredGraph > BottomUpSubStructuredGraphSharedPtr
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
STL namespace.
Array< OneD, const NekDouble > GetLocalToGlobalBndSign() const
Retrieve the sign change for all local boundary modes.
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)
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
Array< OneD, int > m_traceToUniversalMap
Integer map of process trace space quadrature points to universal space.
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global) const
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
void UniversalTraceAssemble(Array< OneD, NekDouble > &pGlobal) const
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:309
virtual void v_UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
const SpatialDomains::PointGeomSharedPtr GetGeom() const
Definition: PointExp.h:111
void UniversalAssemble(Array< OneD, NekDouble > &pGlobal) const
boost::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:266
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > & GetElmtToTrace()
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:269
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMapUnique()
void SetUpUniversalTraceMap(const ExpList &locExp, const ExpListSharedPtr trace, const PeriodicMap &perMap=NullPeriodicMap)
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
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
int id
Geometry ID of entity.
virtual const Array< OneD, const int > & v_GetLocalToGlobalMap()
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
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
boost::shared_ptr< PointExp > PointExpSharedPtr
Definition: PointExp.h:131
void SetUpUniversalDGMap(const ExpList &locExp)
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Array< OneD, unsigned int > m_numLocalIntCoeffsPerPatch
The number of int dofs per patch.
Definition: AssemblyMap.h:390
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
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
virtual ~AssemblyMapDG()
Destructor.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:316
Array< OneD, int > m_traceToUniversalMapUnique
Integer map of unique process trace space quadrature points to universal space (signed).
Array< OneD, int > m_bndCondCoeffsToGlobalCoeffsMap
Integer map of bnd cond coeffs to global coefficients.
Definition: AssemblyMap.h:352
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:312
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
Array< OneD, Array< OneD, LocalRegions::ExpansionSharedPtr > > m_elmtToTrace
list of edge expansions for a given element
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
virtual void v_GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
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
void AssembleBnd(const NekVector< NekDouble > &loc, NekVector< NekDouble > &global, int offset) const
virtual const Array< OneD, NekDouble > & v_GetLocalToGlobalSign() const
virtual const Array< OneD, const int > & v_GetGlobalToUniversalMap()
SpatialDomains::Geometry1DSharedPtr GetGeom1D() const
Definition: Expansion1D.h:169
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
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:342
void GlobalToLocalBnd(const NekVector< NekDouble > &global, NekVector< NekDouble > &loc, int offset) const
Lagrange for SEM basis .
Definition: BasisType.h:53
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
boost::shared_ptr< Geometry > GeometrySharedPtr
Definition: Geometry.h:53
void RealignTraceElement(Array< OneD, int > &toAlign, StdRegions::Orientation orient, int nquad1, int nquad2=0)
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
int GetNumDirichletBndPhys()
Return the number of boundary segments on which Dirichlet boundary conditions are imposed...
int m_numPatches
The number of patches (~elements) in the current level.
Definition: AssemblyMap.h:386
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49
void NoReordering(const BoostGraph &graph, Array< OneD, int > &perm, Array< OneD, int > &iperm)
int m_numDirichletBndPhys
Number of physical dirichlet boundary values in trace.
Definition: AssemblyMapDG.h:97
virtual int v_GetFullSystemBandWidth() const