Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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 = 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 = 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  eid = 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  eid = i;
588  locExpansion = locExpVector[eid];
589  nDim = locExpansion->GetShapeDimension();
590 
591  // Populate mapping for each edge of the element.
592  if (nDim == 1)
593  {
594  int nverts = locExpansion->GetNverts();
595  for(j = 0; j < nverts; ++j)
596  {
597  LocalRegions::PointExpSharedPtr locPointExp =
598  m_elmtToTrace[eid][j]->as<LocalRegions::PointExp>();
599  id = locPointExp->GetGeom()->GetGlobalID();
600  vGlobalId = m_localToGlobalBndMap[cnt+j];
601  m_globalToUniversalBndMap[vGlobalId]
602  = id * maxDof + j + 1;
603  }
604  cnt += nverts;
605  }
606  else if (nDim == 2)
607  {
608  for(j = 0; j < locExpansion->GetNedges(); ++j)
609  {
611  m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
612 
613  id = locSegExp->GetGeom1D()->GetEid();
614  order_e = locExpansion->GetEdgeNcoeffs(j);
615 
616  map<int,int> orientMap;
617  Array<OneD, unsigned int> map1(order_e), map2(order_e);
618  Array<OneD, int> sign1(order_e), sign2(order_e);
619 
620  locExpansion->GetEdgeToElementMap(j, StdRegions::eForwards, map1, sign1);
621  locExpansion->GetEdgeToElementMap(j, locExpansion->GetEorient(j), map2, sign2);
622 
623  for (k = 0; k < map1.num_elements(); ++k)
624  {
625  // Find the elemental co-efficient in the original
626  // mapping.
627  int idx = -1;
628  for (int l = 0; l < map2.num_elements(); ++l)
629  {
630  if (map1[k] == map2[l])
631  {
632  idx = l;
633  break;
634  }
635  }
636 
637  ASSERTL2(idx != -1, "Problem with face to element map!");
638  orientMap[k] = idx;
639  }
640 
641  for(k = 0; k < order_e; ++k)
642  {
643  vGlobalId = m_localToGlobalBndMap[k+cnt];
644  m_globalToUniversalBndMap[vGlobalId]
645  = id * maxDof + orientMap[k] + 1;
646  }
647  cnt += order_e;
648  }
649  }
650  else if (nDim == 3)
651  {
652  for(j = 0; j < locExpansion->GetNfaces(); ++j)
653  {
655  m_elmtToTrace[eid][j]
657 
658  id = locFaceExp->GetGeom2D()->GetFid();
659  order_e = locExpansion->GetFaceNcoeffs(j);
660 
661  map<int,int> orientMap;
662  Array<OneD, unsigned int> map1(order_e), map2(order_e);
663  Array<OneD, int> sign1(order_e), sign2(order_e);
664 
665  locExpansion->GetFaceToElementMap(j, StdRegions::eDir1FwdDir1_Dir2FwdDir2, map1, sign1);
666  locExpansion->GetFaceToElementMap(j, locExpansion->GetForient(j), map2, sign2);
667 
668  for (k = 0; k < map1.num_elements(); ++k)
669  {
670  // Find the elemental co-efficient in the original
671  // mapping.
672  int idx = -1;
673  for (int l = 0; l < map2.num_elements(); ++l)
674  {
675  if (map1[k] == map2[l])
676  {
677  idx = l;
678  break;
679  }
680  }
681 
682  ASSERTL2(idx != -1, "Problem with face to element map!");
683  orientMap[k] = idx;
684  }
685 
686  for(k = 0; k < order_e; ++k)
687  {
688  vGlobalId = m_localToGlobalBndMap[k+cnt];
689  m_globalToUniversalBndMap[vGlobalId]
690  = id * maxDof + orientMap[k] + 1;
691  }
692  cnt += order_e;
693  }
694  }
695  }
696 
697  // Initialise GSlib and populate the unique map.
698  Array<OneD, long> tmp(m_globalToUniversalBndMap.num_elements());
699  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
700  {
701  tmp[i] = m_globalToUniversalBndMap[i];
702  }
703  m_bndGsh = m_gsh = Gs::Init(tmp, m_comm);
704  Gs::Unique(tmp, m_comm);
705  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
706  {
707  m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
708  }
709  }
710 
712  const ExpList &locExp,
713  const ExpListSharedPtr trace,
714  const PeriodicMap &perMap)
715  {
716  Array<OneD, int> tmp;
718  int i;
719  int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
720 
721  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
722 
723  int nTracePhys = trace->GetTotPoints();
724 
725  // Initialise the trace to universal maps.
727  Nektar::Array<OneD, int>(nTracePhys, -1);
729  Nektar::Array<OneD, int>(nTracePhys, -1);
730 
731  // Assume that each element of the expansion is of the same
732  // dimension.
733  nDim = locExpVector[0]->GetShapeDimension();
734 
735  if (nDim == 1)
736  {
737  maxQuad = (1 > maxQuad ? 1 : maxQuad);
738  }
739  else
740  {
741  for (i = 0; i < trace->GetExpSize(); ++i)
742  {
743  quad = trace->GetExp(i)->GetTotPoints();
744  if (quad > maxQuad)
745  {
746  maxQuad = quad;
747  }
748  }
749  }
750  m_comm->AllReduce(maxQuad, LibUtilities::ReduceMax);
751 
752  if (nDim == 1)
753  {
754  for (int i = 0; i < trace->GetExpSize(); ++i)
755  {
756  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
757  offset = trace->GetPhys_Offset(i);
758 
759  // Check to see if this vert is periodic. If it is, then we
760  // need use the unique eid of the two points
761  PeriodicMap::const_iterator it = perMap.find(eid);
762  if (perMap.count(eid) > 0)
763  {
764  PeriodicEntity ent = it->second[0];
765  if (ent.isLocal == false) // Not sure if true in 1D
766  {
767  eid = min(eid, ent.id);
768  }
769  }
770 
771  m_traceToUniversalMap[offset] = eid*maxQuad+1;
772  }
773  }
774  else
775  {
776  for (int i = 0; i < trace->GetExpSize(); ++i)
777  {
778  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
779  offset = trace->GetPhys_Offset(i);
780  quad = trace->GetExp(i)->GetTotPoints();
781 
782  // Check to see if this edge is periodic. If it is, then we
783  // need to reverse the trace order of one edge only in the
784  // universal map so that the data are reversed w.r.t each
785  // other. We do this by using the minimum of the two IDs.
786  PeriodicMap::const_iterator it = perMap.find(eid);
787  bool realign = false;
788  if (perMap.count(eid) > 0)
789  {
790  PeriodicEntity ent = it->second[0];
791  if (ent.isLocal == false)
792  {
793  realign = eid == min(eid, ent.id);
794  eid = min(eid, ent.id);
795  }
796  }
797 
798  for (int j = 0; j < quad; ++j)
799  {
800  m_traceToUniversalMap[j+offset] = eid*maxQuad+j+1;
801  }
802 
803  if (realign)
804  {
805  if (nDim == 2)
806  {
808  tmp = m_traceToUniversalMap+offset,
809  it->second[0].orient, quad);
810  }
811  else
812  {
814  tmp = m_traceToUniversalMap+offset,
815  it->second[0].orient,
816  trace->GetExp(i)->GetNumPoints(0),
817  trace->GetExp(i)->GetNumPoints(1));
818  }
819  }
820  }
821  }
822 
823  Array<OneD, long> tmp2(nTracePhys);
824  for (int i = 0; i < nTracePhys; ++i)
825  {
826  tmp2[i] = m_traceToUniversalMap[i];
827  }
828  m_traceGsh = Gs::Init(tmp2, m_comm);
829  Gs::Unique(tmp2, m_comm);
830  for (int i = 0; i < nTracePhys; ++i)
831  {
832  m_traceToUniversalMapUnique[i] = tmp2[i];
833  }
834  }
835 
837  Array<OneD, int> &toAlign,
839  int nquad1,
840  int nquad2)
841  {
842  if (orient == StdRegions::eBackwards)
843  {
844  ASSERTL1(nquad2 == 0, "nquad2 != 0 for reorienation");
845  for (int i = 0; i < nquad1/2; ++i)
846  {
847  swap(toAlign[i], toAlign[nquad1-1-i]);
848  }
849  }
850  else if (orient != StdRegions::eForwards)
851  {
852  ASSERTL1(nquad2 != 0, "nquad2 == 0 for reorienation");
853 
854  Array<OneD, int> tmp(nquad1*nquad2);
855 
856  // Copy transpose.
857  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
861  {
862  for (int i = 0; i < nquad2; ++i)
863  {
864  for (int j = 0; j < nquad1; ++j)
865  {
866  tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
867  }
868  }
869  }
870  else
871  {
872  for (int i = 0; i < nquad2; ++i)
873  {
874  for (int j = 0; j < nquad1; ++j)
875  {
876  tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
877  }
878  }
879  }
880 
881  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
885  {
886  // Reverse x direction
887  for (int i = 0; i < nquad2; ++i)
888  {
889  for (int j = 0; j < nquad1/2; ++j)
890  {
891  swap(tmp[i*nquad1 + j],
892  tmp[i*nquad1 + nquad1-j-1]);
893  }
894  }
895  }
896 
897  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
901  {
902  // Reverse y direction
903  for (int j = 0; j < nquad1; ++j)
904  {
905  for (int i = 0; i < nquad2/2; ++i)
906  {
907  swap(tmp[i*nquad1 + j],
908  tmp[(nquad2-i-1)*nquad1 + j]);
909  }
910  }
911  }
912  Vmath::Vcopy(nquad1*nquad2, tmp, 1, toAlign, 1);
913  }
914  }
915 
917  Array<OneD, NekDouble> &pGlobal) const
918  {
919  Gs::Gather(pGlobal, Gs::gs_add, m_traceGsh);
920  }
921 
923  {
924  return m_localToGlobalBndMap[i];
925  }
926 
928  {
929  return m_globalToUniversalBndMap[i];
930  }
931 
933  {
935  }
936 
938  {
939  return m_localToGlobalBndMap;
940  }
941 
943  {
945  }
946 
948  {
950  }
951 
953  const int i) const
954  {
955  return GetLocalToGlobalBndSign(i);
956  }
957 
959  const Array<OneD, const NekDouble>& loc,
960  Array<OneD, NekDouble>& global,
961  bool useComm ) const
962  {
963  AssembleBnd(loc,global);
964  }
965 
967  const NekVector<NekDouble>& loc,
968  NekVector< NekDouble>& global,
969  bool useComm) const
970  {
971  AssembleBnd(loc,global);
972  }
973 
975  const Array<OneD, const NekDouble>& global,
976  Array<OneD, NekDouble>& loc) const
977  {
978  GlobalToLocalBnd(global,loc);
979  }
980 
982  const NekVector<NekDouble>& global,
983  NekVector< NekDouble>& loc) const
984  {
985  GlobalToLocalBnd(global,loc);
986  }
987 
989  const Array<OneD, const NekDouble> &loc,
990  Array<OneD, NekDouble> &global) const
991  {
992  AssembleBnd(loc,global);
993  }
994 
996  const NekVector<NekDouble>& loc,
997  NekVector< NekDouble>& global) const
998  {
999  AssembleBnd(loc,global);
1000  }
1001 
1003  Array<OneD, NekDouble>& pGlobal) const
1004  {
1005  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
1006  }
1007 
1009  NekVector< NekDouble>& pGlobal) const
1010  {
1011  UniversalAssemble(pGlobal.GetPtr());
1012  }
1013 
1015  {
1016  return GetBndSystemBandWidth();
1017  }
1018 
1020  {
1021  return m_traceToUniversalMap[i];
1022  }
1023 
1025  {
1026  return m_traceToUniversalMapUnique[i];
1027  }
1028 
1030  {
1031  return m_numDirichletBndPhys;
1032  }
1033 
1036  {
1037  ASSERTL1(i >= 0 && i < m_elmtToTrace.num_elements(),
1038  "i is out of range");
1039  return m_elmtToTrace[i];
1040  }
1041 
1044  {
1045  return m_elmtToTrace;
1046  }
1047  } //namespace
1048 } // namespace
AssemblyMapDG()
Default constructor.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
bool m_signChange
Flag indicating if modes require sign reversal.
Definition: AssemblyMap.h:347
int m_numGlobalBndCoeffs
Total number of global boundary coefficients.
Definition: AssemblyMap.h:316
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:308
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:239
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:358
bool isLocal
Flag specifying if this entity is local to this partition.
int GetBndSystemBandWidth() const
Returns the bandwidth of the boundary system.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2067
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:333
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
AssemblyMapSharedPtr m_nextLevelLocalToGlobalMap
Map from the patches of the previous level to the patches of the current level.
Definition: AssemblyMap.h:397
Base class for constructing local to global mapping of degrees of freedom.
Definition: AssemblyMap.h:59
static gs_data * Init(const Nektar::Array< OneD, long > pId, const LibUtilities::CommSharedPtr &pComm, bool verbose=true)
Initialise Gather-Scatter map.
Definition: GsLib.hpp:166
void UniversalTraceAssemble(Array< OneD, NekDouble > &pGlobal) const
size_t m_hash
Hash for map.
Definition: AssemblyMap.h:311
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:270
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:390
GlobalSysSolnType m_solnType
The solution type of the global system.
Definition: AssemblyMap.h:365
int m_numGlobalDirBndCoeffs
Number of Global Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:320
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:200
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:392
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
int m_lowestStaticCondLevel
Lowest static condensation level.
Definition: AssemblyMap.h:399
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:350
virtual ~AssemblyMapDG()
Destructor.
int m_numLocalDirBndCoeffs
Number of Local Dirichlet Boundary Coefficients.
Definition: AssemblyMap.h:318
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:354
int m_numLocalBndCoeffs
Number of local boundary coefficients.
Definition: AssemblyMap.h:314
int m_staticCondLevel
The level of recursion in the case of multi-level static condensation.
Definition: AssemblyMap.h:386
Array< OneD, NekDouble > m_localToGlobalBndSign
Integer sign of local boundary coeffs to global space.
Definition: AssemblyMap.h:352
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:250
virtual void v_GlobalToLocal(const Array< OneD, const NekDouble > &global, Array< OneD, NekDouble > &loc) const
virtual void v_LocalToGlobal(const Array< OneD, const NekDouble > &loc, Array< OneD, NekDouble > &global, bool useComm) const
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:360
Array< OneD, int > m_globalToUniversalBndMapUnique
Integer map of unique process coeffs to universal space (signed)
Definition: AssemblyMap.h:362
int m_numGlobalCoeffs
Total number of global coefficients.
Definition: AssemblyMap.h:344
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:228
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:1061
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:388
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