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 = 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  // Order list according to m_offset_elmt_id details in Exp
589  // so that triangules are listed first and then quads
590  eid = locExp.GetOffset_Elmt_Id(i);
591  locExpansion = locExpVector[eid];
592  nDim = locExpansion->GetShapeDimension();
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 = locExpansion->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  locExpansion->GetEdgeToElementMap(j, StdRegions::eForwards, map1, sign1);
624  locExpansion->GetEdgeToElementMap(j, locExpansion->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 = locExpansion->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  locExpansion->GetFaceToElementMap(j, StdRegions::eDir1FwdDir1_Dir2FwdDir2, map1, sign1);
669  locExpansion->GetFaceToElementMap(j, locExpansion->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,
964  bool useComm ) const
965  {
966  AssembleBnd(loc,global);
967  }
968 
970  const NekVector<NekDouble>& loc,
971  NekVector< NekDouble>& global,
972  bool useComm) const
973  {
974  AssembleBnd(loc,global);
975  }
976 
978  const Array<OneD, const NekDouble>& global,
979  Array<OneD, NekDouble>& loc) const
980  {
981  GlobalToLocalBnd(global,loc);
982  }
983 
985  const NekVector<NekDouble>& global,
986  NekVector< NekDouble>& loc) const
987  {
988  GlobalToLocalBnd(global,loc);
989  }
990 
992  const Array<OneD, const NekDouble> &loc,
993  Array<OneD, NekDouble> &global) const
994  {
995  AssembleBnd(loc,global);
996  }
997 
999  const NekVector<NekDouble>& loc,
1000  NekVector< NekDouble>& global) const
1001  {
1002  AssembleBnd(loc,global);
1003  }
1004 
1006  Array<OneD, NekDouble>& pGlobal) const
1007  {
1008  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
1009  }
1010 
1012  NekVector< NekDouble>& pGlobal) const
1013  {
1014  UniversalAssemble(pGlobal.GetPtr());
1015  }
1016 
1018  {
1019  return GetBndSystemBandWidth();
1020  }
1021 
1023  {
1024  return m_traceToUniversalMap[i];
1025  }
1026 
1028  {
1029  return m_traceToUniversalMapUnique[i];
1030  }
1031 
1033  {
1034  return m_numDirichletBndPhys;
1035  }
1036 
1039  {
1040  ASSERTL1(i >= 0 && i < m_elmtToTrace.num_elements(),
1041  "i is out of range");
1042  return m_elmtToTrace[i];
1043  }
1044 
1047  {
1048  return m_elmtToTrace;
1049  }
1050  } //namespace
1051 } // 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:2075
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
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:2100
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