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