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 
36 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h>
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>
46 #include <LocalRegions/Expansion2D.h>
47 #include <LocalRegions/Expansion3D.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  {
60  AssemblyMapDG::AssemblyMapDG():
61  m_numDirichletBndPhys(0)
62  {
63  }
64 
65  AssemblyMapDG::~AssemblyMapDG()
66  {
67  }
68 
69 
70  /**
71  *
72  */
73  AssemblyMapDG::AssemblyMapDG(
74  const LibUtilities::SessionReaderSharedPtr &pSession,
75  const SpatialDomains::MeshGraphSharedPtr &graph1D,
76  const ExpList0DSharedPtr &trace,
77  const ExpList &locExp,
78  const Array<OneD, const MultiRegions::ExpListSharedPtr> &bndCondExp,
79  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &bndCond,
80  const PeriodicMap &periodicVerts,
81  const std::string variable)
82  : AssemblyMap(pSession,variable)
83  {
84  int i,j;
85  int cnt, vid, gid;
86  int ntrace_exp = trace->GetExpSize();
87  int nbnd = bndCondExp.num_elements();
88 
89  int nel = locExp.GetExp()->size();
90  map<int, int> MeshVertId;
91 
92  // determine mapping from geometry edges to trace
93  for(i = 0; i < ntrace_exp; ++i)
94  {
95  int id = trace->GetExp(i)->GetGeom()->GetVid(0);
96  MeshVertId[id] = i;
97  }
98 
99 
100  Array<OneD, StdRegions::StdExpansionSharedPtr> vertmap(2*nel);
101  m_elmtToTrace = Array<OneD, Array<OneD,StdRegions::StdExpansionSharedPtr> >(nel);
102 
103  // set up vert expansions links;
104  cnt = 0;
105  for(i = 0; i < nel; ++i)
106  {
107  m_elmtToTrace[i] = vertmap + cnt;
108 
109  for(j = 0; j < locExp.GetExp(i)->GetNverts(); ++j)
110  {
111  int id = ((locExp.GetExp(i)->GetGeom())->GetVertex(j))->GetVid();
112 
113  if(MeshVertId.count(id) > 0)
114  {
115  m_elmtToTrace[i][j] = (*trace).GetExp(MeshVertId.find(id)->second)->as<LocalRegions::PointExp>();
116 
117  }
118  else
119  {
120  ASSERTL0(false,"Failed to find edge map");
121  }
122  }
123  cnt += 2;
124  }
125 
126  // set up Local to Continuous mapping
127  Array<OneD,unsigned int> vmap;
128  LocalRegions::SegExpSharedPtr locSegExp;
129 
130  const boost::shared_ptr<LocalRegions::ExpansionVector> exp1D = locExp.GetExp();
131 
132  m_numGlobalBndCoeffs = exp1D->size()+1;
133  m_numGlobalCoeffs = m_numGlobalBndCoeffs;
134  m_numLocalBndCoeffs = 2*exp1D->size();
135  m_numLocalCoeffs = m_numLocalBndCoeffs;
136  m_localToGlobalBndMap = Array<OneD, int>(m_numLocalBndCoeffs,-1);
137  m_localToGlobalBndSign = Array<OneD, NekDouble>(m_numLocalBndCoeffs,1.0);
138  m_signChange = true;
139  m_staticCondLevel = 0;
140  m_numPatches = exp1D->size();
141  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
142  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(m_numPatches);
143  for(i = 0; i < m_numPatches; ++i)
144  {
145  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) (*exp1D)[i]->NumDGBndryCoeffs();
146  m_numLocalIntCoeffsPerPatch[i] = (unsigned int) 0;
147  }
148 
149  map<int, int> MeshVertToLocalVert;
150 
151  // Order the Dirichlet vertices first.
152  gid = 0;
153  for(i = 0; i < nbnd; i++)
154  {
155  if(bndCond[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
156  {
157  m_numDirichletBndPhys++;
158  vid = bndCondExp[i]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
159 
160  MeshVertToLocalVert[vid] = gid++;
161  }
162  }
163 
164  // set up simple map based on vertex and edge id's
165  cnt = 0;
166  for(i = 0; i < exp1D->size(); ++i)
167  {
168  if((locSegExp = (*exp1D)[i]->as<LocalRegions::SegExp>()))
169  {
170  locSegExp->GetBoundaryMap(vmap);
171 
172  for(j = 0; j < locSegExp->GetNverts(); ++j)
173  {
174  vid = (locSegExp->GetGeom1D())->GetVid(j);
175 
176  if(MeshVertToLocalVert.count(vid) == 0)
177  {
178  MeshVertToLocalVert[vid] = gid++;
179  }
180 
181  m_localToGlobalBndMap[cnt+j] =
182  MeshVertToLocalVert.find(vid)->second;
183  }
184  cnt += locSegExp->NumBndryCoeffs();
185  }
186  else
187  {
188  ASSERTL0(false,"dynamic cast to a segment expansion failed");
189  }
190  }
191 
192  // Set up boundary mapping
193  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD, int >(nbnd);
194  m_bndCondCoeffsToGlobalCoeffsSign = Array<OneD, NekDouble >(nbnd,1.0);
195  m_numLocalDirBndCoeffs = 0;
196  m_numDirichletBndPhys = 0;
197 
198  for(i = 0; i < nbnd; ++i)
199  {
200  vid = bndCondExp[i]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
201  m_bndCondCoeffsToGlobalCoeffsMap[i] = MeshVertToLocalVert.find(vid)->second;
202 
203  if(bndCond[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
204  {
205  m_numLocalDirBndCoeffs += 1;
206  m_numDirichletBndPhys += 1;
207  }
208  }
209 
210  m_numGlobalDirBndCoeffs = m_numLocalDirBndCoeffs;
211  CalculateBndSystemBandWidth();
212 
213 
214  m_bndCondTraceToGlobalTraceMap = Array<OneD, int>(nbnd);
215  for(i = 0; i < bndCondExp.num_elements(); ++i)
216  {
217  int id = bndCondExp[i]->GetExp(0)->GetGeom()->GetVid(0);
218  m_bndCondTraceToGlobalTraceMap[i] =
219  MeshVertId.find(id)->second;
220  }
221 
222  // Now set up mapping from global coefficients to universal.
223 #if 1 // Routines need debugging -> Currently causing crash when turned on
224  ExpListSharedPtr tr = boost::dynamic_pointer_cast<ExpList>(trace);
225  SetUpUniversalDGMap (locExp);
226  SetUpUniversalTraceMap(locExp, tr, periodicVerts);
227 #endif
228  m_hash = boost::hash_range(m_localToGlobalBndMap.begin(),
229  m_localToGlobalBndMap.end());
230 
231 
232  }
233 
234 
235  /**
236  *
237  */
238  AssemblyMapDG::AssemblyMapDG(const LibUtilities::SessionReaderSharedPtr &pSession,
239  const SpatialDomains::MeshGraphSharedPtr &graph2D,
240  const ExpList1DSharedPtr &trace,
241  const ExpList &locExp,
242  const Array<OneD, MultiRegions::ExpListSharedPtr> &bndCondExp,
243  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> &bndCond,
244  const PeriodicMap &periodicEdges,
245  const std::string variable) :
246  AssemblyMap(pSession,variable)
247  {
248 
249  int i,j,k,cnt,eid, id, id1, order_e,gid;
250  int ntrace_exp = trace->GetExpSize();
251  int nbnd = bndCondExp.num_elements();
252  LocalRegions::SegExpSharedPtr locSegExp,locSegExp1;
253  LocalRegions::QuadExpSharedPtr locQuadExp;
254  LocalRegions::TriExpSharedPtr locTriExp;
255  SpatialDomains::Geometry1DSharedPtr SegGeom;
256 
257  const boost::shared_ptr<LocalRegions::ExpansionVector> exp2D = locExp.GetExp();
258  int nel = exp2D->size();
259 
260  map<int, int> MeshEdgeId;
261 
262  m_signChange = true;
263 
264  // determine mapping from geometry edges to trace
265  for(i = 0; i < ntrace_exp; ++i)
266  {
267  if((locSegExp = trace->GetExp(i)->as<LocalRegions::SegExp>()))
268  {
269  id = (locSegExp->GetGeom1D())->GetEid();
270  MeshEdgeId[id] = i;
271  }
272  else
273  {
274  ASSERTL0(false,"Dynamics cast to segment expansion failed");
275  }
276  }
277 
278  // Count total number of edges
279  cnt = 0;
280  for(i = 0; i < nel; ++i)
281  {
282  cnt += (*exp2D)[i]->GetNedges();
283  }
284 
285  Array<OneD, StdRegions::StdExpansionSharedPtr> edgemap(cnt);
286  m_elmtToTrace = Array<OneD, Array<OneD,StdRegions::StdExpansionSharedPtr> >(nel);
287 
288  // set up edge expansions links;
289  cnt = 0;
290  for(i = 0; i < nel; ++i)
291  {
292  m_elmtToTrace[i] = edgemap + cnt;
293 
294  if((locQuadExp = (*exp2D)[i]->as<LocalRegions::QuadExp>()))
295  {
296  for(j = 0; j < locQuadExp->GetNedges(); ++j)
297  {
298  SegGeom = (locQuadExp->GetGeom2D())->GetEdge(j);
299 
300  id = SegGeom->GetEid();
301 
302  if(MeshEdgeId.count(id) > 0)
303  {
304  m_elmtToTrace[i][j] = (*trace).GetExp(MeshEdgeId
305  .find(id)->second)->as<LocalRegions::SegExp>();
306 
307  }
308  else
309  {
310  ASSERTL0(false,"Failed to find edge map");
311  }
312  }
313  }
314  else if((locTriExp = (*exp2D)[i]->as<LocalRegions::TriExp>()))
315  {
316  for(j = 0; j < locTriExp->GetNedges(); ++j)
317  {
318  SegGeom = (locTriExp->GetGeom2D())->GetEdge(j);
319 
320  id = SegGeom->GetEid();
321 
322  if(MeshEdgeId.count(id) > 0)
323  {
324  m_elmtToTrace[i][j] = (*trace).GetExp((MeshEdgeId
325  .find(id))->second)->as<LocalRegions::SegExp>();
326 
327  }
328  else
329  {
330  ASSERTL0(false,"Failed to find edge map");
331  }
332  }
333 
334  }
335  else
336  {
337  ASSERTL0(false,"dynamic cast to a local 2D expansion failed");
338  }
339  cnt += (*exp2D)[i]->GetNedges();
340  }
341 
342  // Set up boundary mapping
343  cnt = 0;
344  for(i = 0; i < nbnd; ++i)
345  {
346  cnt += bndCondExp[i]->GetExpSize();
347  }
348 
349  m_numLocalDirBndCoeffs = 0;
350  m_numDirichletBndPhys = 0;
351 
352  cnt = 0;
353  for(i = 0; i < bndCondExp.num_elements(); ++i)
354  {
355  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
356  {
357  if((locSegExp = bndCondExp[i]->GetExp(j)
358  ->as<LocalRegions::SegExp>()))
359  {
360  SegGeom = locSegExp->GetGeom1D();
361  id = SegGeom->GetEid();
362  }
363  else
364  {
365  ASSERTL0(false,"dynamic cast to a local Segment expansion failed");
366  }
367 
368  if(bndCond[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
369  {
370  m_numLocalDirBndCoeffs += locSegExp->GetNcoeffs();
371  m_numDirichletBndPhys += locSegExp->GetTotPoints();
372  }
373 
374  }
375  cnt += j;
376  }
377 
378  // Set up integer mapping array and sign change for each
379  // degree of freedom + initialise some more data members
380  m_staticCondLevel = 0;
381  m_lowestStaticCondLevel = 0;
382  m_numPatches = nel;
383  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
384  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
385  int nbndry = 0;
386  for(i = 0; i < nel; ++i) // count number of elements in array
387  {
388  eid = locExp.GetOffset_Elmt_Id(i);
389  nbndry += (*exp2D)[eid]->NumDGBndryCoeffs();
390  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) (*exp2D)[eid]->NumDGBndryCoeffs();
391  m_numLocalIntCoeffsPerPatch[i] = (unsigned int) 0;
392  }
393 
394  m_numGlobalDirBndCoeffs = m_numLocalDirBndCoeffs;
395 
396  m_numLocalBndCoeffs = nbndry;
397  m_numLocalCoeffs = nbndry;
398  m_localToGlobalBndMap = Array<OneD, int > (nbndry);
399  m_localToGlobalBndSign = Array<OneD, NekDouble > (nbndry,1);
400 
401  // Set up array for potential mesh optimsation
402  Array<OneD,int> TraceElmtGid(ntrace_exp,-1);
403  int nDir = 0;
404  cnt = 0;
405 
406  // We are now going to construct a graph of the mesh
407  // which can be reordered depending on the type of solver we would
408  // like to use.
409  typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
410 
411  BoostGraph boostGraphObj;
412  int trace_id,trace_id1;
413 
414  // make trace edge renumbering map where first solved
415  // edge starts at 0 so we can set up graph.
416  for(i = 0; i < ntrace_exp; ++i)
417  {
418  if(trace->GetCoeff_Offset(i) >= m_numLocalDirBndCoeffs)
419  {
420  // Initial put in element ordering (starting
421  // from zero) into TraceElmtGid
422  boost::add_vertex(boostGraphObj);
423  TraceElmtGid[i] = cnt++;
424  }
425  else
426  {
427  // Use existing offset for Dirichlet edges
428  TraceElmtGid[i] = trace->GetCoeff_Offset(i);
429  nDir++;
430  }
431  }
432 
433  // Set up boost Graph
434  for(i = 0; i < nel; ++i)
435  {
436  eid = locExp.GetOffset_Elmt_Id(i);
437  nbndry += (*exp2D)[eid]->NumDGBndryCoeffs();
438 
439  for(j = 0; j < (*exp2D)[eid]->GetNedges(); ++j)
440  {
441  locSegExp = m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
442  SegGeom = locSegExp->GetGeom1D();
443 
444  // Add edge to boost graph for non-Dirichlet Boundary
445  id = SegGeom->GetEid();
446  trace_id = MeshEdgeId.find(id)->second;
447  if(trace->GetCoeff_Offset(trace_id) >= m_numLocalDirBndCoeffs)
448  {
449  for(k = j+1; k < (*exp2D)[eid]->GetNedges(); ++k)
450  {
451  locSegExp1 = m_elmtToTrace[eid][k]->as<LocalRegions::SegExp>();
452  SegGeom = locSegExp1->GetGeom1D();
453 
454  id1 = SegGeom->GetEid();
455  trace_id1 = MeshEdgeId.find(id1)->second;
456  if(trace->GetCoeff_Offset(trace_id1)
457  >= m_numLocalDirBndCoeffs)
458  {
459  boost::add_edge( (size_t) TraceElmtGid[trace_id], (size_t) TraceElmtGid[trace_id1], boostGraphObj);
460  }
461  }
462  }
463  }
464  }
465 
466 
467  int nGraphVerts = ntrace_exp-nDir;
468  Array<OneD, int> perm(nGraphVerts);
469  Array<OneD, int> iperm(nGraphVerts);
470  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
471  Array<OneD, int> vwgts(nGraphVerts);
472  for(i = 0; i < nGraphVerts; ++i)
473  {
474  vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
475  }
476 
477  if(nGraphVerts)
478  {
479  switch(m_solnType)
480  {
481  case eDirectFullMatrix:
482  case eIterativeFull:
483  case ePETScFullMatrix:
484  case ePETScStaticCond:
485  case eIterativeStaticCond:
486  {
487  NoReordering(boostGraphObj,perm,iperm);
488  break;
489  }
490  case eDirectStaticCond:
491  {
492  CuthillMckeeReordering(boostGraphObj,perm,iperm);
493  break;
494  }
495  case eDirectMultiLevelStaticCond:
496  case eIterativeMultiLevelStaticCond:
497  {
498  MultiLevelBisectionReordering(boostGraphObj,perm,iperm,bottomUpGraph);
499  break;
500  }
501  default:
502  {
503  ASSERTL0(false,"Unrecognised solution type");
504  }
505  }
506  }
507 
508  // Recast the permutation so that it can be
509  // used as a map Form old trace edge ID to new trace
510  // edge ID
511  cnt = m_numLocalDirBndCoeffs;
512  for(i = 0; i < ntrace_exp-nDir; ++i)
513  {
514  TraceElmtGid[perm[i]+nDir]=cnt;
515  cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
516  }
517 
518  // Now have trace edges Gid position
519  nbndry = cnt = 0;
520  for(i = 0; i < nel; ++i)
521  {
522  // order list according to m_offset_elmt_id details in
523  // Exp2D so that triangules are listed first and then
524  // quads
525  eid = locExp.GetOffset_Elmt_Id(i);
526  nbndry += (*exp2D)[eid]->NumDGBndryCoeffs();
527 
528  for(j = 0; j < (*exp2D)[eid]->GetNedges(); ++j)
529  {
530  locSegExp = m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
531  SegGeom = locSegExp->GetGeom1D();
532 
533  id = SegGeom->GetEid();
534  gid = TraceElmtGid[MeshEdgeId.find(id)->second];
535 
536  //Peter order_e = locSegExp->GetNcoeffs();
537  order_e = (*exp2D)[eid]->GetEdgeNcoeffs(j);
538 
539  if((*exp2D)[eid]->GetEorient(j) == StdRegions::eForwards)
540  {
541  for(k = 0; k < order_e; ++k)
542  {
543  m_localToGlobalBndMap[k+cnt] = gid + k;
544  }
545  }
546  else // backwards orientated
547  {
548  switch(locSegExp->GetBasisType(0))
549  {
550  case LibUtilities::eModified_A:
551  {
552  // reverse vertex order
553  m_localToGlobalBndMap[cnt] = gid + 1;
554  m_localToGlobalBndMap[cnt+1] = gid;
555  for (k = 2; k < order_e; ++k)
556  {
557  m_localToGlobalBndMap[k+cnt] = gid + k;
558  }
559 
560  // negate odd modes
561  for(k = 3; k < order_e; k+=2)
562  {
563  m_localToGlobalBndSign[cnt+k] = -1.0;
564  }
565  break;
566  }
567  case LibUtilities::eGLL_Lagrange:
568  {
569  // reverse order
570  for(k = 0; k < order_e; ++k)
571  {
572  m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
573  }
574  break;
575  }
576  case LibUtilities::eGauss_Lagrange:
577  {
578  // reverse order
579  for(k = 0; k < order_e; ++k)
580  {
581  m_localToGlobalBndMap[cnt+order_e-k-1] = gid + k;
582  }
583  break;
584  }
585  default:
586  {
587  ASSERTL0(false,"Boundary type not permitted");
588  }
589  }
590  }
591 
592  cnt += order_e;
593  }
594  }
595 
596  // set up m_bndCondCoeffsToGlobalCoeffsMap to align with map
597  cnt = 0;
598  for(i = 0; i < nbnd; ++i)
599  {
600  cnt += bndCondExp[i]->GetNcoeffs();
601  }
602 
603  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD,int >(cnt);
604 
605  // Number of boundary expansions
606  int nbndexp = 0;
607  for(cnt = i = 0; i < nbnd; ++i)
608  {
609  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
610  {
611  if((locSegExp = bndCondExp[i]->GetExp(j)
612  ->as<LocalRegions::SegExp>()))
613  {
614  nbndexp++;
615  SegGeom = locSegExp->GetGeom1D();
616  id = SegGeom->GetEid();
617  gid = TraceElmtGid[MeshEdgeId.find(id)->second];
618 
619  order_e = locSegExp->GetNcoeffs();
620 
621  // Since boundary information is defined to be
622  // aligned with the geometry just use forward
623  // defintiion for gid's
624  for(k = 0; k < order_e; ++k)
625  {
626  m_bndCondCoeffsToGlobalCoeffsMap[cnt++] = gid + k;
627  }
628  }
629  }
630  }
631 
632  m_numGlobalBndCoeffs = trace->GetNcoeffs();
633  m_numGlobalCoeffs = m_numGlobalBndCoeffs;
634 
635  CalculateBndSystemBandWidth();
636 
637  if( m_solnType == eDirectMultiLevelStaticCond && nGraphVerts )
638  {
639  if(m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
640  {
641 
642  Array<OneD, int> vwgts_perm(nGraphVerts);
643 
644  for(int i = 0; i < nGraphVerts; i++)
645  {
646  vwgts_perm[i] = vwgts[perm[i]];
647  }
648 
649  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
650 
651  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
652  AllocateSharedPtr(this,bottomUpGraph);
653  }
654  }
655 
656  cnt = 0;
657  m_bndCondTraceToGlobalTraceMap = Array<OneD, int >(nbndexp);
658  for(i = 0; i < bndCondExp.num_elements(); ++i)
659  {
660  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
661  {
662  if((locSegExp = bndCondExp[i]->GetExp(j)
663  ->as<LocalRegions::SegExp>()))
664  {
665  SegGeom = locSegExp->GetGeom1D();
666  id = SegGeom->GetEid();
667 
668  m_bndCondTraceToGlobalTraceMap[cnt++] = MeshEdgeId.find(id)->second;
669  }
670  }
671  }
672 
673  // Now set up mapping from global coefficients to universal.
674  ExpListSharedPtr tr = boost::dynamic_pointer_cast<ExpList>(trace);
675  SetUpUniversalDGMap (locExp);
676  SetUpUniversalTraceMap(locExp, tr, periodicEdges);
677 
678  m_hash = boost::hash_range(m_localToGlobalBndMap.begin(),
679  m_localToGlobalBndMap.end());
680  }
681 
682  /**
683  * Constructor for trace map for three-dimensional expansion.
684  */
685  AssemblyMapDG::AssemblyMapDG(
686  const LibUtilities::SessionReaderSharedPtr &pSession,
687  const SpatialDomains::MeshGraphSharedPtr &graph3D,
688  const ExpList2DSharedPtr &trace,
689  const ExpList &locExp,
690  const Array<OneD, MultiRegions::ExpListSharedPtr> &bndCondExp,
691  const Array<OneD, SpatialDomains::BoundaryConditionShPtr> &bndCond,
692  const PeriodicMap &periodicFaces,
693  const std::string variable):
694  AssemblyMap(pSession,variable)
695  {
696  int i,j,k,cnt,eid, id, id1, order_e,gid;
697  int ntrace_exp = trace->GetExpSize();
698  int nbnd = bndCondExp.num_elements();
699  LocalRegions::QuadExpSharedPtr locQuadExp, locQuadExp1;
700  LocalRegions::TriExpSharedPtr locTriExp, locTriExp1;
701  LocalRegions::HexExpSharedPtr locHexExp;
702  LocalRegions::PrismExpSharedPtr locPrismExp;
703  LocalRegions::PyrExpSharedPtr locPyrExp;
704  LocalRegions::TetExpSharedPtr locTetExp;
705  SpatialDomains::Geometry2DSharedPtr FaceGeom;
706  StdRegions::StdExpansionSharedPtr locBndExp;
707 
708  const boost::shared_ptr<LocalRegions::ExpansionVector> exp3D =
709  locExp.GetExp();
710  int nel = exp3D->size();
711 
712  map<int, int> MeshFaceId;
713 
714  m_signChange = true;
715 
716  // determine mapping from geometry edges to trace
717  for(i = 0; i < ntrace_exp; ++i)
718  {
719  id = trace->GetExp(i)->as<LocalRegions::Expansion2D>()
720  ->GetGeom2D()->GetFid();
721  MeshFaceId[id] = i;
722  }
723 
724  // Count total number of faces
725  cnt = 0;
726  for(i = 0; i < nel; ++i)
727  {
728  cnt += (*exp3D)[i]->GetNfaces();
729  }
730 
731  Array<OneD, StdRegions::StdExpansionSharedPtr> facemap(cnt);
732  m_elmtToTrace = Array<OneD,
733  Array<OneD, StdRegions::StdExpansionSharedPtr> >(nel);
734 
735  // set up face expansions links;
736  cnt = 0;
737  for(i = 0; i < nel; ++i)
738  {
739  m_elmtToTrace[i] = facemap + cnt;
740 
741  for(j = 0; j < (*exp3D)[i]->GetNfaces(); ++j)
742  {
743  id = (*exp3D)[i]->as<LocalRegions::Expansion3D>()
744  ->GetGeom3D()->GetFid(j);
745 
746  if(MeshFaceId.count(id) > 0)
747  {
748  m_elmtToTrace[i][j] =
749  trace->GetExp(MeshFaceId.find(id)->second);
750  }
751  else
752  {
753  ASSERTL0(false,"Failed to find face map");
754  }
755  }
756 
757  cnt += (*exp3D)[i]->GetNfaces();
758  }
759 
760  // Set up boundary mapping
761  cnt = 0;
762  for(i = 0; i < nbnd; ++i)
763  {
764  cnt += bndCondExp[i]->GetExpSize();
765  }
766 
767  set<int> dirFaces;
768 
769  m_numLocalDirBndCoeffs = 0;
770  m_numDirichletBndPhys = 0;
771 
772  cnt = 0;
773  for(i = 0; i < bndCondExp.num_elements(); ++i)
774  {
775  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
776  {
777  locBndExp = bndCondExp[i]->GetExp(j);
778  FaceGeom = locBndExp->as<LocalRegions::Expansion2D>()
779  ->GetGeom2D();
780  id = FaceGeom->GetFid();
781 
782  if(bndCond[i]->GetBoundaryConditionType() ==
783  SpatialDomains::eDirichlet)
784  {
785  m_numLocalDirBndCoeffs += locBndExp->GetNcoeffs();
786  m_numDirichletBndPhys += locBndExp->GetTotPoints();
787  dirFaces.insert(id);
788  }
789  }
790 
791  cnt += j;
792  }
793 
794  // Set up integer mapping array and sign change for each
795  // degree of freedom + initialise some more data members
796  m_staticCondLevel = 0;
797  m_lowestStaticCondLevel = 0;
798  m_numPatches = nel;
799  m_numLocalBndCoeffsPerPatch = Array<OneD, unsigned int>(nel);
800  m_numLocalIntCoeffsPerPatch = Array<OneD, unsigned int>(nel);
801 
802  int nbndry = 0;
803  for(i = 0; i < nel; ++i) // count number of elements in array
804  {
805  eid = locExp.GetOffset_Elmt_Id(i);
806  nbndry += (*exp3D)[eid]->NumDGBndryCoeffs();
807  m_numLocalBndCoeffsPerPatch[i] = (unsigned int) (*exp3D)[eid]->NumDGBndryCoeffs();
808  m_numLocalIntCoeffsPerPatch[i] = (unsigned int) 0;
809  }
810 
811  m_numGlobalDirBndCoeffs = m_numLocalDirBndCoeffs;
812  m_numLocalBndCoeffs = nbndry;
813  m_numLocalCoeffs = nbndry;
814  m_localToGlobalBndMap = Array<OneD, int> (nbndry);
815  m_localToGlobalBndSign = Array<OneD, NekDouble> (nbndry,1);
816 
817  // Set up array for potential mesh optimsation
818  Array<OneD,int> FaceElmtGid(ntrace_exp,-1);
819  int nDir = 0;
820  cnt = 0;
821 
822  // We are now going to construct a graph of the mesh
823  // which can be reordered depending on the type of solver we would
824  // like to use.
825  typedef boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS> BoostGraph;
826 
827  BoostGraph boostGraphObj;
828  int face_id, face_id1;
829  int dirOffset = 0;
830 
831  // make trace face renumbering map where first solved
832  // face starts at 0 so we can set up graph.
833  for(i = 0; i < ntrace_exp; ++i)
834  {
835  id = trace->GetExp(i)->as<LocalRegions::Expansion2D>()
836  ->GetGeom2D()->GetFid();
837 
838  if (dirFaces.count(id) == 0)
839  {
840  // Initial put in element ordering (starting
841  // from zero) into FaceElmtGid
842  boost::add_vertex(boostGraphObj);
843  FaceElmtGid[i] = cnt++;
844  }
845  else
846  {
847  // Use existing offset for Dirichlet edges
848  FaceElmtGid[i] = dirOffset;
849  dirOffset += trace->GetExp(i)->GetNcoeffs();
850  nDir++;
851  }
852  }
853 
854  // Set up boost Graph
855  for(i = 0; i < nel; ++i)
856  {
857  eid = locExp.GetOffset_Elmt_Id(i);
858 
859  for(j = 0; j < (*exp3D)[eid]->GetNfaces(); ++j)
860  {
861  // Add face to boost graph for non-Dirichlet Boundary
862  FaceGeom = m_elmtToTrace[eid][j]
863  ->as<LocalRegions::Expansion2D>()->GetGeom2D();
864  id = FaceGeom->GetFid();
865  face_id = MeshFaceId.find(id)->second;
866 
867  if(dirFaces.count(id) == 0)
868  {
869  for(k = j+1; k < (*exp3D)[eid]->GetNfaces(); ++k)
870  {
871  FaceGeom = m_elmtToTrace[eid][k]
872  ->as<LocalRegions::Expansion2D>()->GetGeom2D();
873  id1 = FaceGeom->GetFid();
874  face_id1 = MeshFaceId.find(id1)->second;
875 
876  if(dirFaces.count(id1) == 0)
877  {
878  boost::add_edge((size_t) FaceElmtGid[face_id],
879  (size_t) FaceElmtGid[face_id1],
880  boostGraphObj);
881  }
882  }
883  }
884  }
885  }
886 
887  int nGraphVerts = ntrace_exp - nDir;
888  Array<OneD, int> perm (nGraphVerts);
889  Array<OneD, int> iperm(nGraphVerts);
890  Array<OneD, int> vwgts(nGraphVerts);
891  BottomUpSubStructuredGraphSharedPtr bottomUpGraph;
892 
893  for(i = 0; i < nGraphVerts; ++i)
894  {
895  vwgts[i] = trace->GetExp(i+nDir)->GetNcoeffs();
896  }
897 
898  if(nGraphVerts)
899  {
900  switch(m_solnType)
901  {
902  case eDirectFullMatrix:
903  case eIterativeFull:
904  case eIterativeStaticCond:
905  case ePETScFullMatrix:
906  case ePETScStaticCond:
907  {
908  NoReordering(boostGraphObj,perm,iperm);
909  break;
910  }
911  case eDirectStaticCond:
912  {
913  CuthillMckeeReordering(boostGraphObj,perm,iperm);
914  break;
915  }
916  case eDirectMultiLevelStaticCond:
917  case eIterativeMultiLevelStaticCond:
918  {
919  MultiLevelBisectionReordering(boostGraphObj,perm,iperm,
920  bottomUpGraph);
921  break;
922  }
923  default:
924  {
925  ASSERTL0(false,"Unrecognised solution type");
926  }
927  }
928  }
929 
930  // Recast the permutation so that it can be used as a map Form old
931  // trace face ID to new trace face ID
932  cnt = m_numLocalDirBndCoeffs;
933  for(i = 0; i < ntrace_exp - nDir; ++i)
934  {
935  FaceElmtGid[perm[i]+nDir] = cnt;
936  cnt += trace->GetExp(perm[i]+nDir)->GetNcoeffs();
937  }
938 
939  // Now have trace edges Gid position
940  cnt = 0;
941  for(i = 0; i < nel; ++i)
942  {
943  // order list according to m_offset_elmt_id details in Exp3D
944  eid = locExp.GetOffset_Elmt_Id(i);
945 
946  for(j = 0; j < (*exp3D)[eid]->GetNfaces(); ++j)
947  {
948  FaceGeom = m_elmtToTrace[eid][j]
949  ->as<LocalRegions::Expansion2D>()->GetGeom2D();
950  id = FaceGeom->GetFid();
951  gid = FaceElmtGid[MeshFaceId.find(id)->second];
952  order_e = (*exp3D)[eid]->GetFaceNcoeffs(j);
953 
954  std::map<int,int> orientMap;
955 
956  Array<OneD, unsigned int> elmMap1 (order_e);
957  Array<OneD, int> elmSign1(order_e);
958  Array<OneD, unsigned int> elmMap2 (order_e);
959  Array<OneD, int> elmSign2(order_e);
960 
961  StdRegions::Orientation fo = (*exp3D)[eid]->GetFaceOrient(j);
962 
963  // Construct mapping which will permute global IDs
964  // according to face orientations.
965  (*exp3D)[eid]->GetFaceToElementMap(j,fo,elmMap1,elmSign1);
966  (*exp3D)[eid]->GetFaceToElementMap(
967  j,StdRegions::eDir1FwdDir1_Dir2FwdDir2,elmMap2,elmSign2);
968 
969  for (k = 0; k < elmMap1.num_elements(); ++k)
970  {
971  // Find the elemental co-efficient in the original
972  // mapping.
973  int idx = -1;
974  for (int l = 0; l < elmMap2.num_elements(); ++l)
975  {
976  if (elmMap1[k] == elmMap2[l])
977  {
978  idx = l;
979  break;
980  }
981  }
982 
983  ASSERTL2(idx != -1, "Problem with face to element map!");
984  orientMap[k] = idx;
985  }
986 
987  for(k = 0; k < order_e; ++k)
988  {
989  m_localToGlobalBndMap [k+cnt] = gid + orientMap[k];
990  m_localToGlobalBndSign[k+cnt] = elmSign2[orientMap[k]];
991  }
992 
993  cnt += order_e;
994  }
995  }
996 
997  // set up m_bndCondCoeffsToGlobalCoeffsMap to align with map
998  cnt = 0;
999  for(i = 0; i < nbnd; ++i)
1000  {
1001  cnt += bndCondExp[i]->GetNcoeffs();
1002  }
1003 
1004  m_bndCondCoeffsToGlobalCoeffsMap = Array<OneD,int>(cnt);
1005 
1006  // Number of boundary expansions
1007  int nbndexp = 0, bndOffset, bndTotal = 0;
1008  for(cnt = i = 0; i < nbnd; ++i)
1009  {
1010  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
1011  {
1012  locBndExp = bndCondExp[i]->GetExp(j);
1013  id = locBndExp->as<LocalRegions::Expansion2D>()
1014  ->GetGeom2D()->GetFid();
1015  gid = FaceElmtGid[MeshFaceId.find(id)->second];
1016  bndOffset = bndCondExp[i]->GetCoeff_Offset(j) + bndTotal;
1017 
1018  // Since boundary information is defined to be aligned with
1019  // the geometry just use forward/forward (both coordinate
1020  // directions) defintiion for gid's
1021  for(k = 0; k < locBndExp->GetNcoeffs(); ++k)
1022  {
1023  m_bndCondCoeffsToGlobalCoeffsMap[bndOffset+k] = gid + k;
1024  }
1025  }
1026 
1027  nbndexp += bndCondExp[i]->GetExpSize();
1028  bndTotal += bndCondExp[i]->GetNcoeffs();
1029  }
1030 
1031  m_numGlobalBndCoeffs = trace->GetNcoeffs();
1032  m_numGlobalCoeffs = m_numGlobalBndCoeffs;
1033 
1034  CalculateBndSystemBandWidth();
1035 
1036  if( m_solnType == eDirectMultiLevelStaticCond && nGraphVerts )
1037  {
1038  if(m_staticCondLevel < (bottomUpGraph->GetNlevels()-1))
1039  {
1040  Array<OneD, int> vwgts_perm(nGraphVerts);
1041 
1042  for(int i = 0; i < nGraphVerts; i++)
1043  {
1044  vwgts_perm[i] = vwgts[perm[i]];
1045  }
1046 
1047  bottomUpGraph->ExpandGraphWithVertexWeights(vwgts_perm);
1048 
1049  m_nextLevelLocalToGlobalMap = MemoryManager<AssemblyMap>::
1050  AllocateSharedPtr(this,bottomUpGraph);
1051  }
1052  }
1053 
1054  cnt = 0;
1055  m_bndCondTraceToGlobalTraceMap = Array<OneD, int>(nbndexp);
1056  for(i = 0; i < bndCondExp.num_elements(); ++i)
1057  {
1058  for(j = 0; j < bndCondExp[i]->GetExpSize(); ++j)
1059  {
1060  locBndExp = bndCondExp[i]->GetExp(j);
1061  FaceGeom = locBndExp->as<LocalRegions::Expansion2D>()
1062  ->GetGeom2D();
1063  id = FaceGeom->GetFid();
1064  m_bndCondTraceToGlobalTraceMap[cnt++] =
1065  MeshFaceId.find(id)->second;
1066  }
1067  }
1068 
1069  // Now set up mapping from global coefficients to universal.
1070  ExpListSharedPtr tr = boost::dynamic_pointer_cast<ExpList>(trace);
1071  SetUpUniversalDGMap (locExp);
1072  SetUpUniversalTraceMap(locExp, tr, periodicFaces);
1073 
1074  m_hash = boost::hash_range(m_localToGlobalBndMap.begin(),
1075  m_localToGlobalBndMap.end());
1076  }
1077 
1078  /**
1079  * Constructs a mapping between the process-local global numbering and
1080  * a universal numbering of the trace space expansion. The universal
1081  * numbering is defined by the mesh edge IDs to enforce consistency
1082  * across processes.
1083  *
1084  * @param locExp List of local elemental expansions.
1085  */
1086  void AssemblyMapDG::SetUpUniversalDGMap(const ExpList &locExp)
1087  {
1088  StdRegions::StdExpansionSharedPtr locExpansion;
1089  int eid = 0;
1090  int cnt = 0;
1091  int id = 0;
1092  int order_e = 0;
1093  int vGlobalId = 0;
1094  int maxDof = 0;
1095  int dof = 0;
1096  int nDim = 0;
1097  int i,j,k;
1098 
1099  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1100 
1101  // Initialise the global to universal maps.
1102  m_globalToUniversalBndMap = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
1103  m_globalToUniversalBndMapUnique = Nektar::Array<OneD, int>(m_numGlobalBndCoeffs, -1);
1104 
1105  // Loop over all the elements in the domain and compute max
1106  // DOF. Reduce across all processes to get universal maximum.
1107  for(i = 0; i < locExpVector.size(); ++i)
1108  {
1109  locExpansion = boost::dynamic_pointer_cast<
1110  StdRegions::StdExpansion>(locExpVector[i]);
1111  nDim = locExpansion->GetShapeDimension();
1112 
1113  // Loop over all edges of element i
1114  if (nDim == 1)
1115  {
1116  maxDof = (1 > maxDof ? 1 : maxDof);
1117  }
1118  else if (nDim == 2)
1119  {
1120  for (j = 0; j < locExpansion->GetNedges(); ++j)
1121  {
1122  dof = locExpansion->GetEdgeNcoeffs(j);
1123  maxDof = (dof > maxDof ? dof : maxDof);
1124  }
1125  }
1126  else if (nDim == 3)
1127  {
1128  for (j = 0; j < locExpansion->GetNfaces(); ++j)
1129  {
1130  dof = locExpansion->GetFaceNcoeffs(j);
1131  maxDof = (dof > maxDof ? dof : maxDof);
1132  }
1133  }
1134  }
1135  m_comm->AllReduce(maxDof, LibUtilities::ReduceMax);
1136 
1137  // Now have trace edges Gid position
1138  cnt = 0;
1139  for(i = 0; i < locExpVector.size(); ++i)
1140  {
1141  locExpansion = boost::dynamic_pointer_cast<
1142  StdRegions::StdExpansion>(locExpVector[i]);
1143  nDim = locExpansion->GetShapeDimension();
1144 
1145  // Order list according to m_offset_elmt_id details in Exp
1146  // so that triangules are listed first and then quads
1147  eid = locExp.GetOffset_Elmt_Id(i);
1148 
1149  // Populate mapping for each edge of the element.
1150  if (nDim == 1)
1151  {
1152  int nverts = locExpansion->GetNverts();
1153  for(j = 0; j < nverts; ++j)
1154  {
1155  LocalRegions::PointExpSharedPtr locPointExp =
1156  m_elmtToTrace[eid][j]->as<LocalRegions::PointExp>();
1157  id = locPointExp->GetGeom()->GetGlobalID();
1158  vGlobalId = m_localToGlobalBndMap[cnt+j];
1159  m_globalToUniversalBndMap[vGlobalId]
1160  = id * maxDof + j + 1;
1161  }
1162  cnt += nverts;
1163  }
1164  else if (nDim == 2)
1165  {
1166  for(j = 0; j < locExpansion->GetNedges(); ++j)
1167  {
1168  LocalRegions::SegExpSharedPtr locSegExp =
1169  m_elmtToTrace[eid][j]->as<LocalRegions::SegExp>();
1170 
1171  id = locSegExp->GetGeom1D()->GetEid();
1172  order_e = locExpVector[eid]->GetEdgeNcoeffs(j);
1173 
1174  map<int,int> orientMap;
1175  Array<OneD, unsigned int> map1(order_e), map2(order_e);
1176  Array<OneD, int> sign1(order_e), sign2(order_e);
1177 
1178  locExpVector[eid]->GetEdgeToElementMap(j, StdRegions::eForwards, map1, sign1);
1179  locExpVector[eid]->GetEdgeToElementMap(j, locExpVector[eid]->GetEorient(j), map2, sign2);
1180 
1181  for (k = 0; k < map1.num_elements(); ++k)
1182  {
1183  // Find the elemental co-efficient in the original
1184  // mapping.
1185  int idx = -1;
1186  for (int l = 0; l < map2.num_elements(); ++l)
1187  {
1188  if (map1[k] == map2[l])
1189  {
1190  idx = l;
1191  break;
1192  }
1193  }
1194 
1195  ASSERTL2(idx != -1, "Problem with face to element map!");
1196  orientMap[k] = idx;
1197  }
1198 
1199  for(k = 0; k < order_e; ++k)
1200  {
1201  vGlobalId = m_localToGlobalBndMap[k+cnt];
1202  m_globalToUniversalBndMap[vGlobalId]
1203  = id * maxDof + orientMap[k] + 1;
1204  }
1205  cnt += order_e;
1206  }
1207  }
1208  else if (nDim == 3)
1209  {
1210  for(j = 0; j < locExpansion->GetNfaces(); ++j)
1211  {
1212  LocalRegions::Expansion2DSharedPtr locFaceExp =
1213  m_elmtToTrace[eid][j]
1214  ->as<LocalRegions::Expansion2D>();
1215 
1216  id = locFaceExp->GetGeom2D()->GetFid();
1217  order_e = locExpVector[eid]->GetFaceNcoeffs(j);
1218 
1219  map<int,int> orientMap;
1220  Array<OneD, unsigned int> map1(order_e), map2(order_e);
1221  Array<OneD, int> sign1(order_e), sign2(order_e);
1222 
1223  locExpVector[eid]->GetFaceToElementMap(j, StdRegions::eDir1FwdDir1_Dir2FwdDir2, map1, sign1);
1224  locExpVector[eid]->GetFaceToElementMap(j, locExpVector[eid]->GetFaceOrient(j), map2, sign2);
1225 
1226  for (k = 0; k < map1.num_elements(); ++k)
1227  {
1228  // Find the elemental co-efficient in the original
1229  // mapping.
1230  int idx = -1;
1231  for (int l = 0; l < map2.num_elements(); ++l)
1232  {
1233  if (map1[k] == map2[l])
1234  {
1235  idx = l;
1236  break;
1237  }
1238  }
1239 
1240  ASSERTL2(idx != -1, "Problem with face to element map!");
1241  orientMap[k] = idx;
1242  }
1243 
1244  for(k = 0; k < order_e; ++k)
1245  {
1246  vGlobalId = m_localToGlobalBndMap[k+cnt];
1247  m_globalToUniversalBndMap[vGlobalId]
1248  = id * maxDof + orientMap[k] + 1;
1249  }
1250  cnt += order_e;
1251  }
1252  }
1253  }
1254 
1255  // Initialise GSlib and populate the unique map.
1256  Array<OneD, long> tmp(m_globalToUniversalBndMap.num_elements());
1257  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
1258  {
1259  tmp[i] = m_globalToUniversalBndMap[i];
1260  }
1261  m_bndGsh = m_gsh = Gs::Init(tmp, m_comm);
1262  Gs::Unique(tmp, m_comm);
1263  for (i = 0; i < m_globalToUniversalBndMap.num_elements(); ++i)
1264  {
1265  m_globalToUniversalBndMapUnique[i] = (tmp[i] >= 0 ? 1 : 0);
1266  }
1267  }
1268 
1269  void AssemblyMapDG::SetUpUniversalTraceMap(
1270  const ExpList &locExp,
1271  const ExpListSharedPtr trace,
1272  const PeriodicMap &perMap)
1273  {
1274  Array<OneD, int> tmp;
1275  StdRegions::StdExpansionSharedPtr locExpansion;
1276  int i;
1277  int maxQuad = 0, quad = 0, nDim = 0, eid = 0, offset = 0;
1278 
1279  const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
1280 
1281  int nTracePhys = trace->GetTotPoints();
1282 
1283  // Initialise the trace to universal maps.
1284  m_traceToUniversalMap =
1285  Nektar::Array<OneD, int>(nTracePhys, -1);
1286  m_traceToUniversalMapUnique =
1287  Nektar::Array<OneD, int>(nTracePhys, -1);
1288 
1289  // Assume that each element of the expansion is of the same
1290  // dimension.
1291  nDim = locExpVector[0]->GetShapeDimension();
1292 
1293  if (nDim == 1)
1294  {
1295  maxQuad = (1 > maxQuad ? 1 : maxQuad);
1296  }
1297  else
1298  {
1299  for (i = 0; i < trace->GetExpSize(); ++i)
1300  {
1301  quad = trace->GetExp(i)->GetTotPoints();
1302  if (quad > maxQuad)
1303  {
1304  maxQuad = quad;
1305  }
1306  }
1307  }
1308  m_comm->AllReduce(maxQuad, LibUtilities::ReduceMax);
1309 
1310  if (nDim == 1)
1311  {
1312  for (int i = 0; i < trace->GetExpSize(); ++i)
1313  {
1314  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
1315  offset = trace->GetPhys_Offset(i);
1316 
1317 #if 1
1318  // Check to see if this vert is periodic. If it is, then we
1319  // need use the unique eid of the two points
1320  PeriodicMap::const_iterator it = perMap.find(eid);
1321  if (perMap.count(eid) > 0)
1322  {
1323  PeriodicEntity ent = it->second[0];
1324  if (ent.isLocal == false) // Not sure if true in 1D
1325  {
1326  eid = min(eid, ent.id);
1327  }
1328  }
1329 #endif
1330 
1331  m_traceToUniversalMap[offset] = eid*maxQuad+1;
1332  }
1333  }
1334  else
1335  {
1336  for (int i = 0; i < trace->GetExpSize(); ++i)
1337  {
1338  eid = trace->GetExp(i)->GetGeom()->GetGlobalID();
1339  offset = trace->GetPhys_Offset(i);
1340  quad = trace->GetExp(i)->GetTotPoints();
1341 
1342  // Check to see if this edge is periodic. If it is, then we
1343  // need to reverse the trace order of one edge only in the
1344  // universal map so that the data are reversed w.r.t each
1345  // other. We do this by using the minimum of the two IDs.
1346  PeriodicMap::const_iterator it = perMap.find(eid);
1347  bool realign = false;
1348  if (perMap.count(eid) > 0)
1349  {
1350  PeriodicEntity ent = it->second[0];
1351  if (ent.isLocal == false)
1352  {
1353  realign = eid == min(eid, ent.id);
1354  eid = min(eid, ent.id);
1355  }
1356  }
1357 
1358  for (int j = 0; j < quad; ++j)
1359  {
1360  m_traceToUniversalMap[j+offset] = eid*maxQuad+j+1;
1361  }
1362 
1363  if (realign)
1364  {
1365  if (nDim == 2)
1366  {
1367  RealignTraceElement(
1368  tmp = m_traceToUniversalMap+offset,
1369  it->second[0].orient, quad);
1370  }
1371  else
1372  {
1373  RealignTraceElement(
1374  tmp = m_traceToUniversalMap+offset,
1375  it->second[0].orient,
1376  trace->GetExp(i)->GetNumPoints(0),
1377  trace->GetExp(i)->GetNumPoints(1));
1378  }
1379  }
1380  }
1381  }
1382 
1383  Array<OneD, long> tmp2(nTracePhys);
1384  for (int i = 0; i < nTracePhys; ++i)
1385  {
1386  tmp2[i] = m_traceToUniversalMap[i];
1387  }
1388  m_traceGsh = Gs::Init(tmp2, m_comm);
1389  Gs::Unique(tmp2, m_comm);
1390  for (int i = 0; i < nTracePhys; ++i)
1391  {
1392  m_traceToUniversalMapUnique[i] = tmp2[i];
1393  }
1394  }
1395 
1396  void AssemblyMapDG::RealignTraceElement(
1397  Array<OneD, int> &toAlign,
1398  StdRegions::Orientation orient,
1399  int nquad1,
1400  int nquad2)
1401  {
1402  if (orient == StdRegions::eBackwards)
1403  {
1404  ASSERTL1(nquad2 == 0, "nquad2 != 0 for reorienation");
1405  for (int i = 0; i < nquad1/2; ++i)
1406  {
1407  swap(toAlign[i], toAlign[nquad1-1-i]);
1408  }
1409  }
1410  else if (orient != StdRegions::eForwards)
1411  {
1412  ASSERTL1(nquad2 != 0, "nquad2 == 0 for reorienation");
1413 
1414  Array<OneD, int> tmp(nquad1*nquad2);
1415 
1416  // Copy transpose.
1417  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
1418  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
1419  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
1420  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1421  {
1422  for (int i = 0; i < nquad2; ++i)
1423  {
1424  for (int j = 0; j < nquad1; ++j)
1425  {
1426  tmp[i*nquad1 + j] = toAlign[j*nquad2 + i];
1427  }
1428  }
1429  }
1430  else
1431  {
1432  for (int i = 0; i < nquad2; ++i)
1433  {
1434  for (int j = 0; j < nquad1; ++j)
1435  {
1436  tmp[i*nquad1 + j] = toAlign[i*nquad1 + j];
1437  }
1438  }
1439  }
1440 
1441  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
1442  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
1443  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
1444  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1445  {
1446  // Reverse x direction
1447  for (int i = 0; i < nquad2; ++i)
1448  {
1449  for (int j = 0; j < nquad1/2; ++j)
1450  {
1451  swap(tmp[i*nquad1 + j],
1452  tmp[i*nquad1 + nquad1-j-1]);
1453  }
1454  }
1455  }
1456 
1457  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
1458  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
1459  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
1460  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
1461  {
1462  // Reverse y direction
1463  for (int j = 0; j < nquad1; ++j)
1464  {
1465  for (int i = 0; i < nquad2/2; ++i)
1466  {
1467  swap(tmp[i*nquad1 + j],
1468  tmp[(nquad2-i-1)*nquad1 + j]);
1469  }
1470  }
1471  }
1472  Vmath::Vcopy(nquad1*nquad2, tmp, 1, toAlign, 1);
1473  }
1474  }
1475 
1476  void AssemblyMapDG::UniversalTraceAssemble(
1477  Array<OneD, NekDouble> &pGlobal) const
1478  {
1479  Gs::Gather(pGlobal, Gs::gs_add, m_traceGsh);
1480  }
1481 
1482  int AssemblyMapDG::v_GetLocalToGlobalMap(const int i) const
1483  {
1484  return m_localToGlobalBndMap[i];
1485  }
1486 
1487  int AssemblyMapDG::v_GetGlobalToUniversalMap(const int i) const
1488  {
1489  return m_globalToUniversalBndMap[i];
1490  }
1491 
1492  int AssemblyMapDG::v_GetGlobalToUniversalMapUnique(const int i) const
1493  {
1494  return m_globalToUniversalBndMapUnique[i];
1495  }
1496 
1497  const Array<OneD,const int>& AssemblyMapDG::v_GetLocalToGlobalMap()
1498  {
1499  return m_localToGlobalBndMap;
1500  }
1501 
1502  const Array<OneD,const int>& AssemblyMapDG::v_GetGlobalToUniversalMap()
1503  {
1504  return m_globalToUniversalBndMap;
1505  }
1506 
1507  const Array<OneD,const int>& AssemblyMapDG::v_GetGlobalToUniversalMapUnique()
1508  {
1509  return m_globalToUniversalBndMapUnique;
1510  }
1511 
1512  NekDouble AssemblyMapDG::v_GetLocalToGlobalSign(
1513  const int i) const
1514  {
1515  return GetLocalToGlobalBndSign(i);
1516  }
1517 
1518  void AssemblyMapDG::v_LocalToGlobal(
1519  const Array<OneD, const NekDouble>& loc,
1520  Array<OneD, NekDouble>& global) const
1521  {
1522  AssembleBnd(loc,global);
1523  }
1524 
1525  void AssemblyMapDG::v_LocalToGlobal(
1526  const NekVector<NekDouble>& loc,
1527  NekVector< NekDouble>& global) const
1528  {
1529  AssembleBnd(loc,global);
1530  }
1531 
1532  void AssemblyMapDG::v_GlobalToLocal(
1533  const Array<OneD, const NekDouble>& global,
1534  Array<OneD, NekDouble>& loc) const
1535  {
1536  GlobalToLocalBnd(global,loc);
1537  }
1538 
1539  void AssemblyMapDG::v_GlobalToLocal(
1540  const NekVector<NekDouble>& global,
1541  NekVector< NekDouble>& loc) const
1542  {
1543  GlobalToLocalBnd(global,loc);
1544  }
1545 
1546  void AssemblyMapDG::v_Assemble(
1547  const Array<OneD, const NekDouble> &loc,
1548  Array<OneD, NekDouble> &global) const
1549  {
1550  AssembleBnd(loc,global);
1551  }
1552 
1553  void AssemblyMapDG::v_Assemble(
1554  const NekVector<NekDouble>& loc,
1555  NekVector< NekDouble>& global) const
1556  {
1557  AssembleBnd(loc,global);
1558  }
1559 
1560  void AssemblyMapDG::v_UniversalAssemble(
1561  Array<OneD, NekDouble>& pGlobal) const
1562  {
1563  Gs::Gather(pGlobal, Gs::gs_add, m_gsh);
1564  }
1565 
1566  void AssemblyMapDG::v_UniversalAssemble(
1567  NekVector< NekDouble>& pGlobal) const
1568  {
1569  UniversalAssemble(pGlobal.GetPtr());
1570  }
1571 
1572  int AssemblyMapDG::v_GetFullSystemBandWidth() const
1573  {
1574  return GetBndSystemBandWidth();
1575  }
1576 
1577  int AssemblyMapDG::GetTraceToUniversalMap(int i)
1578  {
1579  return m_traceToUniversalMap[i];
1580  }
1581 
1582  int AssemblyMapDG::GetTraceToUniversalMapUnique(int i)
1583  {
1584  return m_traceToUniversalMapUnique[i];
1585  }
1586 
1587  int AssemblyMapDG::GetNumDirichletBndPhys()
1588  {
1589  return m_numDirichletBndPhys;
1590  }
1591 
1592  Array<OneD, StdRegions::StdExpansionSharedPtr>&
1593  AssemblyMapDG::GetElmtToTrace(const int i)
1594  {
1595  ASSERTL1(i >= 0 && i < m_elmtToTrace.num_elements(),
1596  "i is out of range");
1597  return m_elmtToTrace[i];
1598  }
1599 
1600  Array<OneD, Array< OneD, StdRegions::StdExpansionSharedPtr> >&
1601  AssemblyMapDG::GetElmtToTrace()
1602  {
1603  return m_elmtToTrace;
1604  }
1605 
1606 
1607  } //namespace
1608 } // namespace
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2