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