Nektar++
DisContField3D.cpp
Go to the documentation of this file.
1  //////////////////////////////////////////////////////////////////////////////
2  //
3  // File DisContField3D.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: Field definition for 3D domain with boundary
33  // conditions using LDG flux
34  //
35  ///////////////////////////////////////////////////////////////////////////////
36 
41  #include <LocalRegions/HexExp.h>
42  #include <LocalRegions/TetExp.h>
43  #include <LocalRegions/PrismExp.h>
44 
45  #include <boost/assign/std/vector.hpp>
46 
47  using namespace boost::assign;
48 
49  namespace Nektar
50  {
51  namespace MultiRegions
52  {
53  /**
54  * @class DisContField3D
55  * Abstraction of a global discontinuous three-dimensional spectral/hp
56  * element expansion which approximates the solution of a set of
57  * partial differential equations.
58  */
59 
60  /**
61  * @brief Default constructor.
62  */
63  DisContField3D::DisContField3D() :
64  ExpList3D (),
65  m_bndCondExpansions (),
66  m_bndConditions (),
67  m_trace(NullExpListSharedPtr)
68  {
69  }
70 
71  /**
72  * @brief Constructs a global discontinuous field based on an input mesh
73  * with boundary conditions.
74  */
78  const std::string &variable,
79  const bool SetUpJustDG)
80  : ExpList3D (pSession,graph3D,variable),
81  m_bndCondExpansions(),
82  m_bndConditions (),
83  m_trace(NullExpListSharedPtr)
84  {
85  if(variable.compare("DefaultVar") != 0) // do not set up BCs if default variable
86  {
88 
89  GenerateBoundaryConditionExpansion(graph3D,bcs,variable);
90  EvaluateBoundaryConditions(0.0, variable);
91 
92  // Find periodic edges for this variable.
93  FindPeriodicFaces(bcs, variable);
94  }
95 
96  if(SetUpJustDG)
97  {
98  SetUpDG();
99  }
100  else
101  {
102  // Set element edges to point to Robin BC edges if required.
103  int i,cnt,f;
104  Array<OneD, int> ElmtID, FaceID;
105  GetBoundaryToElmtMap(ElmtID, FaceID);
106 
107  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
108  {
110  locExpList = m_bndCondExpansions[i];
111 
112  for(f = 0; f < locExpList->GetExpSize(); ++f)
113  {
115  = (*m_exp)[ElmtID[cnt+f]]->
116  as<LocalRegions::Expansion3D>();
118  = locExpList->GetExp(f)->
119  as<LocalRegions::Expansion2D>();
120 
121  exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
122  exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
123  }
124  cnt += m_bndCondExpansions[i]->GetExpSize();
125  }
126  }
127  }
128 
129  /*
130  * @brief Copy type constructor which declares new boundary conditions
131  * and re-uses mapping info and trace space if possible
132  */
134  const DisContField3D &In,
135  const SpatialDomains::MeshGraphSharedPtr &graph3D,
136  const std::string &variable,
137  const bool SetUpJustDG)
138  : ExpList3D(In),
139  m_trace(NullExpListSharedPtr)
140  {
142 
143  GenerateBoundaryConditionExpansion(graph3D,bcs,variable);
144  EvaluateBoundaryConditions(0.0, variable);
145  ApplyGeomInfo();
146 
148  {
149  // Find periodic edges for this variable.
150  FindPeriodicFaces(bcs, variable);
151 
152  if (SetUpJustDG)
153  {
154  SetUpDG(variable);
155  }
156  else
157  {
158  int i,cnt,f;
159  Array<OneD, int> ElmtID,FaceID;
160  GetBoundaryToElmtMap(ElmtID,FaceID);
161 
162  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
163  {
165  locExpList = m_bndCondExpansions[i];
166 
167  for(f = 0; f < locExpList->GetExpSize(); ++f)
168  {
170  = (*m_exp)[ElmtID[cnt+f]]->
171  as<LocalRegions::Expansion3D>();
173  = locExpList->GetExp(f)->
174  as<LocalRegions::Expansion2D>();
175 
176  exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
177  exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
178  }
179 
180  cnt += m_bndCondExpansions[i]->GetExpSize();
181  }
183  }
184 
185  }
186  //else if we have the same boundary condition
187  else
188  {
189  if(SetUpJustDG)
190  {
192  m_trace = In.m_trace;
193  m_traceMap = In.m_traceMap;
197  }
198  else
199  {
201  m_trace = In.m_trace;
202  m_traceMap = In.m_traceMap;
206 
207  int i,cnt,f;
208  Array<OneD, int> ElmtID,FaceID;
209  GetBoundaryToElmtMap(ElmtID,FaceID);
210 
211  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
212  {
214  locExpList = m_bndCondExpansions[i];
215 
216  for(f = 0; f < locExpList->GetExpSize(); ++f)
217  {
219  = (*m_exp)[ElmtID[cnt+f]]->
220  as<LocalRegions::Expansion3D>();
222  = locExpList->GetExp(f)->
223  as<LocalRegions::Expansion2D>();
224 
225  exp3d->SetFaceExp(FaceID[cnt+f],exp2d);
226  exp2d->SetAdjacentElementExp(FaceID[cnt+f],exp3d);
227  }
228 
229  cnt += m_bndCondExpansions[i]->GetExpSize();
230  }
231 
232  if(m_session->DefinesSolverInfo("PROJECTION"))
233  {
234  std::string ProjectStr =
235  m_session->GetSolverInfo("PROJECTION");
236  if (ProjectStr == "MixedCGDG" ||
237  ProjectStr == "Mixed_CG_Discontinuous")
238  {
239  SetUpDG(variable);
240  }
241  else
242  {
244  }
245  }
246  else
247  {
249  }
250  }
251  }
252  }
253 
254  /**
255  *
256  */
258  ExpList3D(In),
259  m_bndCondExpansions (In.m_bndCondExpansions),
260  m_bndConditions (In.m_bndConditions),
261  m_globalBndMat (In.m_globalBndMat),
262  m_trace (In.m_trace),
263  m_traceMap (In.m_traceMap),
264  m_periodicFaces (In.m_periodicFaces),
265  m_periodicEdges (In.m_periodicEdges),
266  m_periodicVerts (In.m_periodicVerts)
267  {
268  }
269 
270  /**
271  * @brief Destructor.
272  */
274  {
275  }
276 
278  const GlobalLinSysKey &mkey)
279  {
281  "Routine currently only tested for HybridDGHelmholtz");
282  ASSERTL1(mkey.GetGlobalSysSolnType() ==
283  m_traceMap->GetGlobalSysSolnType(),
284  "The local to global map is not set up for the requested "
285  "solution type");
286 
287  GlobalLinSysSharedPtr glo_matrix;
288  GlobalLinSysMap::iterator matrixIter = m_globalBndMat->find(mkey);
289 
290  if (matrixIter == m_globalBndMat->end())
291  {
292  glo_matrix = GenGlobalBndLinSys(mkey,m_traceMap);
293  (*m_globalBndMat)[mkey] = glo_matrix;
294  }
295  else
296  {
297  glo_matrix = matrixIter->second;
298  }
299 
300  return glo_matrix;
301  }
302 
303  /**
304  * @brief Set up all DG member variables and maps.
305  */
306  void DisContField3D::SetUpDG(const std::string variable)
307  {
309  {
310  return;
311  }
312 
313  ExpList2DSharedPtr trace;
314 
316  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph3D>(
317  m_graph);
318 
319  // Set up matrix map
322 
323  // Set up Trace space
324  bool UseGenSegExp = true;
327  *m_exp,graph3D, m_periodicFaces, UseGenSegExp);
328 
329  m_trace = trace;
331  m_session,graph3D,trace,*this,m_bndCondExpansions,
332  m_bndConditions, m_periodicFaces,variable);
333 
335  &elmtToTrace = m_traceMap->GetElmtToTrace();
336 
337  // Scatter trace segments to 3D elements. For each element, we find
338  // the trace segment associated to each edge. The element then
339  // retains a pointer to the trace space segments, to ensure
340  // uniqueness of normals when retrieving from two adjoining elements
341  // which do not lie in a plane.
342  for (int i = 0; i < m_exp->size(); ++i)
343  {
344  for (int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j)
345  {
347  (*m_exp)[i]->as<LocalRegions::Expansion3D>();
349  elmtToTrace[i][j]->as<LocalRegions::Expansion2D>();
350  exp3d->SetFaceExp (j, exp2d);
351  exp2d->SetAdjacentElementExp(j, exp3d);
352  }
353  }
354 
355  // Set up physical normals
357 
358  // Set up information for parallel jobs.
359  for (int i = 0; i < m_trace->GetExpSize(); ++i)
360  {
362  m_trace->GetExp(i)->as<LocalRegions::Expansion2D>();
363 
364  int offset = m_trace->GetPhys_Offset(i);
365  int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
367  traceGeomId);
368 
369  if (pIt != m_periodicFaces.end() && !pIt->second[0].isLocal)
370  {
371  if (traceGeomId != min(pIt->second[0].id, traceGeomId))
372  {
373  traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
374  traceEl->GetLeftAdjacentElementFace());
375  }
376  }
377  else if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
378  {
379  traceEl->GetLeftAdjacentElementExp()->NegateFaceNormal(
380  traceEl->GetLeftAdjacentElementFace());
381  }
382  }
383 
384  int cnt, n, e;
385 
386  // Identify boundary faces
387  for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
388  {
389  if (m_bndConditions[n]->GetBoundaryConditionType() !=
391  {
392  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
393  {
394  m_boundaryFaces.insert(
395  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
396  }
397  }
398  cnt += m_bndCondExpansions[n]->GetExpSize();
399  }
400 
401  // Set up information for periodic boundary conditions.
402  boost::unordered_map<int,pair<int,int> > perFaceToExpMap;
403  boost::unordered_map<int,pair<int,int> >::iterator it2;
404  cnt = 0;
406  for (int n = 0; n < m_exp->size(); ++n)
407  {
408  exp3d = (*m_exp)[n]->as<LocalRegions::Expansion3D>();
409  for (int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
410  {
412  exp3d->GetGeom3D()->GetFid(e));
413 
414  if (it != m_periodicFaces.end())
415  {
416  perFaceToExpMap[it->first] = make_pair(n, e);
417  }
418  }
419  }
420 
421  // Set up left-adjacent face list.
422  m_leftAdjacentFaces.resize(cnt);
423  cnt = 0;
424  for (int i = 0; i < m_exp->size(); ++i)
425  {
426  for (int j = 0; j < (*m_exp)[i]->GetNfaces(); ++j, ++cnt)
427  {
429  }
430  }
431 
432  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
433  cnt = 0;
434  for (int n = 0; n < m_exp->size(); ++n)
435  {
436  exp3d = (*m_exp)[n]->as<LocalRegions::Expansion3D>();
437  for (int e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
438  {
439  int faceGeomId = exp3d->GetGeom3D()->GetFid(e);
440  int offset = m_trace->GetPhys_Offset(
441  elmtToTrace[n][e]->GetElmtId());
442 
443  // Check to see if this face is periodic.
444  PeriodicMap::iterator it = m_periodicFaces.find(faceGeomId);
445 
446  if (it != m_periodicFaces.end())
447  {
448  const PeriodicEntity &ent = it->second[0];
449  it2 = perFaceToExpMap.find(ent.id);
450 
451  if (it2 == perFaceToExpMap.end())
452  {
453  if (m_session->GetComm()->GetSize() > 1 &&
454  !ent.isLocal)
455  {
456  continue;
457  }
458  else
459  {
460  ASSERTL1(false, "Periodic edge not found!");
461  }
462  }
463 
465  "Periodic face in non-forward space?");
466 
467  int offset2 = m_trace->GetPhys_Offset(
468  elmtToTrace[it2->second.first][it2->second.second]->
469  GetElmtId());
470 
471  // Calculate relative orientations between faces to
472  // calculate copying map.
473  int nquad1 = elmtToTrace[n][e]->GetNumPoints(0);
474  int nquad2 = elmtToTrace[n][e]->GetNumPoints(1);
475 
476  vector<int> tmpBwd(nquad1*nquad2);
477  vector<int> tmpFwd(nquad1*nquad2);
478 
483  {
484  for (int i = 0; i < nquad2; ++i)
485  {
486  for (int j = 0; j < nquad1; ++j)
487  {
488  tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
489  tmpFwd[i*nquad1+j] = offset + j*nquad2+i;
490  }
491  }
492  }
493  else
494  {
495  for (int i = 0; i < nquad2; ++i)
496  {
497  for (int j = 0; j < nquad1; ++j)
498  {
499  tmpBwd[i*nquad1+j] = offset2 + i*nquad1+j;
500  tmpFwd[i*nquad1+j] = offset + i*nquad1+j;
501  }
502  }
503  }
504 
509  {
510  // Reverse x direction
511  for (int i = 0; i < nquad2; ++i)
512  {
513  for (int j = 0; j < nquad1/2; ++j)
514  {
515  swap(tmpFwd[i*nquad1 + j],
516  tmpFwd[i*nquad1 + nquad1-j-1]);
517  }
518  }
519  }
520 
525  {
526  // Reverse y direction
527  for (int j = 0; j < nquad1; ++j)
528  {
529  for (int i = 0; i < nquad2/2; ++i)
530  {
531  swap(tmpFwd[i*nquad1 + j],
532  tmpFwd[(nquad2-i-1)*nquad1 + j]);
533  }
534  }
535  }
536 
537  for (int i = 0; i < nquad1*nquad2; ++i)
538  {
539  m_periodicFwdCopy.push_back(tmpFwd[i]);
540  m_periodicBwdCopy.push_back(tmpBwd[i]);
541  }
542  }
543  }
544  }
545  }
546 
547  /**
548  * For each boundary region, checks that the types and number of
549  * boundary expansions in that region match.
550  *
551  * @param In ContField3D to compare with.
552  * @return true if boundary conditions match.
553  */
555  const DisContField3D &In)
556  {
557  int i;
558  bool returnval = true;
559 
560  for(i = 0; i < m_bndConditions.num_elements(); ++i)
561  {
562 
563  // check to see if boundary condition type is the same
564  // and there are the same number of boundary
565  // conditions in the boundary definition.
566  if((m_bndConditions[i]->GetBoundaryConditionType()
567  != In.m_bndConditions[i]->GetBoundaryConditionType())||
569  != In.m_bndCondExpansions[i]->GetExpSize()))
570  {
571  returnval = false;
572  break;
573  }
574  }
575 
576  // Compare with all other processes. Return true only if all
577  // processes report having the same boundary conditions.
578  int vSame = returnval ? 1 : 0;
579  m_comm->AllReduce(vSame, LibUtilities::ReduceMin);
580 
581  return (vSame == 1);
582  }
583 
584  /**
585  * According to their boundary region, the separate segmental boundary
586  * expansions are bundled together in an object of the class
587  * MultiRegions#ExpList2D.
588  *
589  * \param graph3D A mesh, containing information about the domain and
590  * the spectral/hp element expansion.
591  * \param bcs An entity containing information about the boundaries and
592  * boundary conditions.
593  * \param variable An optional parameter to indicate for which variable
594  * the boundary conditions should be discretised.
595  */
597  const SpatialDomains::MeshGraphSharedPtr &graph3D,
599  const std::string &variable)
600  {
601  int cnt = 0;
604 
606  bcs.GetBoundaryRegions();
607  const SpatialDomains::BoundaryConditionCollection &bconditions =
608  bcs.GetBoundaryConditions();
609  SpatialDomains::BoundaryRegionCollection::const_iterator it;
610 
611  // count the number of non-periodic boundary regions
612  for (it = bregions.begin(); it != bregions.end(); ++it)
613  {
614  SpatialDomains::BoundaryConditionShPtr boundaryCondition =
615  GetBoundaryCondition(bconditions, it->first, variable);
616  if (boundaryCondition->GetBoundaryConditionType() !=
618  {
619  cnt++;
620  }
621  }
622 
625 
626  cnt = 0;
627 
628  // list Dirichlet boundaries first
629  for (it = bregions.begin(); it != bregions.end(); ++it)
630  {
631  locBCond = GetBoundaryCondition(
632  bconditions, it->first, variable);
633 
634  if(locBCond->GetBoundaryConditionType()
636  {
638  ::AllocateSharedPtr(m_session, *(it->second),
639  graph3D, variable);
640 
641  // Set up normals on non-Dirichlet boundary conditions
642  if(locBCond->GetBoundaryConditionType() !=
644  {
646  }
647 
648  m_bndCondExpansions[cnt] = locExpList;
649  m_bndConditions[cnt++] = locBCond;
650  }
651  }
652  }
653 
654  /**
655  * @brief Determine the periodic faces, edges and vertices for the given
656  * graph.
657  *
658  * @param bcs Information about the boundary conditions.
659  * @param variable Specifies the field.
660  */
663  const std::string &variable)
664  {
666  = bcs.GetBoundaryRegions();
668  = bcs.GetBoundaryConditions();
670  = boost::dynamic_pointer_cast<
672  SpatialDomains::BoundaryRegionCollection::const_iterator it;
673 
675  m_session->GetComm()->GetRowComm();
677  m_session->GetCompositeOrdering();
678  LibUtilities::BndRegionOrdering bndRegOrder =
679  m_session->GetBndRegionOrdering();
681  m_graph->GetComposites();
682 
683  // perComps: Stores a unique collection of pairs of periodic
684  // composites (i.e. if composites 1 and 2 are periodic then this map
685  // will contain either the pair (1,2) or (2,1) but not both).
686  //
687  // The four maps allVerts, allCoord, allEdges and allOrient map a
688  // periodic face to a vector containing the vertex ids of the face;
689  // their coordinates; the edge ids of the face; and their
690  // orientation within that face respectively.
691  //
692  // Finally the three sets locVerts, locEdges and locFaces store any
693  // vertices, edges and faces that belong to a periodic composite and
694  // lie on this process.
695  map<int,int> perComps;
696  map<int,vector<int> > allVerts;
697  map<int,SpatialDomains::PointGeomVector> allCoord;
698  map<int,vector<int> > allEdges;
699  map<int,vector<StdRegions::Orientation> > allOrient;
700  set<int> locVerts;
701  set<int> locEdges;
702  set<int> locFaces;
703 
704  int region1ID, region2ID, i, j, k, cnt;
706 
707  // Set up a set of all local verts and edges.
708  for(i = 0; i < (*m_exp).size(); ++i)
709  {
710  for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
711  {
712  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
713  locVerts.insert(id);
714  }
715 
716  for(j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
717  {
718  int id = (*m_exp)[i]->GetGeom()->GetEid(j);
719  locEdges.insert(id);
720  }
721  }
722 
723  // Begin by populating the perComps map. We loop over all periodic
724  // boundary conditions and determine the composite associated with
725  // it, then fill out the all* maps.
726  for (it = bregions.begin(); it != bregions.end(); ++it)
727  {
728  locBCond = GetBoundaryCondition(
729  bconditions, it->first, variable);
730 
731  if (locBCond->GetBoundaryConditionType()
733  {
734  continue;
735  }
736 
737  // Identify periodic boundary region IDs.
738  region1ID = it->first;
739  region2ID = boost::static_pointer_cast<
741  locBCond)->m_connectedBoundaryRegion;
742 
743  // Check the region only contains a single composite.
744  ASSERTL0(it->second->size() == 1,
745  "Boundary region "+boost::lexical_cast<string>(
746  region1ID)+" should only contain 1 composite.");
747 
748  // From this identify composites by looking at the original
749  // boundary region ordering. Note that in serial the mesh
750  // partitioner is not run, so this map will be empty and
751  // therefore needs to be populated by using the corresponding
752  // boundary region.
753  int cId1, cId2;
754  if (vComm->GetSize() == 1)
755  {
756  cId1 = it->second->begin()->first;
757  cId2 = bregions.find(region2ID)->second->begin()->first;
758  }
759  else
760  {
761  cId1 = bndRegOrder.find(region1ID)->second[0];
762  cId2 = bndRegOrder.find(region2ID)->second[0];
763  }
764 
765  SpatialDomains::Composite c = it->second->begin()->second;
766  vector<unsigned int> tmpOrder;
767 
768  // From the composite, we now construct the allVerts, allEdges
769  // and allCoord map so that they can be transferred across
770  // processors. We also populate the locFaces set to store a
771  // record of all faces local to this process.
772  for (i = 0; i < c->size(); ++i)
773  {
775  boost::dynamic_pointer_cast<
776  SpatialDomains::Geometry2D>((*c)[i]);
777  ASSERTL1(faceGeom, "Unable to cast to shared ptr");
778 
779  // Get geometry ID of this face and store in locFaces.
780  int faceId = (*c)[i]->GetGlobalID();
781  locFaces.insert(faceId);
782 
783  // In serial, mesh partitioning will not have occurred so
784  // need to fill composite ordering map manually.
785  if (vComm->GetSize() == 1)
786  {
787  tmpOrder.push_back((*c)[i]->GetGlobalID());
788  }
789 
790  // Loop over vertices and edges of the face to populate
791  // allVerts, allEdges and allCoord maps.
792  vector<int> vertList, edgeList;
794  vector<StdRegions::Orientation> orientVec;
795  for (j = 0; j < faceGeom->GetNumVerts(); ++j)
796  {
797  vertList .push_back(faceGeom->GetVid (j));
798  edgeList .push_back(faceGeom->GetEid (j));
799  coordVec .push_back(faceGeom->GetVertex(j));
800  orientVec.push_back(faceGeom->GetEorient(j));
801  }
802 
803  allVerts [faceId] = vertList;
804  allEdges [faceId] = edgeList;
805  allCoord [faceId] = coordVec;
806  allOrient[faceId] = orientVec;
807  }
808 
809  // In serial, record the composite ordering in compOrder for
810  // later in the routine.
811  if (vComm->GetSize() == 1)
812  {
813  compOrder[it->second->begin()->first] = tmpOrder;
814  }
815 
816  // See if we already have either region1 or region2 stored in
817  // perComps map already and do a sanity check to ensure regions
818  // are mutually periodic.
819  if (perComps.count(cId1) == 0)
820  {
821  if (perComps.count(cId2) == 0)
822  {
823  perComps[cId1] = cId2;
824  }
825  else
826  {
827  std::stringstream ss;
828  ss << "Boundary region " << cId2 << " should be "
829  << "periodic with " << perComps[cId2] << " but "
830  << "found " << cId1 << " instead!";
831  ASSERTL0(perComps[cId2] == cId1, ss.str());
832  }
833  }
834  else
835  {
836  std::stringstream ss;
837  ss << "Boundary region " << cId1 << " should be "
838  << "periodic with " << perComps[cId1] << " but "
839  << "found " << cId2 << " instead!";
840  ASSERTL0(perComps[cId1] == cId1, ss.str());
841  }
842  }
843 
844  // The next routines process local face lists to exchange vertices,
845  // edges and faces.
846  int n = vComm->GetSize();
847  int p = vComm->GetRank();
848  int totFaces;
849  Array<OneD, int> facecounts(n,0);
850  Array<OneD, int> vertcounts(n,0);
851  Array<OneD, int> faceoffset(n,0);
852  Array<OneD, int> vertoffset(n,0);
853 
854  // First exchange the number of faces on each process.
855  facecounts[p] = locFaces.size();
856  vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
857 
858  // Set up an offset map to allow us to distribute face IDs to all
859  // processors.
860  faceoffset[0] = 0;
861  for (i = 1; i < n; ++i)
862  {
863  faceoffset[i] = faceoffset[i-1] + facecounts[i-1];
864  }
865 
866  // Calculate total number of faces.
867  totFaces = Vmath::Vsum(n, facecounts, 1);
868 
869  // faceIds holds face IDs for each periodic face. faceVerts holds
870  // the number of vertices in this face.
871  Array<OneD, int> faceIds (totFaces, 0);
872  Array<OneD, int> faceVerts(totFaces, 0);
873 
874  // Process p writes IDs of its faces into position faceoffset[p] of
875  // faceIds which allows us to perform an AllReduce to distribute
876  // information amongst processors.
877  set<int>::iterator sIt;
878  for (i = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
879  {
880  faceIds [faceoffset[p] + i ] = *sIt;
881  faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
882  }
883 
884  vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
885  vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
886 
887  // procVerts holds number of vertices (and also edges since each
888  // face is 2D) on each process.
889  Array<OneD, int> procVerts(n,0);
890  int nTotVerts;
891 
892  // Note if there are no periodic faces at all calling Vsum will
893  // cause a segfault.
894  if (totFaces > 0)
895  {
896  // Calculate number of vertices on each processor.
897  nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
898  }
899  else
900  {
901  nTotVerts = 0;
902  }
903 
904  for (i = 0; i < n; ++i)
905  {
906  if (facecounts[i] > 0)
907  {
908  procVerts[i] = Vmath::Vsum(
909  facecounts[i], faceVerts + faceoffset[i], 1);
910  }
911  else
912  {
913  procVerts[i] = 0;
914  }
915  }
916 
917  // vertoffset is defined in the same manner as edgeoffset
918  // beforehand.
919  vertoffset[0] = 0;
920  for (i = 1; i < n; ++i)
921  {
922  vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
923  }
924 
925  // At this point we exchange all vertex IDs, edge IDs and vertex
926  // coordinates for each face. The coordinates are necessary because
927  // we need to calculate relative face orientations between periodic
928  // faces to determined edge and vertex connectivity.
929  Array<OneD, int> vertIds(nTotVerts, 0);
930  Array<OneD, int> edgeIds(nTotVerts, 0);
931  Array<OneD, int> edgeOrt(nTotVerts, 0);
932  Array<OneD, NekDouble> vertX (nTotVerts, 0.0);
933  Array<OneD, NekDouble> vertY (nTotVerts, 0.0);
934  Array<OneD, NekDouble> vertZ (nTotVerts, 0.0);
935 
936  for (cnt = 0, sIt = locFaces.begin();
937  sIt != locFaces.end(); ++sIt)
938  {
939  for (j = 0; j < allVerts[*sIt].size(); ++j)
940  {
941  int vertId = allVerts[*sIt][j];
942  vertIds[vertoffset[p] + cnt ] = vertId;
943  vertX [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(0);
944  vertY [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(1);
945  vertZ [vertoffset[p] + cnt ] = (*allCoord[*sIt][j])(2);
946  edgeIds[vertoffset[p] + cnt ] = allEdges [*sIt][j];
947  edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
948  }
949  }
950 
951  vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
952  vComm->AllReduce(vertX, LibUtilities::ReduceSum);
953  vComm->AllReduce(vertY, LibUtilities::ReduceSum);
954  vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
955  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
956  vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
957 
958  // Finally now we have all of this information, we construct maps
959  // which make accessing the information easier. These are
960  // conceptually the same as all* maps at the beginning of the
961  // routine, but now hold information for all periodic vertices.
962  map<int, vector<int> > vertMap;
963  map<int, vector<int> > edgeMap;
964  map<int, SpatialDomains::PointGeomVector> coordMap;
965 
966  // These final two maps are required for determining the relative
967  // orientation of periodic edges. vCoMap associates vertex IDs with
968  // their coordinates, and eIdMap maps an edge ID to the two vertices
969  // which construct it.
970  map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
971  map<int, pair<int, int> > eIdMap;
972 
973  for (cnt = i = 0; i < totFaces; ++i)
974  {
975  vector<int> edges(faceVerts[i]);
976  vector<int> verts(faceVerts[i]);
977  SpatialDomains::PointGeomVector coord(faceVerts[i]);
978 
979  // Keep track of cnt to enable correct edge vertices to be
980  // inserted into eIdMap.
981  int tmp = cnt;
982  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
983  {
984  edges[j] = edgeIds[cnt];
985  verts[j] = vertIds[cnt];
988  3, verts[j], vertX[cnt], vertY[cnt], vertZ[cnt]);
989  vCoMap[vertIds[cnt]] = coord[j];
990 
991  // Try to insert edge into the eIdMap to avoid re-inserting.
992  pair<map<int, pair<int, int> >::iterator, bool> testIns =
993  eIdMap.insert(make_pair(
994  edgeIds[cnt],
995  make_pair(vertIds[tmp+j],
996  vertIds[tmp+((j+1) % faceVerts[i])])));
997 
998  if (testIns.second == false)
999  {
1000  continue;
1001  }
1002 
1003  // If the edge is reversed with respect to the face, then
1004  // swap the edges so that we have the original ordering of
1005  // the edge in the 3D element. This is necessary to properly
1006  // determine edge orientation.
1007  if ((StdRegions::Orientation)edgeOrt[cnt]
1009  {
1010  swap(testIns.first->second.first,
1011  testIns.first->second.second);
1012  }
1013  }
1014 
1015  vertMap [faceIds[i]] = verts;
1016  edgeMap [faceIds[i]] = edges;
1017  coordMap[faceIds[i]] = coord;
1018  }
1019 
1020  // Go through list of composites and figure out which edges are
1021  // parallel from original ordering in session file. This includes
1022  // composites which are not necessarily on this process.
1023  map<int,int>::iterator cIt, pIt;
1024  map<int,int>::const_iterator oIt;
1025 
1026  // Store temporary map of periodic vertices which will hold all
1027  // periodic vertices on the entire mesh so that doubly periodic
1028  // vertices/edges can be counted properly across partitions. Local
1029  // vertices/edges are copied into m_periodicVerts and
1030  // m_periodicEdges at the end of the function.
1031  PeriodicMap periodicVerts;
1032  PeriodicMap periodicEdges;
1033 
1034  // Construct two maps which determine how vertices and edges of
1035  // faces connect given a specific face orientation. The key of the
1036  // map is the number of vertices in the face, used to determine
1037  // difference between tris and quads.
1038  map<int, map<StdRegions::Orientation, vector<int> > > vmap;
1039  map<int, map<StdRegions::Orientation, vector<int> > > emap;
1040 
1041  map<StdRegions::Orientation, vector<int> > quadVertMap;
1042  quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] += 0,1,2,3;
1043  quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] += 3,2,1,0;
1044  quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] += 1,0,3,2;
1045  quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] += 2,3,0,1;
1046  quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] += 0,3,2,1;
1047  quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] += 1,2,3,0;
1048  quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] += 3,0,1,2;
1049  quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] += 2,1,0,3;
1050 
1051  map<StdRegions::Orientation, vector<int> > quadEdgeMap;
1052  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] += 0,1,2,3;
1053  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] += 2,1,0,3;
1054  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] += 0,3,2,1;
1055  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] += 2,3,0,1;
1056  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] += 3,2,1,0;
1057  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] += 1,2,3,0;
1058  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] += 3,0,1,2;
1059  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] += 1,0,3,2;
1060 
1061  map<StdRegions::Orientation, vector<int> > triVertMap;
1062  triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] += 0,1,2;
1063  triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] += 1,0,2;
1064 
1065  map<StdRegions::Orientation, vector<int> > triEdgeMap;
1066  triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] += 0,1,2;
1067  triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] += 0,2,1;
1068 
1069  vmap[3] = triVertMap;
1070  vmap[4] = quadVertMap;
1071  emap[3] = triEdgeMap;
1072  emap[4] = quadEdgeMap;
1073 
1074  map<int,int> allCompPairs;
1075 
1076  // Finally we have enough information to populate the periodic
1077  // vertex, edge and face maps. Begin by looping over all pairs of
1078  // periodic composites to determine pairs of periodic faces.
1079  for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
1080  {
1082  const int id1 = cIt->first;
1083  const int id2 = cIt->second;
1084  std::string id1s = boost::lexical_cast<string>(id1);
1085  std::string id2s = boost::lexical_cast<string>(id2);
1086 
1087  if (compMap.count(id1) > 0)
1088  {
1089  c[0] = compMap[id1];
1090  }
1091 
1092  if (compMap.count(id2) > 0)
1093  {
1094  c[1] = compMap[id2];
1095  }
1096 
1097  ASSERTL0(c[0] || c[1],
1098  "Neither composite not found on this process!");
1099 
1100  // Loop over composite ordering to construct list of all
1101  // periodic faces, regardless of whether they are on this
1102  // process.
1103  map<int,int> compPairs;
1104 
1105  ASSERTL0(compOrder.count(id1) > 0,
1106  "Unable to find composite "+id1s+" in order map.");
1107  ASSERTL0(compOrder.count(id2) > 0,
1108  "Unable to find composite "+id2s+" in order map.");
1109  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1110  "Periodic composites "+id1s+" and "+id2s+
1111  " should have the same number of elements.");
1112  ASSERTL0(compOrder[id1].size() > 0,
1113  "Periodic composites "+id1s+" and "+id2s+
1114  " are empty!");
1115 
1116  // Look up composite ordering to determine pairs.
1117  for (i = 0; i < compOrder[id1].size(); ++i)
1118  {
1119  int eId1 = compOrder[id1][i];
1120  int eId2 = compOrder[id2][i];
1121 
1122  ASSERTL0(compPairs.count(eId1) == 0,
1123  "Already paired.");
1124 
1125  // Sanity check that the faces are mutually periodic.
1126  if (compPairs.count(eId2) != 0)
1127  {
1128  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
1129  }
1130  compPairs[eId1] = eId2;
1131  }
1132 
1133  // Now that we have all pairs of periodic faces, loop over the
1134  // ones local to this process and populate face/edge/vertex
1135  // maps.
1136  for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1137  {
1138  int ids [2] = {pIt->first, pIt->second};
1139  bool local[2] = {locFaces.count(pIt->first) > 0,
1140  locFaces.count(pIt->second) > 0};
1141 
1142  ASSERTL0(coordMap.count(ids[0]) > 0 &&
1143  coordMap.count(ids[1]) > 0,
1144  "Unable to find face in coordinate map");
1145 
1146  allCompPairs[pIt->first ] = pIt->second;
1147  allCompPairs[pIt->second] = pIt->first;
1148 
1149  // Loop up coordinates of the faces, check they have the
1150  // same number of vertices.
1152  = { coordMap[ids[0]], coordMap[ids[1]] };
1153 
1154  ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
1155  "Two periodic faces have different number "
1156  "of vertices!");
1157 
1158  // o will store relative orientation of faces. Note that in
1159  // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
1160  // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
1161  // different going from face1->face2 instead of face2->face1
1162  // (check this).
1164 
1165  // Record periodic faces.
1166  for (i = 0; i < 2; ++i)
1167  {
1168  if (!local[i])
1169  {
1170  continue;
1171  }
1172 
1173  // Reference to the other face.
1174  int other = (i+1) % 2;
1175 
1176  // Calculate relative face orientation.
1177  if (tmpVec[0].size() == 3)
1178  {
1180  tmpVec[i], tmpVec[other]);
1181  }
1182  else
1183  {
1185  tmpVec[i], tmpVec[other]);
1186  }
1187 
1188  // Record face ID, orientation and whether other face is
1189  // local.
1190  PeriodicEntity ent(ids [other], o,
1191  local[other]);
1192  m_periodicFaces[ids[i]].push_back(ent);
1193  }
1194 
1195  int nFaceVerts = vertMap[ids[0]].size();
1196 
1197  // Determine periodic vertices.
1198  for (i = 0; i < 2; ++i)
1199  {
1200  int other = (i+1) % 2;
1201 
1202  // Calculate relative face orientation.
1203  if (tmpVec[0].size() == 3)
1204  {
1206  tmpVec[i], tmpVec[other]);
1207  }
1208  else
1209  {
1211  tmpVec[i], tmpVec[other]);
1212  }
1213 
1214  if (nFaceVerts == 3)
1215  {
1216  ASSERTL0(
1219  "Unsupported face orientation for face "+
1220  boost::lexical_cast<string>(ids[i]));
1221  }
1222 
1223  // Look up vertices for this face.
1224  vector<int> per1 = vertMap[ids[i]];
1225  vector<int> per2 = vertMap[ids[other]];
1226 
1227  // tmpMap will hold the pairs of vertices which are
1228  // periodic.
1229  map<int, pair<int, bool> > tmpMap;
1230  map<int, pair<int, bool> >::iterator mIt;
1231 
1232  // Use vmap to determine which vertices connect given
1233  // the orientation o.
1234  for (j = 0; j < nFaceVerts; ++j)
1235  {
1236  int v = vmap[nFaceVerts][o][j];
1237  tmpMap[per1[j]] = make_pair(
1238  per2[v], locVerts.count(per2[v]) > 0);
1239  }
1240 
1241  // Now loop over tmpMap to associate periodic vertices.
1242  for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1243  {
1244  PeriodicEntity ent2(mIt->second.first,
1246  mIt->second.second);
1247 
1248  // See if this vertex has been recorded already.
1249  PeriodicMap::iterator perIt = periodicVerts.find(
1250  mIt->first);
1251 
1252  if (perIt == periodicVerts.end())
1253  {
1254  // Vertex is new - just record this entity as
1255  // usual.
1256  periodicVerts[mIt->first].push_back(ent2);
1257  perIt = periodicVerts.find(mIt->first);
1258  }
1259  else
1260  {
1261  // Vertex is known - loop over the vertices
1262  // inside the record and potentially add vertex
1263  // mIt->second to the list.
1264  for (k = 0; k < perIt->second.size(); ++k)
1265  {
1266  if (perIt->second[k].id == mIt->second.first)
1267  {
1268  break;
1269  }
1270  }
1271 
1272  if (k == perIt->second.size())
1273  {
1274  perIt->second.push_back(ent2);
1275  }
1276  }
1277  }
1278  }
1279 
1280  // Determine periodic edges. Logic is the same as above, and
1281  // perhaps should be condensed to avoid replication.
1282  for (i = 0; i < 2; ++i)
1283  {
1284  int other = (i+1) % 2;
1285 
1286  if (tmpVec[0].size() == 3)
1287  {
1289  tmpVec[i], tmpVec[other]);
1290  }
1291  else
1292  {
1294  tmpVec[i], tmpVec[other]);
1295  }
1296 
1297  vector<int> per1 = edgeMap[ids[i]];
1298  vector<int> per2 = edgeMap[ids[other]];
1299 
1300  map<int, pair<int, bool> > tmpMap;
1301  map<int, pair<int, bool> >::iterator mIt;
1302 
1303  for (j = 0; j < nFaceVerts; ++j)
1304  {
1305  int e = emap[nFaceVerts][o][j];
1306  tmpMap[per1[j]] = make_pair(
1307  per2[e], locEdges.count(per2[e]) > 0);
1308  }
1309 
1310  for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1311  {
1312  // Note we assume orientation of edges is forwards -
1313  // this may be reversed later.
1314  PeriodicEntity ent2(mIt->second.first,
1316  mIt->second.second);
1317  PeriodicMap::iterator perIt = periodicEdges.find(
1318  mIt->first);
1319 
1320  if (perIt == periodicEdges.end())
1321  {
1322  periodicEdges[mIt->first].push_back(ent2);
1323  perIt = periodicEdges.find(mIt->first);
1324  }
1325  else
1326  {
1327  for (k = 0; k < perIt->second.size(); ++k)
1328  {
1329  if (perIt->second[k].id == mIt->second.first)
1330  {
1331  break;
1332  }
1333  }
1334 
1335  if (k == perIt->second.size())
1336  {
1337  perIt->second.push_back(ent2);
1338  }
1339  }
1340  }
1341  }
1342  }
1343  }
1344 
1345  Array<OneD, int> pairSizes(n, 0);
1346  pairSizes[p] = allCompPairs.size();
1347  vComm->AllReduce(pairSizes, LibUtilities::ReduceSum);
1348 
1349  int totPairSizes = Vmath::Vsum(n, pairSizes, 1);
1350 
1351  Array<OneD, int> pairOffsets(n, 0);
1352  pairOffsets[0] = 0;
1353 
1354  for (i = 1; i < n; ++i)
1355  {
1356  pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1357  }
1358 
1359  Array<OneD, int> first (totPairSizes, 0);
1360  Array<OneD, int> second(totPairSizes, 0);
1361 
1362  cnt = pairOffsets[p];
1363 
1364  for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1365  {
1366  first [cnt ] = pIt->first;
1367  second[cnt++] = pIt->second;
1368  }
1369 
1370  vComm->AllReduce(first, LibUtilities::ReduceSum);
1371  vComm->AllReduce(second, LibUtilities::ReduceSum);
1372 
1373  allCompPairs.clear();
1374 
1375  for(cnt = 0; cnt < totPairSizes; ++cnt)
1376  {
1377  allCompPairs[first[cnt]] = second[cnt];
1378  }
1379 
1380  // Search for periodic vertices and edges which are not in
1381  // a periodic composite but lie in this process. First,
1382  // loop over all information we have from other
1383  // processors.
1384  for (cnt = i = 0; i < totFaces; ++i)
1385  {
1386  int faceId = faceIds[i];
1387 
1388  ASSERTL0(allCompPairs.count(faceId) > 0,
1389  "Unable to find matching periodic face.");
1390 
1391  int perFaceId = allCompPairs[faceId];
1392 
1393  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1394  {
1395  int vId = vertIds[cnt];
1396 
1397  PeriodicMap::iterator perId = periodicVerts.find(vId);
1398 
1399  if (perId == periodicVerts.end())
1400  {
1401 
1402  // This vertex is not included in the map. Figure out which
1403  // vertex it is supposed to be periodic with. perFaceId is
1404  // the face ID which is periodic with faceId. The logic is
1405  // much the same as the loop above.
1407  = { coordMap[faceId], coordMap[perFaceId] };
1408 
1409  int nFaceVerts = tmpVec[0].size();
1410  StdRegions::Orientation o = nFaceVerts == 3 ?
1412  tmpVec[0], tmpVec[1]) :
1414  tmpVec[0], tmpVec[1]);
1415 
1416  // Use vmap to determine which vertex of the other face
1417  // should be periodic with this one.
1418  int perVertexId = vertMap[perFaceId][vmap[nFaceVerts][o][j]];
1419 
1420 
1421  PeriodicEntity ent(perVertexId,
1423  locVerts.count(perVertexId) > 0);
1424 
1425  periodicVerts[vId].push_back(ent);
1426  }
1427 
1428  int eId = edgeIds[cnt];
1429 
1430  perId = periodicEdges.find(eId);
1431 
1432  if (perId == periodicEdges.end())
1433  {
1434  // This edge is not included in the map. Figure
1435  // out which edge it is supposed to be periodic
1436  // with. perFaceId is the face ID which is
1437  // periodic with faceId. The logic is much the
1438  // same as the loop above.
1440  = { coordMap[faceId], coordMap[perFaceId] };
1441 
1442  int nFaceEdges = tmpVec[0].size();
1443  StdRegions::Orientation o = nFaceEdges == 3 ?
1445  tmpVec[0], tmpVec[1]) :
1447  tmpVec[0], tmpVec[1]);
1448 
1449  // Use emap to determine which edge of the other
1450  // face should be periodic with this one.
1451  int perEdgeId = edgeMap[perFaceId][emap[nFaceEdges][o][j]];
1452 
1453 
1454  PeriodicEntity ent(perEdgeId,
1456  locEdges.count(perEdgeId) > 0);
1457 
1458  periodicEdges[eId].push_back(ent);
1459  }
1460  }
1461  }
1462 
1463  // Finally, we must loop over the periodicVerts and periodicEdges
1464  // map to complete connectivity information.
1465  PeriodicMap::iterator perIt, perIt2;
1466  for (perIt = periodicVerts.begin();
1467  perIt != periodicVerts.end(); ++perIt)
1468  {
1469  // For each vertex that is periodic with this one...
1470  for (i = 0; i < perIt->second.size(); ++i)
1471  {
1472  // Find the vertex in the periodicVerts map...
1473  perIt2 = periodicVerts.find(perIt->second[i].id);
1474  ASSERTL0(perIt2 != periodicVerts.end(),
1475  "Couldn't find periodic vertex.");
1476 
1477  // Now search through this vertex's list and make sure that
1478  // we have a record of any vertices which aren't in the
1479  // original list.
1480  for (j = 0; j < perIt2->second.size(); ++j)
1481  {
1482  if (perIt2->second[j].id == perIt->first)
1483  {
1484  continue;
1485  }
1486 
1487  for (k = 0; k < perIt->second.size(); ++k)
1488  {
1489  if (perIt2->second[j].id == perIt->second[k].id)
1490  {
1491  break;
1492  }
1493  }
1494 
1495  if (k == perIt->second.size())
1496  {
1497  perIt->second.push_back(perIt2->second[j]);
1498  }
1499  }
1500  }
1501  }
1502 
1503  for (perIt = periodicEdges.begin();
1504  perIt != periodicEdges.end(); ++perIt)
1505  {
1506  for (i = 0; i < perIt->second.size(); ++i)
1507  {
1508  perIt2 = periodicEdges.find(perIt->second[i].id);
1509  ASSERTL0(perIt2 != periodicEdges.end(),
1510  "Couldn't find periodic edge.");
1511 
1512  for (j = 0; j < perIt2->second.size(); ++j)
1513  {
1514  if (perIt2->second[j].id == perIt->first)
1515  {
1516  continue;
1517  }
1518 
1519  for (k = 0; k < perIt->second.size(); ++k)
1520  {
1521  if (perIt2->second[j].id == perIt->second[k].id)
1522  {
1523  break;
1524  }
1525  }
1526 
1527  if (k == perIt->second.size())
1528  {
1529  perIt->second.push_back(perIt2->second[j]);
1530  }
1531  }
1532  }
1533  }
1534 
1535  // Loop over periodic edges to determine relative edge orientations.
1536  for (perIt = periodicEdges.begin();
1537  perIt != periodicEdges.end(); perIt++)
1538  {
1539  // Find edge coordinates
1540  map<int, pair<int, int> >::iterator eIt
1541  = eIdMap.find(perIt->first);
1542  SpatialDomains::PointGeom v[2] = {
1543  *vCoMap[eIt->second.first],
1544  *vCoMap[eIt->second.second]
1545  };
1546 
1547  // Loop over each edge, and construct a vector that takes us
1548  // from one vertex to another. Use this to figure out which
1549  // vertex maps to which.
1550  for (i = 0; i < perIt->second.size(); ++i)
1551  {
1552  eIt = eIdMap.find(perIt->second[i].id);
1553 
1554  SpatialDomains::PointGeom w[2] = {
1555  *vCoMap[eIt->second.first],
1556  *vCoMap[eIt->second.second]
1557  };
1558 
1559  NekDouble cx = 0.5*(w[0](0)-v[0](0)+w[1](0)-v[1](0));
1560  NekDouble cy = 0.5*(w[0](1)-v[0](1)+w[1](1)-v[1](1));
1561  NekDouble cz = 0.5*(w[0](2)-v[0](2)+w[1](2)-v[1](2));
1562 
1563  int vMap[2] = {-1,-1};
1564  for (j = 0; j < 2; ++j)
1565  {
1566  NekDouble x = v[j](0);
1567  NekDouble y = v[j](1);
1568  NekDouble z = v[j](2);
1569  for (k = 0; k < 2; ++k)
1570  {
1571  NekDouble x1 = w[k](0)-cx;
1572  NekDouble y1 = w[k](1)-cy;
1573  NekDouble z1 = w[k](2)-cz;
1574 
1575  if (sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z1-z)*(z1-z))
1576  < 1e-8)
1577  {
1578  vMap[k] = j;
1579  break;
1580  }
1581  }
1582  }
1583 
1584  // Sanity check the map.
1585  ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
1586  "Unable to align periodic vertices.");
1587  ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
1588  (vMap[1] == 0 || vMap[1] == 1) &&
1589  (vMap[0] != vMap[1]),
1590  "Unable to align periodic vertices.");
1591 
1592  // If 0 -> 0 then edges are aligned already; otherwise
1593  // reverse the orientation.
1594  if (vMap[0] != 0)
1595  {
1596  perIt->second[i].orient = StdRegions::eBackwards;
1597  }
1598  }
1599  }
1600 
1601  // Do one final loop over periodic vertices/edges to remove
1602  // non-local vertices/edges from map.
1603  for (perIt = periodicVerts.begin();
1604  perIt != periodicVerts.end(); ++perIt)
1605  {
1606  if (locVerts.count(perIt->first) > 0)
1607  {
1608  m_periodicVerts.insert(*perIt);
1609  }
1610  }
1611 
1612  for (perIt = periodicEdges.begin();
1613  perIt != periodicEdges.end(); ++perIt)
1614  {
1615  if (locEdges.count(perIt->first) > 0)
1616  {
1617  m_periodicEdges.insert(*perIt);
1618  }
1619  }
1620  }
1621 
1622  bool DisContField3D::IsLeftAdjacentFace(const int n, const int e)
1623  {
1624  set<int>::iterator it;
1626  m_traceMap->GetElmtToTrace()[n][e]->
1627  as<LocalRegions::Expansion2D>();
1628 
1629  int offset = m_trace->GetPhys_Offset(traceEl->GetElmtId());
1630 
1631  bool fwd = true;
1632  if (traceEl->GetLeftAdjacentElementFace () == -1 ||
1633  traceEl->GetRightAdjacentElementFace() == -1)
1634  {
1635  // Boundary edge (1 connected element). Do nothing in
1636  // serial.
1637  it = m_boundaryFaces.find(traceEl->GetElmtId());
1638 
1639  // If the edge does not have a boundary condition set on
1640  // it, then assume it is a partition edge.
1641  if (it == m_boundaryFaces.end())
1642  {
1643  int traceGeomId = traceEl->GetGeom2D()->GetGlobalID();
1645  traceGeomId);
1646 
1647  if (pIt != m_periodicFaces.end() && !pIt->second[0].isLocal)
1648  {
1649  fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1650  }
1651  else
1652  {
1653  fwd = m_traceMap->
1654  GetTraceToUniversalMapUnique(offset) >= 0;
1655  }
1656  }
1657  }
1658  else if (traceEl->GetLeftAdjacentElementFace () != -1 &&
1659  traceEl->GetRightAdjacentElementFace() != -1)
1660  {
1661  // Non-boundary edge (2 connected elements).
1662  fwd = dynamic_cast<Nektar::StdRegions::StdExpansion*>
1663  (traceEl->GetLeftAdjacentElementExp().get()) ==
1664  (*m_exp)[n].get();
1665  }
1666  else
1667  {
1668  ASSERTL2(false, "Unconnected trace element!");
1669  }
1670 
1671  return fwd;
1672  }
1673 
1674  /**
1675  * \brief This method extracts the "forward" and "backward" trace
1676  * data from the array \a field and puts the data into output
1677  * vectors \a Fwd and \a Bwd.
1678  *
1679  * We first define the convention which defines "forwards" and
1680  * "backwards". First an association is made between the face of
1681  * each element and its corresponding face in the trace space
1682  * using the mapping #m_traceMap. The element can either be
1683  * left-adjacent or right-adjacent to this trace face (see
1684  * Expansion2D::GetLeftAdjacentElementExp). Boundary faces are
1685  * always left-adjacent since left-adjacency is populated first.
1686  *
1687  * If the element is left-adjacent we extract the face trace data
1688  * from \a field into the forward trace space \a Fwd; otherwise,
1689  * we place it in the backwards trace space \a Bwd. In this way,
1690  * we form a unique set of trace normals since these are always
1691  * extracted from left-adjacent elements.
1692  *
1693  * \param field is a NekDouble array which contains the 3D data
1694  * from which we wish to extract the backward and forward
1695  * orientated trace/face arrays.
1696  *
1697  * \return Updates a NekDouble array \a Fwd and \a Bwd
1698  */
1701  {
1702  v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
1703  }
1704 
1705 
1707  const Array<OneD, const NekDouble> &field,
1710  {
1711  // Loop over elements and collect forward and backward expansions.
1712  int nexp = GetExpSize();
1713  int cnt, n, e, npts, offset, phys_offset;
1714  Array<OneD,NekDouble> e_tmp;
1715 
1717  &elmtToTrace = m_traceMap->GetElmtToTrace();
1718 
1719  set<int>::iterator it;
1721  boost::unordered_map<int,pair<int,int> >::iterator it3;
1722 
1723  // Zero vectors.
1724  Vmath::Zero(Fwd.num_elements(), Fwd, 1);
1725  Vmath::Zero(Bwd.num_elements(), Bwd, 1);
1726 
1728  bool fwd;
1729 
1730  for(cnt = n = 0; n < nexp; ++n)
1731  {
1732  exp3d = (*m_exp)[n]->as<LocalRegions::Expansion3D>();
1733  phys_offset = GetPhys_Offset(n);
1734  for(e = 0; e < exp3d->GetNfaces(); ++e, ++cnt)
1735  {
1736  offset = m_trace->GetPhys_Offset(
1737  elmtToTrace[n][e]->GetElmtId());
1738 
1739  fwd = m_leftAdjacentFaces[cnt];
1740  if (fwd)
1741  {
1742  exp3d->GetFacePhysVals(e, elmtToTrace[n][e],
1743  field + phys_offset,
1744  e_tmp = Fwd + offset);
1745  }
1746  else
1747  {
1748  exp3d->GetFacePhysVals(e, elmtToTrace[n][e],
1749  field + phys_offset,
1750  e_tmp = Bwd + offset);
1751  }
1752  }
1753  }
1754 
1755  // fill boundary conditions into missing elements
1756  int id1,id2 = 0;
1757  cnt = 0;
1758 
1759  for(n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1760  {
1761  if(m_bndConditions[n]->GetBoundaryConditionType() ==
1763  {
1764  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
1765  {
1766  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
1767  id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
1768  id2 = m_trace->GetPhys_Offset(
1769  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1770  Vmath::Vcopy(npts,
1771  &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
1772  &Bwd[id2], 1);
1773  }
1774 
1775  cnt += e;
1776  }
1777  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
1779  m_bndConditions[n]->GetBoundaryConditionType() ==
1781  {
1782  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
1783  {
1784  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
1785  id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
1786  id2 = m_trace->GetPhys_Offset(
1787  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1788 
1789  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[id1] == 0.0,
1790  "method not set up for non-zero Neumann "
1791  "boundary condition");
1792 
1793  Vmath::Vcopy(npts,&Fwd[id2],1,&Bwd[id2],1);
1794  }
1795 
1796  cnt += e;
1797  }
1798  else
1799  {
1800  ASSERTL0(false, "Method only set up for Dirichlet, Neumann "
1801  "and Robin conditions.");
1802  }
1803  }
1804 
1805  // Copy any periodic boundary conditions.
1806  for (n = 0; n < m_periodicFwdCopy.size(); ++n)
1807  {
1808  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
1809  }
1810 
1811  // Do parallel exchange for forwards/backwards spaces.
1812  m_traceMap->UniversalTraceAssemble(Fwd);
1813  m_traceMap->UniversalTraceAssemble(Bwd);
1814  }
1815 
1817  Array<OneD, NekDouble> &outarray)
1818  {
1819  ASSERTL1(m_physState == true,
1820  "Field is not in physical space.");
1821 
1822  v_ExtractTracePhys(m_phys, outarray);
1823  }
1824 
1826  const Array<OneD, const NekDouble> &inarray,
1827  Array<OneD, NekDouble> &outarray)
1828  {
1829  // Loop over elemente and collect forward expansion
1830  int nexp = GetExpSize();
1831  int n,e,offset,phys_offset;
1832  Array<OneD,NekDouble> e_tmp;
1834  &elmtToTrace = m_traceMap->GetElmtToTrace();
1835 
1836  ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
1837  "input array is of insufficient length");
1838 
1839  // use m_trace tmp space in element to fill values
1840  for(n = 0; n < nexp; ++n)
1841  {
1842  phys_offset = GetPhys_Offset(n);
1843 
1844  for(e = 0; e < (*m_exp)[n]->GetNfaces(); ++e)
1845  {
1846  offset = m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
1847  (*m_exp)[n]->GetFacePhysVals(e, elmtToTrace[n][e],
1848  inarray + phys_offset,
1849  e_tmp = outarray + offset);
1850  }
1851  }
1852  }
1853 
1854  /**
1855  * @brief Add trace contributions into elemental coefficient spaces.
1856  *
1857  * Given some quantity \f$ \vec{Fn} \f$, which conatins this
1858  * routine calculates the integral
1859  *
1860  * \f[
1861  * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
1862  * \f]
1863  *
1864  * and adds this to the coefficient space provided by outarray.
1865  *
1866  * @see Expansion3D::AddFaceNormBoundaryInt
1867  *
1868  * @param Fn The trace quantities.
1869  * @param outarray Resulting 3D coefficient space.
1870  *
1871  */
1873  const Array<OneD, const NekDouble> &Fn,
1874  Array<OneD, NekDouble> &outarray)
1875  {
1876  int e,n,offset, t_offset;
1877  Array<OneD, NekDouble> e_outarray;
1879  &elmtToTrace = m_traceMap->GetElmtToTrace();
1880 
1881  for(n = 0; n < GetExpSize(); ++n)
1882  {
1883  offset = GetCoeff_Offset(n);
1884  e_outarray = outarray+offset;
1885  for(e = 0; e < (*m_exp)[n]->GetNfaces(); ++e)
1886  {
1887  t_offset = m_trace->GetPhys_Offset(
1888  elmtToTrace[n][e]->GetElmtId());
1889  (*m_exp)[n]->AddFaceNormBoundaryInt(e,
1890  elmtToTrace[n][e],
1891  Fn + t_offset,
1892  e_outarray);
1893  }
1894  }
1895  }
1896  /**
1897  * @brief Add trace contributions into elemental coefficient spaces.
1898  *
1899  * Given some quantity \f$ \vec{Fn} \f$, which conatins this
1900  * routine calculates the integral
1901  *
1902  * \f[
1903  * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
1904  * \f]
1905  *
1906  * and adds this to the coefficient space provided by
1907  * outarray. The value of q is determined from the routine
1908  * IsLeftAdjacentFace() which if true we use Fwd else we use
1909  * Bwd
1910  *
1911  * @see Expansion3D::AddFaceNormBoundaryInt
1912  *
1913  * @param Fwd The trace quantities associated with left (fwd)
1914  * adjancent elmt.
1915  * @param Bwd The trace quantities associated with right (bwd)
1916  * adjacent elet.
1917  * @param outarray Resulting 3D coefficient space.
1918  */
1920  const Array<OneD, const NekDouble> &Fwd,
1921  const Array<OneD, const NekDouble> &Bwd,
1922  Array<OneD, NekDouble> &outarray)
1923  {
1924  int e,n,offset, t_offset;
1925  Array<OneD, NekDouble> e_outarray;
1927  &elmtToTrace = m_traceMap->GetElmtToTrace();
1928 
1929  for(n = 0; n < GetExpSize(); ++n)
1930  {
1931  offset = GetCoeff_Offset(n);
1932  for(e = 0; e < (*m_exp)[n]->GetNfaces(); ++e)
1933  {
1934  t_offset = m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
1935 
1936  // Evaluate upwind flux less local edge
1937  if(IsLeftAdjacentFace(n,e))
1938  {
1939  (*m_exp)[n]->AddFaceNormBoundaryInt(
1940  e, elmtToTrace[n][e], Fwd + t_offset,
1941  e_outarray = outarray+offset);
1942  }
1943  else
1944  {
1945  (*m_exp)[n]->AddFaceNormBoundaryInt(
1946  e, elmtToTrace[n][e], Bwd + t_offset,
1947  e_outarray = outarray+offset);
1948  }
1949  }
1950  }
1951  }
1952 
1953  /**
1954  * @brief Set up a list of elemeent IDs and edge IDs that link to the
1955  * boundary conditions.
1956  */
1958  Array<OneD, int> &ElmtID,
1959  Array<OneD, int> &FaceID)
1960  {
1961  map<int,int> globalIdMap;
1962  int i, n;
1963  int cnt;
1964  int nbcs = 0;
1965 
1967  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph3D>(
1968  m_graph);
1969 
1970  // Populate global ID map (takes global geometry ID to local
1971  // expansion list ID).
1973  for (i = 0; i < GetExpSize(); ++i)
1974  {
1975  exp3d = (*m_exp)[i]->as<LocalRegions::Expansion3D>();
1976  globalIdMap[exp3d->GetGeom3D()->GetGlobalID()] = i;
1977  }
1978 
1979  // Determine number of boundary condition expansions.
1980  for(i = 0; i < m_bndConditions.num_elements(); ++i)
1981  {
1982  nbcs += m_bndCondExpansions[i]->GetExpSize();
1983  }
1984 
1985  // make sure arrays are of sufficient length
1986  if(ElmtID.num_elements() != nbcs)
1987  {
1988  ElmtID = Array<OneD, int>(nbcs);
1989  }
1990 
1991  if(FaceID.num_elements() != nbcs)
1992  {
1993  FaceID = Array<OneD, int>(nbcs);
1994  }
1995 
1997  for(cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1998  {
1999  for(i = 0; i < m_bndCondExpansions[n]->GetExpSize(); ++i, ++cnt)
2000  {
2001  exp2d = m_bndCondExpansions[n]->GetExp(i)->
2002  as<LocalRegions::Expansion2D>();
2003  // Use face to element map from MeshGraph3D.
2005  graph3D->GetElementsFromFace(exp2d->GetGeom2D());
2006 
2007  ElmtID[cnt] = globalIdMap[(*tmp)[0]->
2008  m_Element->GetGlobalID()];
2009  FaceID[cnt] = (*tmp)[0]->m_FaceIndx;
2010  }
2011  }
2012  }
2013 
2014  /**
2015  * @brief Reset this field, so that geometry information can be updated.
2016  */
2018  {
2019  ExpList::v_Reset();
2020 
2021  // Reset boundary condition expansions.
2022  for (int n = 0; n < m_bndCondExpansions.num_elements(); ++n)
2023  {
2024  m_bndCondExpansions[n]->Reset();
2025  }
2026  }
2027 
2028  /**
2029  * Solving Helmholtz Equation in 3D
2030  */
2032  const Array<OneD, const NekDouble> &inarray,
2033  Array<OneD, NekDouble> &outarray,
2034  const FlagList &flags,
2035  const StdRegions::ConstFactorMap &factors,
2036  const StdRegions::VarCoeffMap &varcoeff,
2037  const Array<OneD, const NekDouble> &dirForcing)
2038  {
2039  int i,j,n,cnt,cnt1,nbndry;
2040  int nexp = GetExpSize();
2042 
2044  DNekVec F(m_ncoeffs,f,eWrapper);
2045  Array<OneD,NekDouble> e_f, e_l;
2046 
2047  //----------------------------------
2048  // Setup RHS Inner product
2049  //----------------------------------
2050  IProductWRTBase(inarray,f);
2051  Vmath::Neg(m_ncoeffs,f,1);
2052 
2053  //----------------------------------
2054  // Solve continuous flux System
2055  //----------------------------------
2056  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
2057  int NumDirichlet = m_traceMap->GetNumLocalDirBndCoeffs();
2058  int e_ncoeffs,id;
2059 
2060  // Retrieve block matrix of U^e
2062  const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
2063 
2064  // Retrieve global trace space storage, \Lambda, from trace expansion
2065  Array<OneD,NekDouble> BndSol = m_trace->UpdateCoeffs();
2066 
2067  // Create trace space forcing, F
2068  Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
2069 
2070  // Zero \Lambda
2071  Vmath::Zero(GloBndDofs,BndSol,1);
2072 
2073  // Retrieve number of local trace space coefficients N_{\lambda},
2074  // and set up local elemental trace solution \lambda^e.
2075  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
2076  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
2077  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
2078 
2079  //----------------------------------
2080  // Evaluate Trace Forcing vector F
2081  // Kirby et al, 2010, P23, Step 5.
2082  //----------------------------------
2083  // Loop over all expansions in the domain
2084  for(cnt = cnt1 = n = 0; n < nexp; ++n)
2085  {
2086  nbndry = (*m_exp)[m_offset_elmt_id[n]]->NumDGBndryCoeffs();
2087 
2088  e_ncoeffs = (*m_exp)[m_offset_elmt_id[n]]->GetNcoeffs();
2089  e_f = f + cnt;
2090  e_l = loc_lambda + cnt1;
2091 
2092  // Local trace space \lambda^e
2093  DNekVec Floc (nbndry, e_l, eWrapper);
2094  // Local forcing f^e
2095  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
2096  // Compute local (U^e)^{\top} f^e
2097  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
2098 
2099  cnt += e_ncoeffs;
2100  cnt1 += nbndry;
2101  }
2102 
2103  // Assemble local \lambda_e into global \Lambda
2104  m_traceMap->AssembleBnd(loc_lambda,BndRhs);
2105 
2106  // Copy Dirichlet boundary conditions and weak forcing into trace
2107  // space
2108  cnt = 0;
2109  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
2110  {
2111  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eDirichlet)
2112  {
2113  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
2114  {
2115  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2116  BndSol[id] = m_bndCondExpansions[i]->GetCoeffs()[j];
2117  }
2118  }
2119  else
2120  {
2121  //Add weak boundary condition to trace forcing
2122  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
2123  {
2124  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
2125  BndRhs[id] += m_bndCondExpansions[i]->GetCoeffs()[j];
2126  }
2127  }
2128  }
2129 
2130  //----------------------------------
2131  // Solve trace problem: \Lambda = K^{-1} F
2132  // K is the HybridDGHelmBndLam matrix.
2133  //----------------------------------
2134  if(GloBndDofs - NumDirichlet > 0)
2135  {
2137  m_traceMap,factors,varcoeff);
2139  LinSys->Solve(BndRhs,BndSol,m_traceMap);
2140  }
2141 
2142  //----------------------------------
2143  // Internal element solves
2144  //----------------------------------
2146  const DNekScalBlkMatSharedPtr& InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
2147  DNekVec out(m_ncoeffs,outarray,eWrapper);
2148  Vmath::Zero(m_ncoeffs,outarray,1);
2149 
2150  // get local trace solution from BndSol
2151  m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
2152 
2153  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
2154  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
2155  }
2156 
2157  /**
2158  * @brief Calculates the result of the multiplication of a global matrix
2159  * of type specified by @a mkey with a vector given by @a inarray.
2160  *
2161  * @param mkey Key representing desired matrix multiplication.
2162  * @param inarray Input vector.
2163  * @param outarray Resulting multiplication.
2164  */
2166  const GlobalMatrixKey &gkey,
2167  const Array<OneD,const NekDouble> &inarray,
2168  Array<OneD, NekDouble> &outarray,
2169  CoeffState coeffstate)
2170  {
2171  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
2172  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
2173  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
2174  const DNekScalBlkMatSharedPtr& HDGHelm = GetBlockMatrix(gkey);
2175 
2176  m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2177  LocLambda = (*HDGHelm) * LocLambda;
2178  m_traceMap->AssembleBnd(loc_lambda,outarray);
2179  }
2180 
2181  /**
2182  * Search through the edge expansions and identify which ones have
2183  * Robin/Mixed type boundary conditions. If find a Robin boundary then
2184  * store the edge id of the boundary condition and the array of points
2185  * of the physical space boundary condition which are hold the boundary
2186  * condition primitive variable coefficient at the quatrature points
2187  *
2188  * \return std map containing the robin boundary condition
2189  * info using a key of the element id
2190  *
2191  * There is a next member to allow for more than one Robin
2192  * boundary condition per element
2193  */
2194  map<int, RobinBCInfoSharedPtr> DisContField3D::v_GetRobinBCInfo(void)
2195  {
2196  int i,cnt;
2197  map<int, RobinBCInfoSharedPtr> returnval;
2198  Array<OneD, int> ElmtID,FaceID;
2199  GetBoundaryToElmtMap(ElmtID,FaceID);
2200 
2201  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
2202  {
2203  MultiRegions::ExpListSharedPtr locExpList;
2204 
2205  if(m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::eRobin)
2206  {
2207  int e,elmtid;
2208  Array<OneD, NekDouble> Array_tmp;
2209 
2210  locExpList = m_bndCondExpansions[i];
2211 
2212  for(e = 0; e < locExpList->GetExpSize(); ++e)
2213  {
2214  RobinBCInfoSharedPtr rInfo = MemoryManager<RobinBCInfo>::AllocateSharedPtr(FaceID[cnt+e],Array_tmp = locExpList->GetPhys() + locExpList->GetPhys_Offset(e));
2215  elmtid = ElmtID[cnt+e];
2216  // make link list if necessary
2217  if(returnval.count(elmtid) != 0)
2218  {
2219  rInfo->next = returnval.find(elmtid)->second;
2220  }
2221  returnval[elmtid] = rInfo;
2222  }
2223  }
2224  cnt += m_bndCondExpansions[i]->GetExpSize();
2225  }
2226 
2227  return returnval;
2228  }
2229 
2230  /**
2231  * @brief Evaluate HDG post-processing to increase polynomial order of
2232  * solution.
2233  *
2234  * This function takes the solution (assumed to be one order lower) in
2235  * physical space, and postprocesses at the current polynomial order by
2236  * solving the system:
2237  *
2238  * \f[
2239  * \begin{aligned}
2240  * (\nabla w, \nabla u^*) &= (\nabla w, u), \\
2241  * \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
2242  * \end{aligned}
2243  * \f]
2244  *
2245  * where \f$ u \f$ corresponds with the current solution as stored
2246  * inside #m_coeffs.
2247  *
2248  * @param outarray The resulting field \f$ u^* \f$.
2249  */
2251  Array<OneD, NekDouble> &outarray)
2252  {
2253  int i,cnt,f,ncoeff_face;
2254  Array<OneD, NekDouble> force, out_tmp,qrhs,qrhs1;
2256  &elmtToTrace = m_traceMap->GetElmtToTrace();
2257 
2258  int eid,nq_elmt, nm_elmt;
2259  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
2260  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), face_lambda;
2261  Array<OneD, NekDouble> tmp_coeffs;
2262  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
2263 
2264  face_lambda = loc_lambda;
2265 
2266  // Calculate Q using standard DG formulation.
2267  for(i = cnt = 0; i < GetExpSize(); ++i)
2268  {
2270  (*m_exp)[i]->as<LocalRegions::Expansion3D>();
2271 
2272  eid = m_offset_elmt_id[i];
2273  nq_elmt = (*m_exp)[eid]->GetTotPoints();
2274  nm_elmt = (*m_exp)[eid]->GetNcoeffs();
2275  qrhs = Array<OneD, NekDouble>(nq_elmt);
2276  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
2277  force = Array<OneD, NekDouble>(2*nm_elmt);
2278  out_tmp = force + nm_elmt;
2280 
2281  int num_points0 = (*m_exp)[eid]->GetBasis(0)->GetNumPoints();
2282  int num_points1 = (*m_exp)[eid]->GetBasis(1)->GetNumPoints();
2283  int num_points2 = (*m_exp)[eid]->GetBasis(2)->GetNumPoints();
2284  int num_modes0 = (*m_exp)[eid]->GetBasis(0)->GetNumModes();
2285  int num_modes1 = (*m_exp)[eid]->GetBasis(1)->GetNumModes();
2286  int num_modes2 = (*m_exp)[eid]->GetBasis(2)->GetNumModes();
2287 
2288  // Probably a better way of setting up lambda than this. Note
2289  // cannot use PutCoeffsInToElmts since lambda space is mapped
2290  // during the solve.
2291  int nFaces = (*m_exp)[eid]->GetNfaces();
2292  Array<OneD, Array<OneD, NekDouble> > faceCoeffs(nFaces);
2293  for(f = 0; f < nFaces; ++f)
2294  {
2295  ncoeff_face = elmtToTrace[eid][f]->GetNcoeffs();
2296  faceCoeffs[f] = Array<OneD, NekDouble>(ncoeff_face);
2297  Vmath::Vcopy(ncoeff_face, face_lambda, 1, faceCoeffs[f], 1);
2298  exp->SetFaceToGeomOrientation(f, faceCoeffs[f]);
2299  face_lambda = face_lambda + ncoeff_face;
2300  }
2301 
2302  //creating orthogonal expansion (checking if we have quads or triangles)
2303  LibUtilities::ShapeType shape = (*m_exp)[eid]->DetShapeType();
2304  switch(shape)
2305  {
2307  {
2311  LibUtilities::BasisKey BkeyH1(LibUtilities::eOrtho_A, num_modes0, PkeyH1);
2312  LibUtilities::BasisKey BkeyH2(LibUtilities::eOrtho_A, num_modes1, PkeyH2);
2313  LibUtilities::BasisKey BkeyH3(LibUtilities::eOrtho_A, num_modes2, PkeyH3);
2314  SpatialDomains::HexGeomSharedPtr hGeom = boost::dynamic_pointer_cast<SpatialDomains::HexGeom>((*m_exp)[eid]->GetGeom());
2315  ppExp = MemoryManager<LocalRegions::HexExp>::AllocateSharedPtr(BkeyH1, BkeyH2, BkeyH3, hGeom);
2316  }
2317  break;
2319  {
2323  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1);
2324  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes1, PkeyT2);
2325  LibUtilities::BasisKey BkeyT3(LibUtilities::eOrtho_C, num_modes2, PkeyT3);
2326  SpatialDomains::TetGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TetGeom>((*m_exp)[eid]->GetGeom());
2327  ppExp = MemoryManager<LocalRegions::TetExp>::AllocateSharedPtr(BkeyT1, BkeyT2, BkeyT3, tGeom);
2328  }
2329  break;
2330  case LibUtilities::ePrism:
2331  {
2335  LibUtilities::BasisKey BkeyP1(LibUtilities::eOrtho_A, num_modes0, PkeyP1);
2336  LibUtilities::BasisKey BkeyP2(LibUtilities::eOrtho_A, num_modes1, PkeyP2);
2337  LibUtilities::BasisKey BkeyP3(LibUtilities::eOrtho_B, num_modes2, PkeyP3);
2338  SpatialDomains::PrismGeomSharedPtr pGeom = boost::dynamic_pointer_cast<SpatialDomains::PrismGeom>((*m_exp)[eid]->GetGeom());
2339  ppExp = MemoryManager<LocalRegions::PrismExp>::AllocateSharedPtr(BkeyP1, BkeyP2, BkeyP3, pGeom);
2340  }
2341  break;
2342  default:
2343  ASSERTL0(false, "Wrong shape type, HDG postprocessing is not implemented");
2344  };
2345 
2346 
2347  //DGDeriv
2348  // (d/dx w, q_0)
2349  (*m_exp)[eid]->DGDeriv(
2350  0,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
2351  elmtToTrace[eid], faceCoeffs, out_tmp);
2352  (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2353  ppExp->IProductWRTDerivBase(0,qrhs,force);
2354 
2355 
2356  // + (d/dy w, q_1)
2357  (*m_exp)[eid]->DGDeriv(
2358  1,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
2359  elmtToTrace[eid], faceCoeffs, out_tmp);
2360  (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2361  ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2362 
2363  Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
2364 
2365  // + (d/dz w, q_2)
2366  (*m_exp)[eid]->DGDeriv(
2367  2,tmp_coeffs = m_coeffs + m_coeff_offset[eid],
2368  elmtToTrace[eid], faceCoeffs, out_tmp);
2369  (*m_exp)[eid]->BwdTrans(out_tmp,qrhs);
2370  ppExp->IProductWRTDerivBase(2,qrhs,out_tmp);
2371 
2372  Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
2373  // determine force[0] = (1,u)
2374  (*m_exp)[eid]->BwdTrans(
2375  tmp_coeffs = m_coeffs + m_coeff_offset[eid],qrhs);
2376  force[0] = (*m_exp)[eid]->Integral(qrhs);
2377 
2378  // multiply by inverse Laplacian matrix
2379  // get matrix inverse
2380  LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean, ppExp->DetShapeType(), *ppExp);
2381  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
2382 
2383  NekVector<NekDouble> in (nm_elmt, force, eWrapper);
2384  NekVector<NekDouble> out(nm_elmt);
2385 
2386  out = (*lapsys)*in;
2387 
2388  // Transforming back to modified basis
2389  Array<OneD, NekDouble> work(nq_elmt);
2390  ppExp->BwdTrans(out.GetPtr(), work);
2391  (*m_exp)[eid]->FwdTrans(work,
2392  tmp_coeffs = outarray + m_coeff_offset[eid]);
2393  }
2394  }
2395 
2396  /**
2397  * \brief This function evaluates the boundary conditions at a certain
2398  * time-level.
2399  *
2400  * Based on the boundary condition \f$g(\boldsymbol{x},t)\f$ evaluated
2401  * at a given time-level \a t, this function transforms the boundary
2402  * conditions onto the coefficients of the (one-dimensional) boundary
2403  * expansion. Depending on the type of boundary conditions, these
2404  * expansion coefficients are calculated in different ways:
2405  * - <b>Dirichlet boundary conditions</b><BR>
2406  * In order to ensure global \f$C^0\f$ continuity of the spectral/hp
2407  * approximation, the Dirichlet boundary conditions are projected onto
2408  * the boundary expansion by means of a modified \f$C^0\f$ continuous
2409  * Galerkin projection. This projection can be viewed as a collocation
2410  * projection at the vertices, followed by an \f$L^2\f$ projection on
2411  * the interior modes of the edges. The resulting coefficients
2412  * \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ will be stored for the
2413  * boundary expansion.
2414  * - <b>Neumann boundary conditions</b>
2415  * In the discrete Galerkin formulation of the problem to be solved,
2416  * the Neumann boundary conditions appear as the set of surface
2417  * integrals: \f[\boldsymbol{\hat{g}}=\int_{\Gamma}
2418  * \phi^e_n(\boldsymbol{x})g(\boldsymbol{x})d(\boldsymbol{x})\quad
2419  * \forall n \f]
2420  * As a result, it are the coefficients \f$\boldsymbol{\hat{g}}\f$
2421  * that will be stored in the boundary expansion
2422  *
2423  * @param time The time at which the boundary conditions
2424  * should be evaluated.
2425  * @param bndCondExpansions List of boundary conditions.
2426  * @param bndConditions Information about the boundary conditions.
2427  */
2429  const NekDouble time,
2430  const std::string varName,
2431  const NekDouble x2_in,
2432  const NekDouble x3_in)
2433  {
2434  int i;
2435  int npoints;
2436  int nbnd = m_bndCondExpansions.num_elements();
2437  MultiRegions::ExpListSharedPtr locExpList;
2438 
2439  for (i = 0; i < nbnd; ++i)
2440  {
2441  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
2442  {
2443  locExpList = m_bndCondExpansions[i];
2444  npoints = locExpList->GetNpoints();
2445 
2446  Array<OneD, NekDouble> x0(npoints, 0.0);
2447  Array<OneD, NekDouble> x1(npoints, 0.0);
2448  Array<OneD, NekDouble> x2(npoints, 0.0);
2449 
2450  locExpList->GetCoords(x0, x1, x2);
2451 
2452  if (m_bndConditions[i]->GetBoundaryConditionType()
2454  {
2455  string filebcs = boost::static_pointer_cast<
2457  m_bndConditions[i])->m_filename;
2458 
2459  if (filebcs != "")
2460  {
2461  ExtractFileBCs(filebcs, varName, locExpList);
2462  }
2463  else
2464  {
2465  LibUtilities::Equation condition = boost::static_pointer_cast<SpatialDomains::
2467  m_bndConditions[i])->m_dirichletCondition;
2468 
2469  condition.Evaluate(x0, x1, x2, time,
2470  locExpList->UpdatePhys());
2471 
2472  locExpList->FwdTrans_BndConstrained(
2473  locExpList->GetPhys(),
2474  locExpList->UpdateCoeffs());
2475  }
2476  }
2477  else if (m_bndConditions[i]->GetBoundaryConditionType()
2479  {
2480  LibUtilities::Equation condition = boost::
2481  static_pointer_cast<SpatialDomains::
2483  m_bndConditions[i])->m_neumannCondition;
2484 
2485  condition.Evaluate(x0, x1, x2, time,
2486  locExpList->UpdatePhys());
2487 
2488  locExpList->IProductWRTBase(locExpList->GetPhys(),
2489  locExpList->UpdateCoeffs());
2490  }
2491  else if (m_bndConditions[i]->GetBoundaryConditionType()
2493  {
2494  LibUtilities::Equation condition = boost::
2495  static_pointer_cast<SpatialDomains::
2497  m_bndConditions[i])->m_robinFunction;
2498 
2499  LibUtilities::Equation coeff = boost::
2500  static_pointer_cast<SpatialDomains::
2502  m_bndConditions[i])->m_robinPrimitiveCoeff;
2503 
2504  condition.Evaluate(x0, x1, x2, time,
2505  locExpList->UpdatePhys());
2506 
2507  locExpList->IProductWRTBase(locExpList->GetPhys(),
2508  locExpList->UpdateCoeffs());
2509 
2510  // Put primitive coefficient into the physical
2511  // space storage
2512  coeff.Evaluate(x0, x1, x2, time,
2513  locExpList->UpdatePhys());
2514 
2515  }
2516  else
2517  {
2518  ASSERTL0(false, "This type of BC not implemented yet");
2519  }
2520  }
2521  }
2522  }
2523  } // end of namespace
2524 } // end of namespace
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:772
boost::shared_ptr< MeshGraph3D > MeshGraph3DSharedPtr
Definition: MeshGraph3D.h:224
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1394
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
boost::shared_ptr< ElementFaceVector > ElementFaceVectorSharedPtr
Definition: MeshGraph.h:135
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:1867
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1343
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2027
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
std::map< int, vector< PeriodicEntity > > PeriodicMap
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry3D.h:61
void SetAdjacentElementExp(int face, Expansion3DSharedPtr &f)
Definition: Expansion2D.h:215
bool isLocal
Flag specifying if this entity is local to this partition.
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
Definition: ExpList.h:1875
bool IsLeftAdjacentFace(const int n, const int e)
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:235
void ExtractFileBCs(const std::string &fileName, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:1862
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
SpatialDomains::Geometry3DSharedPtr GetGeom3D() const
Definition: Expansion3D.h:149
std::set< int > m_boundaryFaces
A set storing the global IDs of any boundary faces.
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshPartition.h:53
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2080
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:926
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshPartition.h:52
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:1381
boost::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:110
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:909
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:1837
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
This function evaluates the boundary conditions at a certain time-level.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
boost::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:1152
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:57
Principle Orthogonal Functions .
Definition: BasisType.h:47
void FindPeriodicFaces(const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
Determine the periodic faces, edges and vertices for the given graph.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:958
virtual map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2)
Get the orientation of face1.
Definition: QuadGeom.cpp:257
The base class for all shapes.
Definition: StdExpansion.h:69
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:947
DisContField3D()
Default constructor.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:225
boost::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:48
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:935
NekDouble Evaluate() const
Definition: Equation.h:102
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:887
Array< OneD, int > m_offset_elmt_id
Array containing the element id m_offset_elmt_id[n] that the n^th consecutive block of data in m_coef...
Definition: ExpList.h:969
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:204
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
static std::string npts
Definition: InputFld.cpp:43
Principle Orthogonal Functions .
Definition: BasisType.h:48
Principle Orthogonal Functions .
Definition: BasisType.h:46
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:883
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph3D, const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Calculates the result of the multiplication of a global matrix of type specified by mkey with a vecto...
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:880
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
bool SameTypeOfBoundaryConditions(const DisContField3D &In)
Defines a list of flags.
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:213
boost::shared_ptr< ExpList2D > ExpList2DSharedPtr
Shared pointer to an ExpList2D object.
Definition: ExpList2D.h:49
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:111
Describe a linear system.
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:112
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:1828
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:877
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:55
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
2D geometry information
Definition: Geometry2D.h:65
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const Array< OneD, const NekDouble > &dirForcing)
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &FaceID)
Set up a list of elemeent IDs and edge IDs that link to the boundary conditions.
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1507
boost::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:109
vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:58
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1351
boost::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:106
Abstraction of a three-dimensional multi-elemental expansion which is merely a collection of local ex...
Definition: ExpList3D.h:49
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:2631
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2)
Definition: TriGeom.cpp:236
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:225
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:206
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:714
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
Describes the specification for a Basis.
Definition: Basis.h:50
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:50
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:285
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
PeriodicMap m_periodicEdges
A map which identifies groups of periodic edges.
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49