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