Nektar++
DisContField.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File DisContField.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 1D domain with boundary
33 // conditions using LDG-H
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
40 #include <LocalRegions/HexExp.h>
41 #include <LocalRegions/QuadExp.h>
42 #include <LocalRegions/TetExp.h>
43 #include <LocalRegions/TriExp.h>
46 #include <StdRegions/StdSegExp.h>
47 
48 using namespace std;
49 
50 namespace Nektar
51 {
52 namespace MultiRegions
53 {
54 /**
55  * @class DisContField
56  * This class augments the list of local expansions inherited from
57  * ExpList with boundary conditions. Inter-element boundaries are
58  * handled using an discontinuous Galerkin scheme.
59  */
60 
61 /**
62  * Constructs an empty expansion list with no boundary conditions.
63  */
64 DisContField::DisContField()
65  : ExpList(), m_bndCondExpansions(), m_bndConditions(),
66  m_trace(NullExpListSharedPtr)
67 {
68 }
69 
70 /**
71  * An expansion list for the boundary expansions is generated first for
72  * the field. These are subsequently evaluated for time zero. The trace
73  * map is then constructed.
74  * @param graph1D A mesh, containing information about the domain
75  * and the spectral/hp element expansions.
76  * @param bcs Information about the enforced boundary
77  * conditions.
78  * @param variable The session variable associated with the
79  * boundary conditions to enforce.
80  * @param solnType Type of global system to use.
81  */
84  const std::string &variable, const bool SetUpJustDG,
85  const bool DeclareCoeffPhysArrays,
86  const Collections::ImplementationType ImpType,
87  const std::string bcvariable)
88  : ExpList(pSession, graph, DeclareCoeffPhysArrays, variable, ImpType),
89  m_trace(NullExpListSharedPtr)
90 {
91  std::string bcvar;
92  if (bcvariable == "NotSet")
93  {
94  bcvar = variable;
95  }
96  else
97  {
98  bcvar = bcvariable;
99  }
100 
101  if (bcvar.compare("DefaultVar") != 0)
102  {
104 
105  GenerateBoundaryConditionExpansion(graph, bcs, bcvar,
106  DeclareCoeffPhysArrays);
107  if (DeclareCoeffPhysArrays)
108  {
109  EvaluateBoundaryConditions(0.0, bcvar);
110  }
111  ApplyGeomInfo();
112  FindPeriodicTraces(bcs, bcvar);
113  }
114 
115  if (SetUpJustDG)
116  {
117  SetUpDG(variable);
118  }
119  else
120  {
121  int i, cnt;
122  Array<OneD, int> ElmtID, TraceID;
123  GetBoundaryToElmtMap(ElmtID, TraceID);
124 
125  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
126  {
128  locExpList = m_bndCondExpansions[i];
129 
130  for (int e = 0; e < locExpList->GetExpSize(); ++e)
131  {
132  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(TraceID[cnt + e],
133  locExpList->GetExp(e));
134  locExpList->GetExp(e)->SetAdjacentElementExp(
135  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
136  }
137  cnt += m_bndCondExpansions[i]->GetExpSize();
138  }
139 
140  if (m_session->DefinesSolverInfo("PROJECTION"))
141  {
142  std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
143  if ((ProjectStr == "MixedCGDG") ||
144  (ProjectStr == "Mixed_CG_Discontinuous"))
145  {
146  SetUpDG();
147  }
148  else
149  {
151  }
152  }
153  else
154  {
156  }
157  }
158 }
159 
160 /**
161  * @brief Set up all DG member variables and maps.
162  */
163 void DisContField::SetUpDG(const std::string variable)
164 {
165  // Check for multiple calls
167  {
168  return;
169  }
170 
171  // Set up matrix map
173 
174  // Set up trace space
177  m_comm);
178 
179  PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
180  : (m_expType == e2D) ? m_periodicEdges
181  : m_periodicFaces;
182 
185  m_bndConditions, periodicTraces, variable);
186 
187  if (m_session->DefinesCmdLineArgument("verbose"))
188  {
189  m_traceMap->PrintStats(std::cout, variable);
190  }
191 
193  m_traceMap->GetElmtToTrace();
194 
195  // Scatter trace to elements. For each element, we find
196  // the trace associated to each elemental trace. The element
197  // then retains a pointer to the elemental trace, to
198  // ensure uniqueness of normals when retrieving from two
199  // adjoining elements which do not lie in a plane.
200  for (int i = 0; i < m_exp->size(); ++i)
201  {
202  for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
203  {
204  LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];
205  ;
206 
207  exp->SetAdjacentElementExp(j, (*m_exp)[i]);
208  (*m_exp)[i]->SetTraceExp(j, exp);
209  }
210  }
211 
212  // Set up physical normals
214 
215  int cnt, n;
216 
217  // Identify boundary trace
218  for (cnt = 0, n = 0; n < m_bndCondExpansions.size(); ++n)
219  {
220  if (m_bndConditions[n]->GetBoundaryConditionType() !=
222  {
223  for (int v = 0; v < m_bndCondExpansions[n]->GetExpSize(); ++v)
224  {
225  m_boundaryTraces.insert(
226  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + v));
227  }
228  cnt += m_bndCondExpansions[n]->GetExpSize();
229  }
230  }
231 
232  // Set up information for periodic boundary conditions.
233  std::unordered_map<int, pair<int, int>> perTraceToExpMap;
234  for (cnt = n = 0; n < m_exp->size(); ++n)
235  {
236  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
237  {
238  auto it = periodicTraces.find((*m_exp)[n]->GetGeom()->GetTid(v));
239 
240  if (it != periodicTraces.end())
241  {
242  perTraceToExpMap[it->first] = make_pair(n, v);
243  }
244  }
245  }
246 
247  // Set up left-adjacent tracelist.
248  m_leftAdjacentTraces.resize(cnt);
249 
250  // count size of trace
251  for (cnt = n = 0; n < m_exp->size(); ++n)
252  {
253  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
254  {
256  }
257  }
258 
259  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
260  for (cnt = n = 0; n < m_exp->size(); ++n)
261  {
262  for (int v = 0; v < (*m_exp)[n]->GetNtraces(); ++v, ++cnt)
263  {
264  int GeomId = (*m_exp)[n]->GetGeom()->GetTid(v);
265 
266  // Check to see if this trace is periodic.
267  auto it = periodicTraces.find(GeomId);
268 
269  if (it != periodicTraces.end())
270  {
271  const PeriodicEntity &ent = it->second[0];
272  auto it2 = perTraceToExpMap.find(ent.id);
273 
274  if (it2 == perTraceToExpMap.end())
275  {
276  if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
277  !ent.isLocal)
278  {
279  continue;
280  }
281  else
282  {
283  ASSERTL1(false, "Periodic trace not found!");
284  }
285  }
286 
288  "Periodic trace in non-forward space?");
289 
290  int offset = m_trace->GetPhys_Offset(
291  (m_traceMap->GetElmtToTrace())[n][v]->GetElmtId());
292 
293  int offset2 = m_trace->GetPhys_Offset(
294  (m_traceMap->GetElmtToTrace())[it2->second.first]
295  [it2->second.second]
296  ->GetElmtId());
297 
298  switch (m_expType)
299  {
300  case e1D:
301  {
302  m_periodicFwdCopy.push_back(offset);
303  m_periodicBwdCopy.push_back(offset2);
304  }
305  break;
306  case e2D:
307  {
308  // Calculate relative orientations between trace to
309  // calculate copying map.
310  int nquad = elmtToTrace[n][v]->GetNumPoints(0);
311 
312  vector<int> tmpBwd(nquad);
313  vector<int> tmpFwd(nquad);
314 
315  if (ent.orient == StdRegions::eForwards)
316  {
317  for (int i = 0; i < nquad; ++i)
318  {
319  tmpBwd[i] = offset2 + i;
320  tmpFwd[i] = offset + i;
321  }
322  }
323  else
324  {
325  for (int i = 0; i < nquad; ++i)
326  {
327  tmpBwd[i] = offset2 + i;
328  tmpFwd[i] = offset + nquad - i - 1;
329  }
330  }
331 
332  for (int i = 0; i < nquad; ++i)
333  {
334  m_periodicFwdCopy.push_back(tmpFwd[i]);
335  m_periodicBwdCopy.push_back(tmpBwd[i]);
336  }
337  }
338  break;
339  case e3D:
340  {
341  // Calculate relative orientations between faces to
342  // calculate copying map.
343  int nquad1 = elmtToTrace[n][v]->GetNumPoints(0);
344  int nquad2 = elmtToTrace[n][v]->GetNumPoints(1);
345 
346  vector<int> tmpBwd(nquad1 * nquad2);
347  vector<int> tmpFwd(nquad1 * nquad2);
348 
349  if (ent.orient ==
351  ent.orient ==
353  ent.orient ==
356  {
357  for (int i = 0; i < nquad2; ++i)
358  {
359  for (int j = 0; j < nquad1; ++j)
360  {
361  tmpBwd[i * nquad1 + j] =
362  offset2 + i * nquad1 + j;
363  tmpFwd[i * nquad1 + j] =
364  offset + j * nquad2 + i;
365  }
366  }
367  }
368  else
369  {
370  for (int i = 0; i < nquad2; ++i)
371  {
372  for (int j = 0; j < nquad1; ++j)
373  {
374  tmpBwd[i * nquad1 + j] =
375  offset2 + i * nquad1 + j;
376  tmpFwd[i * nquad1 + j] =
377  offset + i * nquad1 + j;
378  }
379  }
380  }
381 
382  if (ent.orient ==
384  ent.orient ==
386  ent.orient ==
389  {
390  // Reverse x direction
391  for (int i = 0; i < nquad2; ++i)
392  {
393  for (int j = 0; j < nquad1 / 2; ++j)
394  {
395  swap(tmpFwd[i * nquad1 + j],
396  tmpFwd[i * nquad1 + nquad1 - j - 1]);
397  }
398  }
399  }
400 
401  if (ent.orient ==
403  ent.orient ==
405  ent.orient ==
408  {
409  // Reverse y direction
410  for (int j = 0; j < nquad1; ++j)
411  {
412  for (int i = 0; i < nquad2 / 2; ++i)
413  {
414  swap(tmpFwd[i * nquad1 + j],
415  tmpFwd[(nquad2 - i - 1) * nquad1 + j]);
416  }
417  }
418  }
419 
420  for (int i = 0; i < nquad1 * nquad2; ++i)
421  {
422  m_periodicFwdCopy.push_back(tmpFwd[i]);
423  m_periodicBwdCopy.push_back(tmpBwd[i]);
424  }
425  }
426  break;
427  default:
428  ASSERTL1(false, "not set up");
429  }
430  }
431  }
432  }
433 
435  *this, m_trace, elmtToTrace, m_leftAdjacentTraces);
436 }
437 
438 /**
439  * @brief This routine determines if an element is to the "left"
440  * of the adjacent trace, which arises from the idea there is a local
441  * normal direction between two elements (i.e. on the trace)
442  * and one elements would then be the left.
443  *
444  * This is typically required since we only define one normal
445  * direction along the trace which goes from the left adjacent
446  * element to the right adjacent element. It is also useful
447  * in DG discretisations to set up a local Riemann problem
448  * where we need to have the concept of a local normal
449  * direction which is unique between two elements.
450  *
451  * There are two cases to be checked:
452  *
453  * 1) First is the trace on a boundary condition (perioidic or
454  * otherwise) or on a partition boundary. If a partition
455  * boundary then trace is always considered to be the left
456  * adjacent trace and the normal is pointing outward of the
457  * soltuion domain. We have to perform an additional case for
458  * a periodic boundary where wer chose the element with the
459  * lowest global id. If the trace is on a parallel partition
460  * we use a member of the traceMap where one side is chosen to
461  * contribute to a unique map and have a value which is not
462  * -1 so this is the identifier for the left adjacent side
463  *
464  * 2) If the element is a elemental boundary on one element
465  * the left adjacent element is defined by a link to the left
466  * element from the trace expansion and this is consistent
467  * with the defitiion of the normal which is determined by the
468  * largest id (in contrast to the periodic boundary definition
469  * !). This reversal of convention does not really matter
470  * providing the normals are defined consistently.
471  */
472 
473 bool DisContField::IsLeftAdjacentTrace(const int n, const int e)
474 {
476  m_traceMap->GetElmtToTrace()[n][e];
477 
478  PeriodicMap periodicTraces = (m_expType == e1D) ? m_periodicVerts
479  : (m_expType == e2D) ? m_periodicEdges
480  : m_periodicFaces;
481 
482  bool fwd = true;
483  if (traceEl->GetLeftAdjacentElementTrace() == -1 ||
484  traceEl->GetRightAdjacentElementTrace() == -1)
485  {
486  // Boundary edge (1 connected element). Do nothing in
487  // serial.
488  auto it = m_boundaryTraces.find(traceEl->GetElmtId());
489 
490  // If the edge does not have a boundary condition set on
491  // it, then assume it is a partition edge or periodic.
492  if (it == m_boundaryTraces.end())
493  {
494  fwd = true;
495  }
496  }
497  else if (traceEl->GetLeftAdjacentElementTrace() != -1 &&
498  traceEl->GetRightAdjacentElementTrace() != -1)
499  {
500  // Non-boundary edge (2 connected elements).
501  fwd = (traceEl->GetLeftAdjacentElementExp().get()) == (*m_exp)[n].get();
502  }
503  else
504  {
505  ASSERTL2(false, "Unconnected trace element!");
506  }
507 
508  return fwd;
509 }
510 
511 // Given all boundary regions for the whole solution determine
512 // which ones (if any) are part of domain and proivde all other
513 // conditions are given as UserDefined Dirichlet.
515  const SpatialDomains::CompositeMap &domain,
517  const std::string &variable)
518 {
519  boost::ignore_unused(variable);
521 
522  returnval =
524 
525  map<int, int> GeometryToRegionsMap;
526 
528  Allbcs.GetBoundaryRegions();
530  Allbcs.GetBoundaryConditions();
531 
532  // Set up a map of all boundary regions
533  for (auto &it : bregions)
534  {
535  for (auto &bregionIt : *it.second)
536  {
537  // can assume that all regions only contain one point in 1D
538  // Really do not need loop above
539  int id = bregionIt.second->m_geomVec[0]->GetGlobalID();
540  GeometryToRegionsMap[id] = it.first;
541  }
542  }
543 
544  map<int, SpatialDomains::GeometrySharedPtr> EndOfDomain;
545 
546  // Now find out which points in domain have only one vertex
547  for (auto &domIt : domain)
548  {
549  SpatialDomains::CompositeSharedPtr geomvector = domIt.second;
550  for (int i = 0; i < geomvector->m_geomVec.size(); ++i)
551  {
552  for (int j = 0; j < 2; ++j)
553  {
554  int vid = geomvector->m_geomVec[i]->GetVid(j);
555  if (EndOfDomain.count(vid) == 0)
556  {
557  EndOfDomain[vid] = geomvector->m_geomVec[i]->GetVertex(j);
558  }
559  else
560  {
561  EndOfDomain.erase(vid);
562  }
563  }
564  }
565  }
566  ASSERTL1(EndOfDomain.size() == 2, "Did not find two ends of domain");
567 
568  for (auto &regIt : EndOfDomain)
569  {
570  if (GeometryToRegionsMap.count(regIt.first) != 0)
571  {
572  // Set up boundary condition up
573  auto iter = GeometryToRegionsMap.find(regIt.first);
574  ASSERTL1(iter != GeometryToRegionsMap.end(),
575  "Failied to find GeometryToRegionMap");
576 
577  int regionId = iter->second;
578  auto bregionsIter = bregions.find(regionId);
579  ASSERTL1(bregionsIter != bregions.end(),
580  "Failed to find boundary region");
581 
582  SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
583  returnval->AddBoundaryRegions(regionId, breg);
584 
585  auto bconditionsIter = bconditions.find(regionId);
586  ASSERTL1(bconditionsIter != bconditions.end(),
587  "Failed to find boundary collection");
588 
590  bconditionsIter->second;
591  returnval->AddBoundaryConditions(regionId, bcond);
592  }
593  }
594 
595  return returnval;
596 }
597 
598 /**
599  * Constructor for use in multidomain computations where a
600  * domain list can be passed instead of graph1D
601  *
602  * @param domain Subdomain specified in the inputfile from
603  * which the DisContField is set up
604  */
606  const SpatialDomains::MeshGraphSharedPtr &graph1D,
607  const SpatialDomains::CompositeMap &domain,
609  const std::string &variable,
610  const LibUtilities::CommSharedPtr &comm,
611  bool SetToOneSpaceDimension,
612  const Collections::ImplementationType ImpType)
613  : ExpList(pSession, domain, graph1D, true, variable, SetToOneSpaceDimension,
614  comm, ImpType)
615 {
616  if (variable.compare("DefaultVar") != 0)
617  {
619  GetDomainBCs(domain, Allbcs, variable);
620 
621  GenerateBoundaryConditionExpansion(m_graph, *DomBCs, variable);
622  EvaluateBoundaryConditions(0.0, variable);
623  ApplyGeomInfo();
624  FindPeriodicTraces(*DomBCs, variable);
625  }
626 
627  SetUpDG(variable);
628 }
629 
630 /**
631  * Constructs a field as a copy of an existing field.
632  * @param In Existing DisContField object to copy.
633  */
635  const bool DeclareCoeffPhysArrays)
636  : ExpList(In, DeclareCoeffPhysArrays),
637  m_bndCondExpansions(In.m_bndCondExpansions),
638  m_bndConditions(In.m_bndConditions), m_globalBndMat(In.m_globalBndMat),
639  m_traceMap(In.m_traceMap), m_boundaryTraces(In.m_boundaryTraces),
640  m_periodicVerts(In.m_periodicVerts),
641  m_periodicFwdCopy(In.m_periodicFwdCopy),
642  m_periodicBwdCopy(In.m_periodicBwdCopy),
643  m_leftAdjacentTraces(In.m_leftAdjacentTraces),
644  m_locTraceToTraceMap(In.m_locTraceToTraceMap)
645 {
646  if (In.m_trace)
647  {
649  *In.m_trace, DeclareCoeffPhysArrays);
650  }
651 }
652 
653 /*
654  * @brief Copy type constructor which declares new boundary conditions
655  * and re-uses mapping info and trace space if possible
656  */
659  const std::string &variable, const bool SetUpJustDG,
660  const bool DeclareCoeffPhysArrays)
661  : ExpList(In, DeclareCoeffPhysArrays)
662 {
663 
665 
666  // Set up boundary conditions for this variable.
667  // Do not set up BCs if default variable
668  if (variable.compare("DefaultVar") != 0)
669  {
671  GenerateBoundaryConditionExpansion(graph, bcs, variable);
672 
673  if (DeclareCoeffPhysArrays)
674  {
675  EvaluateBoundaryConditions(0.0, variable);
676  }
677 
679  {
680  // Find periodic edges for this variable.
681  FindPeriodicTraces(bcs, variable);
682 
683  if (SetUpJustDG)
684  {
685  SetUpDG(variable);
686  }
687  else
688  {
689  // set elmt edges to point to robin bc edges if required
690  int i, cnt = 0;
691  Array<OneD, int> ElmtID, TraceID;
692  GetBoundaryToElmtMap(ElmtID, TraceID);
693 
694  for (i = 0; i < m_bndCondExpansions.size(); ++i)
695  {
697 
698  int e;
699  locExpList = m_bndCondExpansions[i];
700 
701  for (e = 0; e < locExpList->GetExpSize(); ++e)
702  {
703  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
704  TraceID[cnt + e], locExpList->GetExp(e));
705  locExpList->GetExp(e)->SetAdjacentElementExp(
706  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
707  }
708 
709  cnt += m_bndCondExpansions[i]->GetExpSize();
710  }
711 
712  if (m_session->DefinesSolverInfo("PROJECTION"))
713  {
714  std::string ProjectStr =
715  m_session->GetSolverInfo("PROJECTION");
716 
717  if ((ProjectStr == "MixedCGDG") ||
718  (ProjectStr == "Mixed_CG_Discontinuous"))
719  {
720  SetUpDG();
721  }
722  else
723  {
725  }
726  }
727  else
728  {
730  }
731  }
732  }
733  else
734  {
736  m_trace = In.m_trace;
737  m_traceMap = In.m_traceMap;
746 
747  if (SetUpJustDG == false)
748  {
749  // set elmt edges to point to robin bc edges if required
750  int i, cnt = 0;
751  Array<OneD, int> ElmtID, TraceID;
752  GetBoundaryToElmtMap(ElmtID, TraceID);
753 
754  for (i = 0; i < m_bndCondExpansions.size(); ++i)
755  {
757 
758  int e;
759  locExpList = m_bndCondExpansions[i];
760 
761  for (e = 0; e < locExpList->GetExpSize(); ++e)
762  {
763  (*m_exp)[ElmtID[cnt + e]]->SetTraceExp(
764  TraceID[cnt + e], locExpList->GetExp(e));
765  locExpList->GetExp(e)->SetAdjacentElementExp(
766  TraceID[cnt + e], (*m_exp)[ElmtID[cnt + e]]);
767  }
768  cnt += m_bndCondExpansions[i]->GetExpSize();
769  }
770 
772  }
773  }
774  }
775 }
776 
777 /**
778  * Constructs a field as a copy of an existing explist field.
779  * @param In Existing ExpList object to copy.
780  */
782 {
783 }
784 
785 /**
786  *
787  */
789 {
790 }
791 
792 /**
793  * \brief This function discretises the boundary conditions by setting
794  * up a list of one-dimensions lower boundary expansions.
795  *
796  * According to their boundary region, the separate boundary
797  * expansions are bundled together in an object of the class
798  *
799  * @param graph A mesh, containing information about the domain
800  * and the spectral/hp element expansions.
801  * @param bcs Information about the enforced boundary
802  * conditions.
803  * @param variable The session variable associated with the
804  * boundary conditions to enforce.
805  * @param DeclareCoeffPhysArrays bool to identify if array
806  * space should be setup.
807  * Default is true.
808  */
811  const SpatialDomains::BoundaryConditions &bcs, const std::string variable,
812  const bool DeclareCoeffPhysArrays)
813 {
814  int cnt = 0;
818  bcs.GetBoundaryRegions();
820  bcs.GetBoundaryConditions();
821 
826 
827  m_bndCondBndWeight = Array<OneD, NekDouble>{bregions.size(), 0.0};
828 
829  // count the number of non-periodic boundary points
830  for (auto &it : bregions)
831  {
832  bc = GetBoundaryCondition(bconditions, it.first, variable);
833 
835  m_session, *(it.second), graph, DeclareCoeffPhysArrays, variable,
836  false, bc->GetComm());
837 
838  m_bndCondExpansions[cnt] = locExpList;
839  m_bndConditions[cnt] = bc;
840 
841  std::string type = m_bndConditions[cnt]->GetUserDefined();
842 
843  // Set up normals on non-Dirichlet boundary conditions. Second
844  // two conditions ideally should be in local solver setup (when
845  // made into factory)
846  if (bc->GetBoundaryConditionType() != SpatialDomains::eDirichlet ||
847  boost::iequals(type, "I") || boost::iequals(type, "CalcBC"))
848  {
850  }
851  cnt++;
852  }
853 }
854 
855 /**
856  * @brief Determine the periodic faces, edges and vertices for the given
857  * graph.
858  *
859  * @param bcs Information about the boundary conditions.
860  * @param variable Specifies the field.
861  */
863  const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
864 {
866  bcs.GetBoundaryRegions();
868  bcs.GetBoundaryConditions();
869 
870  LibUtilities::CommSharedPtr vComm = m_comm->GetRowComm();
871 
872  switch (m_expType)
873  {
874  case e1D:
875  {
876  int i, region1ID, region2ID;
877 
879  map<int, int> BregionToVertMap;
880 
881  // Construct list of all periodic Region and their
882  // global vertex on this process.
883  for (auto &it : bregions)
884  {
885  locBCond =
886  GetBoundaryCondition(bconditions, it.first, variable);
887 
888  if (locBCond->GetBoundaryConditionType() !=
890  {
891  continue;
892  }
893  int id =
894  it.second->begin()->second->m_geomVec[0]->GetGlobalID();
895 
896  BregionToVertMap[it.first] = id;
897  }
898 
899  set<int> islocal;
900 
901  int n = vComm->GetSize();
902  int p = vComm->GetRank();
903 
904  Array<OneD, int> nregions(n, 0);
905  nregions[p] = BregionToVertMap.size();
906  vComm->AllReduce(nregions, LibUtilities::ReduceSum);
907 
908  int totRegions = Vmath::Vsum(n, nregions, 1);
909 
910  Array<OneD, int> regOffset(n, 0);
911 
912  for (i = 1; i < n; ++i)
913  {
914  regOffset[i] = regOffset[i - 1] + nregions[i - 1];
915  }
916 
917  Array<OneD, int> bregmap(totRegions, 0);
918  Array<OneD, int> bregid(totRegions, 0);
919 
920  i = regOffset[p];
921  for (auto &iit : BregionToVertMap)
922  {
923  bregid[i] = iit.first;
924  bregmap[i++] = iit.second;
925  islocal.insert(iit.first);
926  }
927 
928  vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
929  vComm->AllReduce(bregid, LibUtilities::ReduceSum);
930 
931  for (int i = 0; i < totRegions; ++i)
932  {
933  BregionToVertMap[bregid[i]] = bregmap[i];
934  }
935 
936  // Construct list of all periodic pairs local to this process.
937  for (auto &it : bregions)
938  {
939  locBCond =
940  GetBoundaryCondition(bconditions, it.first, variable);
941 
942  if (locBCond->GetBoundaryConditionType() !=
944  {
945  continue;
946  }
947 
948  // Identify periodic boundary region IDs.
949  region1ID = it.first;
950  region2ID =
951  std::static_pointer_cast<
953  ->m_connectedBoundaryRegion;
954 
955  ASSERTL0(BregionToVertMap.count(region1ID) != 0,
956  "Cannot determine vertex of region1ID");
957 
958  ASSERTL0(BregionToVertMap.count(region2ID) != 0,
959  "Cannot determine vertex of region2ID");
960 
961  PeriodicEntity ent(BregionToVertMap[region2ID],
963  islocal.count(region2ID) != 0);
964 
965  m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
966  }
967  }
968  break;
969  case e2D:
970  {
971  int region1ID, region2ID, i, j, k, cnt;
973 
975  m_graph->GetCompositeOrdering();
977  m_graph->GetBndRegionOrdering();
978  SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
979 
980  // Unique collection of pairs of periodic composites
981  // (i.e. if composites 1 and 2 are periodic then this
982  // map will contain either the pair (1,2) or (2,1) but
983  // not both).
984  map<int, int> perComps;
985  map<int, vector<int>> allVerts;
986  set<int> locVerts;
987  map<int, pair<int, StdRegions::Orientation>> allEdges;
988 
989  // Set up a set of all local verts and edges.
990  for (i = 0; i < (*m_exp).size(); ++i)
991  {
992  for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
993  {
994  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
995  locVerts.insert(id);
996  }
997  }
998 
999  // Construct list of all periodic pairs local to this process.
1000  for (auto &it : bregions)
1001  {
1002  locBCond =
1003  GetBoundaryCondition(bconditions, it.first, variable);
1004 
1005  if (locBCond->GetBoundaryConditionType() !=
1007  {
1008  continue;
1009  }
1010 
1011  // Identify periodic boundary region IDs.
1012  region1ID = it.first;
1013  region2ID =
1014  std::static_pointer_cast<
1016  ->m_connectedBoundaryRegion;
1017 
1018  // From this identify composites. Note that in
1019  // serial this will be an empty map.
1020  int cId1, cId2;
1021  if (vComm->GetSize() == 1)
1022  {
1023  cId1 = it.second->begin()->first;
1024  cId2 = bregions.find(region2ID)->second->begin()->first;
1025  }
1026  else
1027  {
1028  cId1 = bndRegOrder.find(region1ID)->second[0];
1029  cId2 = bndRegOrder.find(region2ID)->second[0];
1030  }
1031 
1032  ASSERTL0(it.second->size() == 1,
1033  "Boundary region " +
1034  boost::lexical_cast<string>(region1ID) +
1035  " should only contain 1 composite.");
1036 
1037  vector<unsigned int> tmpOrder;
1038 
1039  // Construct set containing all periodic edgesd on
1040  // this process
1042  it.second->begin()->second;
1043 
1044  for (i = 0; i < c->m_geomVec.size(); ++i)
1045  {
1047  std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1048  c->m_geomVec[i]);
1049  ASSERTL0(segGeom, "Unable to cast to shared ptr");
1050 
1052  m_graph->GetElementsFromEdge(segGeom);
1053  ASSERTL0(elmt->size() == 1,
1054  "The periodic boundaries belong to "
1055  "more than one element of the mesh");
1056 
1058  std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1059  elmt->at(0).first);
1060 
1061  allEdges[c->m_geomVec[i]->GetGlobalID()] =
1062  make_pair(elmt->at(0).second,
1063  geom->GetEorient(elmt->at(0).second));
1064 
1065  // In serial mesh partitioning will not have occurred so
1066  // need to fill composite ordering map manually.
1067  if (vComm->GetSize() == 1)
1068  {
1069  tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1070  }
1071 
1072  vector<int> vertList(2);
1073  vertList[0] = segGeom->GetVid(0);
1074  vertList[1] = segGeom->GetVid(1);
1075  allVerts[c->m_geomVec[i]->GetGlobalID()] = vertList;
1076  }
1077 
1078  if (vComm->GetSize() == 1)
1079  {
1080  compOrder[it.second->begin()->first] = tmpOrder;
1081  }
1082 
1083  // See if we already have either region1 or
1084  // region2 stored in perComps map.
1085  if (perComps.count(cId1) == 0)
1086  {
1087  if (perComps.count(cId2) == 0)
1088  {
1089  perComps[cId1] = cId2;
1090  }
1091  else
1092  {
1093  std::stringstream ss;
1094  ss << "Boundary region " << cId2 << " should be "
1095  << "periodic with " << perComps[cId2] << " but "
1096  << "found " << cId1 << " instead!";
1097  ASSERTL0(perComps[cId2] == cId1, ss.str());
1098  }
1099  }
1100  else
1101  {
1102  std::stringstream ss;
1103  ss << "Boundary region " << cId1 << " should be "
1104  << "periodic with " << perComps[cId1] << " but "
1105  << "found " << cId2 << " instead!";
1106  ASSERTL0(perComps[cId1] == cId1, ss.str());
1107  }
1108  }
1109 
1110  // In case of periodic partition being split over many composites
1111  // may not have all periodic matches so check this
1112  int idmax = -1;
1113  for (auto &cIt : perComps)
1114  {
1115  idmax = max(idmax, cIt.first);
1116  idmax = max(idmax, cIt.second);
1117  }
1118  vComm->AllReduce(idmax, LibUtilities::ReduceMax);
1119  idmax++;
1120  Array<OneD, int> perid(idmax, -1);
1121  for (auto &cIt : perComps)
1122  {
1123  perid[cIt.first] = cIt.second;
1124  }
1125  vComm->AllReduce(perid, LibUtilities::ReduceMax);
1126  // update all partitions
1127  for (int i = 0; i < idmax; ++i)
1128  {
1129  if (perid[i] > -1)
1130  {
1131  // skip if complemetary relationship has
1132  // already been speficied in list
1133  if (perComps.count(perid[i]))
1134  {
1135  continue;
1136  }
1137  perComps[i] = perid[i];
1138  }
1139  }
1140 
1141  // Process local edge list to obtain relative edge
1142  // orientations.
1143  int n = vComm->GetSize();
1144  int p = vComm->GetRank();
1145  int totEdges;
1146  Array<OneD, int> edgecounts(n, 0);
1147  Array<OneD, int> edgeoffset(n, 0);
1148  Array<OneD, int> vertoffset(n, 0);
1149 
1150  edgecounts[p] = allEdges.size();
1151  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
1152 
1153  edgeoffset[0] = 0;
1154  for (i = 1; i < n; ++i)
1155  {
1156  edgeoffset[i] = edgeoffset[i - 1] + edgecounts[i - 1];
1157  }
1158 
1159  totEdges = Vmath::Vsum(n, edgecounts, 1);
1160  Array<OneD, int> edgeIds(totEdges, 0);
1161  Array<OneD, int> edgeIdx(totEdges, 0);
1162  Array<OneD, int> edgeOrient(totEdges, 0);
1163  Array<OneD, int> edgeVerts(totEdges, 0);
1164 
1165  auto sIt = allEdges.begin();
1166 
1167  for (i = 0; sIt != allEdges.end(); ++sIt)
1168  {
1169  edgeIds[edgeoffset[p] + i] = sIt->first;
1170  edgeIdx[edgeoffset[p] + i] = sIt->second.first;
1171  edgeOrient[edgeoffset[p] + i] = sIt->second.second;
1172  edgeVerts[edgeoffset[p] + i++] = allVerts[sIt->first].size();
1173  }
1174 
1175  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1176  vComm->AllReduce(edgeIdx, LibUtilities::ReduceSum);
1177  vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
1178  vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
1179 
1180  // Calculate number of vertices on each processor.
1181  Array<OneD, int> procVerts(n, 0);
1182  int nTotVerts;
1183 
1184  // Note if there are no periodic edges at all calling Vsum will
1185  // cause a segfault.
1186  if (totEdges > 0)
1187  {
1188  nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
1189  }
1190  else
1191  {
1192  nTotVerts = 0;
1193  }
1194 
1195  for (i = 0; i < n; ++i)
1196  {
1197  if (edgecounts[i] > 0)
1198  {
1199  procVerts[i] = Vmath::Vsum(edgecounts[i],
1200  edgeVerts + edgeoffset[i], 1);
1201  }
1202  else
1203  {
1204  procVerts[i] = 0;
1205  }
1206  }
1207  vertoffset[0] = 0;
1208 
1209  for (i = 1; i < n; ++i)
1210  {
1211  vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1212  }
1213 
1214  Array<OneD, int> traceIds(nTotVerts, 0);
1215  for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
1216  {
1217  for (j = 0; j < allVerts[sIt->first].size(); ++j)
1218  {
1219  traceIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
1220  }
1221  }
1222 
1223  vComm->AllReduce(traceIds, LibUtilities::ReduceSum);
1224 
1225  // For simplicity's sake create a map of edge id -> orientation.
1226  map<int, StdRegions::Orientation> orientMap;
1227  map<int, vector<int>> vertMap;
1228 
1229  for (cnt = i = 0; i < totEdges; ++i)
1230  {
1231  ASSERTL0(orientMap.count(edgeIds[i]) == 0,
1232  "Already found edge in orientation map!");
1233 
1234  // Work out relative orientations. To avoid having
1235  // to exchange vertex locations, we figure out if
1236  // the edges are backwards or forwards orientated
1237  // with respect to a forwards orientation that is
1238  // CCW. Since our local geometries are
1239  // forwards-orientated with respect to the
1240  // Cartesian axes, we need to invert the
1241  // orientation for the top and left edges of a
1242  // quad and the left edge of a triangle.
1244  (StdRegions::Orientation)edgeOrient[i];
1245 
1246  if (edgeIdx[i] > 1)
1247  {
1250  }
1251 
1252  orientMap[edgeIds[i]] = o;
1253 
1254  vector<int> verts(edgeVerts[i]);
1255 
1256  for (j = 0; j < edgeVerts[i]; ++j)
1257  {
1258  verts[j] = traceIds[cnt++];
1259  }
1260  vertMap[edgeIds[i]] = verts;
1261  }
1262 
1263  // Go through list of composites and figure out which
1264  // edges are parallel from original ordering in
1265  // session file. This includes composites which are
1266  // not necessarily on this process.
1267  map<int, int> allCompPairs;
1268 
1269  // Store temporary map of periodic vertices which will hold all
1270  // periodic vertices on the entire mesh so that doubly periodic
1271  // vertices can be counted properly across partitions. Local
1272  // vertices are copied into m_periodicVerts at the end of the
1273  // function.
1274  PeriodicMap periodicVerts;
1275 
1276  for (auto &cIt : perComps)
1277  {
1279  const int id1 = cIt.first;
1280  const int id2 = cIt.second;
1281  std::string id1s = boost::lexical_cast<string>(id1);
1282  std::string id2s = boost::lexical_cast<string>(id2);
1283 
1284  if (compMap.count(id1) > 0)
1285  {
1286  c[0] = compMap[id1];
1287  }
1288 
1289  if (compMap.count(id2) > 0)
1290  {
1291  c[1] = compMap[id2];
1292  }
1293 
1294  // Loop over composite ordering to construct list of all
1295  // periodic edges regardless of whether they are on this
1296  // process.
1297  map<int, int> compPairs;
1298 
1299  ASSERTL0(compOrder.count(id1) > 0,
1300  "Unable to find composite " + id1s + " in order map.");
1301  ASSERTL0(compOrder.count(id2) > 0,
1302  "Unable to find composite " + id2s + " in order map.");
1303  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1304  "Periodic composites " + id1s + " and " + id2s +
1305  " should have the same number of elements.");
1306  ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
1307  id1s + " and " + id2s +
1308  " are empty!");
1309 
1310  // TODO: Add more checks.
1311  for (i = 0; i < compOrder[id1].size(); ++i)
1312  {
1313  int eId1 = compOrder[id1][i];
1314  int eId2 = compOrder[id2][i];
1315 
1316  ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
1317 
1318  if (compPairs.count(eId2) != 0)
1319  {
1320  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
1321  }
1322  compPairs[eId1] = eId2;
1323  }
1324 
1325  // Construct set of all edges that we have
1326  // locally on this processor.
1327  set<int> locEdges;
1328  for (i = 0; i < 2; ++i)
1329  {
1330  if (!c[i])
1331  {
1332  continue;
1333  }
1334 
1335  if (c[i]->m_geomVec.size() > 0)
1336  {
1337  for (j = 0; j < c[i]->m_geomVec.size(); ++j)
1338  {
1339  locEdges.insert(c[i]->m_geomVec[j]->GetGlobalID());
1340  }
1341  }
1342  }
1343 
1344  // Loop over all edges in the geometry composite.
1345  for (auto &pIt : compPairs)
1346  {
1347  int ids[2] = {pIt.first, pIt.second};
1348  bool local[2] = {locEdges.count(pIt.first) > 0,
1349  locEdges.count(pIt.second) > 0};
1350 
1351  ASSERTL0(orientMap.count(ids[0]) > 0 &&
1352  orientMap.count(ids[1]) > 0,
1353  "Unable to find edge in orientation map.");
1354 
1355  allCompPairs[pIt.first] = pIt.second;
1356  allCompPairs[pIt.second] = pIt.first;
1357 
1358  for (i = 0; i < 2; ++i)
1359  {
1360  if (!local[i])
1361  {
1362  continue;
1363  }
1364 
1365  int other = (i + 1) % 2;
1366 
1368  orientMap[ids[i]] == orientMap[ids[other]]
1371 
1372  PeriodicEntity ent(ids[other], o, local[other]);
1373  m_periodicEdges[ids[i]].push_back(ent);
1374  }
1375 
1376  for (i = 0; i < 2; ++i)
1377  {
1378  int other = (i + 1) % 2;
1379 
1381  orientMap[ids[i]] == orientMap[ids[other]]
1384 
1385  // Determine periodic vertices.
1386  vector<int> perVerts1 = vertMap[ids[i]];
1387  vector<int> perVerts2 = vertMap[ids[other]];
1388 
1389  map<int, pair<int, bool>> tmpMap;
1390  if (o == StdRegions::eForwards)
1391  {
1392  tmpMap[perVerts1[0]] = make_pair(
1393  perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1394  tmpMap[perVerts1[1]] = make_pair(
1395  perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1396  }
1397  else
1398  {
1399  tmpMap[perVerts1[0]] = make_pair(
1400  perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1401  tmpMap[perVerts1[1]] = make_pair(
1402  perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1403  }
1404 
1405  for (auto &mIt : tmpMap)
1406  {
1407  // See if this vertex has been recorded already.
1408  PeriodicEntity ent2(mIt.second.first,
1410  mIt.second.second);
1411  auto perIt = periodicVerts.find(mIt.first);
1412 
1413  if (perIt == periodicVerts.end())
1414  {
1415  periodicVerts[mIt.first].push_back(ent2);
1416  perIt = periodicVerts.find(mIt.first);
1417  }
1418  else
1419  {
1420  bool doAdd = true;
1421  for (j = 0; j < perIt->second.size(); ++j)
1422  {
1423  if (perIt->second[j].id == mIt.second.first)
1424  {
1425  doAdd = false;
1426  break;
1427  }
1428  }
1429 
1430  if (doAdd)
1431  {
1432  perIt->second.push_back(ent2);
1433  }
1434  }
1435  }
1436  }
1437  }
1438  }
1439 
1440  // Search for periodic vertices and edges which are not in
1441  // a periodic composite but lie in this process. First, loop
1442  // over all information we have from other processors.
1443  for (cnt = i = 0; i < totEdges; ++i)
1444  {
1445  int edgeId = edgeIds[i];
1446 
1447  ASSERTL0(allCompPairs.count(edgeId) > 0,
1448  "Unable to find matching periodic edge.");
1449 
1450  int perEdgeId = allCompPairs[edgeId];
1451 
1452  for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1453  {
1454  int vId = traceIds[cnt];
1455 
1456  auto perId = periodicVerts.find(vId);
1457 
1458  if (perId == periodicVerts.end())
1459  {
1460  // This vertex is not included in the
1461  // map. Figure out which vertex it is
1462  // supposed to be periodic with. perEdgeId
1463  // is the edge ID which is periodic with
1464  // edgeId. The logic is much the same as
1465  // the loop above.
1466  int perVertexId =
1467  orientMap[edgeId] == orientMap[perEdgeId]
1468  ? vertMap[perEdgeId][(j + 1) % 2]
1469  : vertMap[perEdgeId][j];
1470 
1471  PeriodicEntity ent(perVertexId,
1473  locVerts.count(perVertexId) > 0);
1474 
1475  periodicVerts[vId].push_back(ent);
1476  }
1477  }
1478  }
1479 
1480  // Loop over all periodic vertices to complete connectivity
1481  // information.
1482  for (auto &perIt : periodicVerts)
1483  {
1484  // Loop over associated vertices.
1485  for (i = 0; i < perIt.second.size(); ++i)
1486  {
1487  auto perIt2 = periodicVerts.find(perIt.second[i].id);
1488  ASSERTL0(perIt2 != periodicVerts.end(),
1489  "Couldn't find periodic vertex.");
1490 
1491  for (j = 0; j < perIt2->second.size(); ++j)
1492  {
1493  if (perIt2->second[j].id == perIt.first)
1494  {
1495  continue;
1496  }
1497 
1498  bool doAdd = true;
1499  for (k = 0; k < perIt.second.size(); ++k)
1500  {
1501  if (perIt2->second[j].id == perIt.second[k].id)
1502  {
1503  doAdd = false;
1504  break;
1505  }
1506  }
1507 
1508  if (doAdd)
1509  {
1510  perIt.second.push_back(perIt2->second[j]);
1511  }
1512  }
1513  }
1514  }
1515 
1516  // Do one final loop over periodic vertices to remove non-local
1517  // vertices from map.
1518  for (auto &perIt : periodicVerts)
1519  {
1520  if (locVerts.count(perIt.first) > 0)
1521  {
1522  m_periodicVerts.insert(perIt);
1523  }
1524  }
1525  }
1526  break;
1527  case e3D:
1528  {
1530  m_graph->GetCompositeOrdering();
1531  SpatialDomains::BndRegionOrdering bndRegOrder =
1532  m_graph->GetBndRegionOrdering();
1533  SpatialDomains::CompositeMap compMap = m_graph->GetComposites();
1534 
1535  // perComps: Stores a unique collection of pairs of periodic
1536  // composites (i.e. if composites 1 and 2 are periodic then this map
1537  // will contain either the pair (1,2) or (2,1) but not both).
1538  //
1539  // The four maps allVerts, allCoord, allEdges and allOrient map a
1540  // periodic face to a vector containing the vertex ids of the face;
1541  // their coordinates; the edge ids of the face; and their
1542  // orientation within that face respectively.
1543  //
1544  // Finally the three sets locVerts, locEdges and locFaces store any
1545  // vertices, edges and faces that belong to a periodic composite and
1546  // lie on this process.
1547  map<int, RotPeriodicInfo> rotComp;
1548  map<int, int> perComps;
1549  map<int, vector<int>> allVerts;
1550  map<int, SpatialDomains::PointGeomVector> allCoord;
1551  map<int, vector<int>> allEdges;
1552  map<int, vector<StdRegions::Orientation>> allOrient;
1553  set<int> locVerts;
1554  set<int> locEdges;
1555  set<int> locFaces;
1556 
1557  int region1ID, region2ID, i, j, k, cnt;
1559 
1560  // Set up a set of all local verts and edges.
1561  for (i = 0; i < (*m_exp).size(); ++i)
1562  {
1563  for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
1564  {
1565  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
1566  locVerts.insert(id);
1567  }
1568 
1569  for (j = 0; j < (*m_exp)[i]->GetGeom()->GetNumEdges(); ++j)
1570  {
1571  int id = (*m_exp)[i]->GetGeom()->GetEid(j);
1572  locEdges.insert(id);
1573  }
1574  }
1575 
1576  // Begin by populating the perComps map. We loop over all periodic
1577  // boundary conditions and determine the composite associated with
1578  // it, then fill out the all* maps.
1579  for (auto &it : bregions)
1580  {
1581 
1582  locBCond =
1583  GetBoundaryCondition(bconditions, it.first, variable);
1584 
1585  if (locBCond->GetBoundaryConditionType() !=
1587  {
1588  continue;
1589  }
1590 
1591  // Identify periodic boundary region IDs.
1592  region1ID = it.first;
1593  region2ID =
1594  std::static_pointer_cast<
1596  ->m_connectedBoundaryRegion;
1597 
1598  // Check the region only contains a single composite.
1599  ASSERTL0(it.second->size() == 1,
1600  "Boundary region " +
1601  boost::lexical_cast<string>(region1ID) +
1602  " should only contain 1 composite.");
1603 
1604  // From this identify composites by looking at the original
1605  // boundary region ordering. Note that in serial the mesh
1606  // partitioner is not run, so this map will be empty and
1607  // therefore needs to be populated by using the corresponding
1608  // boundary region.
1609  int cId1, cId2;
1610  if (vComm->GetSize() == 1)
1611  {
1612  cId1 = it.second->begin()->first;
1613  cId2 = bregions.find(region2ID)->second->begin()->first;
1614  }
1615  else
1616  {
1617  cId1 = bndRegOrder.find(region1ID)->second[0];
1618  cId2 = bndRegOrder.find(region2ID)->second[0];
1619  }
1620 
1621  // check to see if boundary is rotationally aligned
1622  if (boost::icontains(locBCond->GetUserDefined(), "Rotated"))
1623  {
1624  vector<string> tmpstr;
1625 
1626  boost::split(tmpstr, locBCond->GetUserDefined(),
1627  boost::is_any_of(":"));
1628 
1629  if (boost::iequals(tmpstr[0], "Rotated"))
1630  {
1631  ASSERTL1(tmpstr.size() > 2,
1632  "Expected Rotated user defined string to "
1633  "contain direction and rotation angle "
1634  "and optionally a tolerance, "
1635  "i.e. Rotated:dir:PI/2:1e-6");
1636 
1637  ASSERTL1((tmpstr[1] == "x") || (tmpstr[1] == "y") ||
1638  (tmpstr[1] == "z"),
1639  "Rotated Dir is "
1640  "not specified as x,y or z");
1641 
1642  RotPeriodicInfo RotInfo;
1643  RotInfo.m_dir = (tmpstr[1] == "x") ? 0
1644  : (tmpstr[1] == "y") ? 1
1645  : 2;
1646 
1647  LibUtilities::Interpreter strEval;
1648  int ExprId = strEval.DefineFunction("", tmpstr[2]);
1649  RotInfo.m_angle = strEval.Evaluate(ExprId);
1650 
1651  if (tmpstr.size() == 4)
1652  {
1653  try
1654  {
1655  RotInfo.m_tol =
1656  boost::lexical_cast<NekDouble>(tmpstr[3]);
1657  }
1658  catch (...)
1659  {
1661  "failed to cast tolerance input "
1662  "to a double value in Rotated"
1663  "boundary information");
1664  }
1665  }
1666  else
1667  {
1668  RotInfo.m_tol = 1e-8;
1669  }
1670  rotComp[cId1] = RotInfo;
1671  }
1672  }
1673 
1675  it.second->begin()->second;
1676 
1677  vector<unsigned int> tmpOrder;
1678 
1679  // store the rotation info of this
1680 
1681  // From the composite, we now construct the allVerts, allEdges
1682  // and allCoord map so that they can be transferred across
1683  // processors. We also populate the locFaces set to store a
1684  // record of all faces local to this process.
1685  for (i = 0; i < c->m_geomVec.size(); ++i)
1686  {
1688  std::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
1689  c->m_geomVec[i]);
1690  ASSERTL1(faceGeom, "Unable to cast to shared ptr");
1691 
1692  // Get geometry ID of this face and store in locFaces.
1693  int faceId = c->m_geomVec[i]->GetGlobalID();
1694  locFaces.insert(faceId);
1695 
1696  // In serial, mesh partitioning will not have occurred so
1697  // need to fill composite ordering map manually.
1698  if (vComm->GetSize() == 1)
1699  {
1700  tmpOrder.push_back(c->m_geomVec[i]->GetGlobalID());
1701  }
1702 
1703  // Loop over vertices and edges of the face to populate
1704  // allVerts, allEdges and allCoord maps.
1705  vector<int> vertList, edgeList;
1707  vector<StdRegions::Orientation> orientVec;
1708  for (j = 0; j < faceGeom->GetNumVerts(); ++j)
1709  {
1710  vertList.push_back(faceGeom->GetVid(j));
1711  edgeList.push_back(faceGeom->GetEid(j));
1712  coordVec.push_back(faceGeom->GetVertex(j));
1713  orientVec.push_back(faceGeom->GetEorient(j));
1714  }
1715 
1716  allVerts[faceId] = vertList;
1717  allEdges[faceId] = edgeList;
1718  allCoord[faceId] = coordVec;
1719  allOrient[faceId] = orientVec;
1720  }
1721 
1722  // In serial, record the composite ordering in compOrder for
1723  // later in the routine.
1724  if (vComm->GetSize() == 1)
1725  {
1726  compOrder[it.second->begin()->first] = tmpOrder;
1727  }
1728 
1729  // See if we already have either region1 or region2 stored in
1730  // perComps map already and do a sanity check to ensure regions
1731  // are mutually periodic.
1732  if (perComps.count(cId1) == 0)
1733  {
1734  if (perComps.count(cId2) == 0)
1735  {
1736  perComps[cId1] = cId2;
1737  }
1738  else
1739  {
1740  std::stringstream ss;
1741  ss << "Boundary region " << cId2 << " should be "
1742  << "periodic with " << perComps[cId2] << " but "
1743  << "found " << cId1 << " instead!";
1744  ASSERTL0(perComps[cId2] == cId1, ss.str());
1745  }
1746  }
1747  else
1748  {
1749  std::stringstream ss;
1750  ss << "Boundary region " << cId1 << " should be "
1751  << "periodic with " << perComps[cId1] << " but "
1752  << "found " << cId2 << " instead!";
1753  ASSERTL0(perComps[cId1] == cId1, ss.str());
1754  }
1755  }
1756 
1757  // In case of periodic partition being split over many composites
1758  // may not have all periodic matches so check this
1759  int idmax = -1;
1760  for (auto &cIt : perComps)
1761  {
1762  idmax = max(idmax, cIt.first);
1763  idmax = max(idmax, cIt.second);
1764  }
1765  vComm->AllReduce(idmax, LibUtilities::ReduceMax);
1766  idmax++;
1767  Array<OneD, int> perid(idmax, -1);
1768  for (auto &cIt : perComps)
1769  {
1770  perid[cIt.first] = cIt.second;
1771  }
1772  vComm->AllReduce(perid, LibUtilities::ReduceMax);
1773  // update all partitions
1774  for (int i = 0; i < idmax; ++i)
1775  {
1776  if (perid[i] > -1)
1777  {
1778  // skip if equivlaent relationship has
1779  // already been speficied in list
1780  if (perComps.count(perid[i]))
1781  {
1782  continue;
1783  }
1784  perComps[i] = perid[i];
1785  }
1786  }
1787 
1788  // The next routines process local face lists to
1789  // exchange vertices,
1790  // edges and faces.
1791  int n = vComm->GetSize();
1792  int p = vComm->GetRank();
1793  int totFaces;
1794  Array<OneD, int> facecounts(n, 0);
1795  Array<OneD, int> vertcounts(n, 0);
1796  Array<OneD, int> faceoffset(n, 0);
1797  Array<OneD, int> vertoffset(n, 0);
1798 
1799  Array<OneD, int> rotcounts(n, 0);
1800  Array<OneD, int> rotoffset(n, 0);
1801 
1802  rotcounts[p] = rotComp.size();
1803  vComm->AllReduce(rotcounts, LibUtilities::ReduceSum);
1804  int totrot = Vmath::Vsum(n, rotcounts, 1);
1805 
1806  if (totrot)
1807  {
1808  for (i = 1; i < n; ++i)
1809  {
1810  rotoffset[i] = rotoffset[i - 1] + rotcounts[i - 1];
1811  }
1812 
1813  Array<OneD, int> compid(totrot, 0);
1814  Array<OneD, int> rotdir(totrot, 0);
1815  Array<OneD, NekDouble> rotangle(totrot, 0.0);
1816  Array<OneD, NekDouble> rottol(totrot, 0.0);
1817 
1818  // fill in rotational informaiton
1819  auto rIt = rotComp.begin();
1820 
1821  for (i = 0; rIt != rotComp.end(); ++rIt)
1822  {
1823  compid[rotoffset[p] + i] = rIt->first;
1824  rotdir[rotoffset[p] + i] = rIt->second.m_dir;
1825  rotangle[rotoffset[p] + i] = rIt->second.m_angle;
1826  rottol[rotoffset[p] + i++] = rIt->second.m_tol;
1827  }
1828 
1829  vComm->AllReduce(compid, LibUtilities::ReduceSum);
1830  vComm->AllReduce(rotdir, LibUtilities::ReduceSum);
1831  vComm->AllReduce(rotangle, LibUtilities::ReduceSum);
1832  vComm->AllReduce(rottol, LibUtilities::ReduceSum);
1833 
1834  // Fill in full rotational composite list
1835  for (i = 0; i < totrot; ++i)
1836  {
1837  RotPeriodicInfo rinfo(rotdir[i], rotangle[i], rottol[i]);
1838 
1839  rotComp[compid[i]] = rinfo;
1840  }
1841  }
1842 
1843  // First exchange the number of faces on each process.
1844  facecounts[p] = locFaces.size();
1845  vComm->AllReduce(facecounts, LibUtilities::ReduceSum);
1846 
1847  // Set up an offset map to allow us to distribute face IDs to all
1848  // processors.
1849  faceoffset[0] = 0;
1850  for (i = 1; i < n; ++i)
1851  {
1852  faceoffset[i] = faceoffset[i - 1] + facecounts[i - 1];
1853  }
1854 
1855  // Calculate total number of faces.
1856  totFaces = Vmath::Vsum(n, facecounts, 1);
1857 
1858  // faceIds holds face IDs for each periodic face. faceVerts holds
1859  // the number of vertices in this face.
1860  Array<OneD, int> faceIds(totFaces, 0);
1861  Array<OneD, int> faceVerts(totFaces, 0);
1862 
1863  // Process p writes IDs of its faces into position faceoffset[p] of
1864  // faceIds which allows us to perform an AllReduce to distribute
1865  // information amongst processors.
1866  auto sIt = locFaces.begin();
1867  for (i = 0; sIt != locFaces.end(); ++sIt)
1868  {
1869  faceIds[faceoffset[p] + i] = *sIt;
1870  faceVerts[faceoffset[p] + i++] = allVerts[*sIt].size();
1871  }
1872 
1873  vComm->AllReduce(faceIds, LibUtilities::ReduceSum);
1874  vComm->AllReduce(faceVerts, LibUtilities::ReduceSum);
1875 
1876  // procVerts holds number of vertices (and also edges since each
1877  // face is 2D) on each process.
1878  Array<OneD, int> procVerts(n, 0);
1879  int nTotVerts;
1880 
1881  // Note if there are no periodic faces at all calling Vsum will
1882  // cause a segfault.
1883  if (totFaces > 0)
1884  {
1885  // Calculate number of vertices on each processor.
1886  nTotVerts = Vmath::Vsum(totFaces, faceVerts, 1);
1887  }
1888  else
1889  {
1890  nTotVerts = 0;
1891  }
1892 
1893  for (i = 0; i < n; ++i)
1894  {
1895  if (facecounts[i] > 0)
1896  {
1897  procVerts[i] = Vmath::Vsum(facecounts[i],
1898  faceVerts + faceoffset[i], 1);
1899  }
1900  else
1901  {
1902  procVerts[i] = 0;
1903  }
1904  }
1905 
1906  // vertoffset is defined in the same manner as edgeoffset
1907  // beforehand.
1908  vertoffset[0] = 0;
1909  for (i = 1; i < n; ++i)
1910  {
1911  vertoffset[i] = vertoffset[i - 1] + procVerts[i - 1];
1912  }
1913 
1914  // At this point we exchange all vertex IDs, edge IDs and vertex
1915  // coordinates for each face. The coordinates are necessary because
1916  // we need to calculate relative face orientations between periodic
1917  // faces to determined edge and vertex connectivity.
1918  Array<OneD, int> vertIds(nTotVerts, 0);
1919  Array<OneD, int> edgeIds(nTotVerts, 0);
1920  Array<OneD, int> edgeOrt(nTotVerts, 0);
1921  Array<OneD, NekDouble> vertX(nTotVerts, 0.0);
1922  Array<OneD, NekDouble> vertY(nTotVerts, 0.0);
1923  Array<OneD, NekDouble> vertZ(nTotVerts, 0.0);
1924 
1925  for (cnt = 0, sIt = locFaces.begin(); sIt != locFaces.end(); ++sIt)
1926  {
1927  for (j = 0; j < allVerts[*sIt].size(); ++j)
1928  {
1929  int vertId = allVerts[*sIt][j];
1930  vertIds[vertoffset[p] + cnt] = vertId;
1931  vertX[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(0);
1932  vertY[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(1);
1933  vertZ[vertoffset[p] + cnt] = (*allCoord[*sIt][j])(2);
1934  edgeIds[vertoffset[p] + cnt] = allEdges[*sIt][j];
1935  edgeOrt[vertoffset[p] + cnt++] = allOrient[*sIt][j];
1936  }
1937  }
1938 
1939  vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
1940  vComm->AllReduce(vertX, LibUtilities::ReduceSum);
1941  vComm->AllReduce(vertY, LibUtilities::ReduceSum);
1942  vComm->AllReduce(vertZ, LibUtilities::ReduceSum);
1943  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
1944  vComm->AllReduce(edgeOrt, LibUtilities::ReduceSum);
1945 
1946  // Finally now we have all of this information, we construct maps
1947  // which make accessing the information easier. These are
1948  // conceptually the same as all* maps at the beginning of the
1949  // routine, but now hold information for all periodic vertices.
1950  map<int, vector<int>> vertMap;
1951  map<int, vector<int>> edgeMap;
1952  map<int, SpatialDomains::PointGeomVector> coordMap;
1953 
1954  // These final two maps are required for determining the relative
1955  // orientation of periodic edges. vCoMap associates vertex IDs with
1956  // their coordinates, and eIdMap maps an edge ID to the two vertices
1957  // which construct it.
1958  map<int, SpatialDomains::PointGeomSharedPtr> vCoMap;
1959  map<int, pair<int, int>> eIdMap;
1960 
1961  for (cnt = i = 0; i < totFaces; ++i)
1962  {
1963  vector<int> edges(faceVerts[i]);
1964  vector<int> verts(faceVerts[i]);
1965  SpatialDomains::PointGeomVector coord(faceVerts[i]);
1966 
1967  // Keep track of cnt to enable correct edge vertices to be
1968  // inserted into eIdMap.
1969  int tmp = cnt;
1970  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
1971  {
1972  edges[j] = edgeIds[cnt];
1973  verts[j] = vertIds[cnt];
1975  AllocateSharedPtr(3, verts[j], vertX[cnt], vertY[cnt],
1976  vertZ[cnt]);
1977  vCoMap[vertIds[cnt]] = coord[j];
1978 
1979  // Try to insert edge into the eIdMap to avoid re-inserting.
1980  auto testIns = eIdMap.insert(make_pair(
1981  edgeIds[cnt],
1982  make_pair(vertIds[tmp + j],
1983  vertIds[tmp + ((j + 1) % faceVerts[i])])));
1984 
1985  if (testIns.second == false)
1986  {
1987  continue;
1988  }
1989 
1990  // If the edge is reversed with respect to the face, then
1991  // swap the edges so that we have the original ordering of
1992  // the edge in the 3D element. This is necessary to properly
1993  // determine edge orientation. Note that the logic relies on
1994  // the fact that all edge forward directions are CCW
1995  // orientated: we use a tensor product ordering for 2D
1996  // elements so need to reverse this for edge IDs 2 and 3.
1997  StdRegions::Orientation edgeOrient =
1998  static_cast<StdRegions::Orientation>(edgeOrt[cnt]);
1999  if (j > 1)
2000  {
2001  edgeOrient = edgeOrient == StdRegions::eForwards
2004  }
2005 
2006  if (edgeOrient == StdRegions::eBackwards)
2007  {
2008  swap(testIns.first->second.first,
2009  testIns.first->second.second);
2010  }
2011  }
2012 
2013  vertMap[faceIds[i]] = verts;
2014  edgeMap[faceIds[i]] = edges;
2015  coordMap[faceIds[i]] = coord;
2016  }
2017 
2018  // Go through list of composites and figure out which edges are
2019  // parallel from original ordering in session file. This includes
2020  // composites which are not necessarily on this process.
2021 
2022  // Store temporary map of periodic vertices which will hold all
2023  // periodic vertices on the entire mesh so that doubly periodic
2024  // vertices/edges can be counted properly across partitions. Local
2025  // vertices/edges are copied into m_periodicVerts and
2026  // m_periodicEdges at the end of the function.
2027  PeriodicMap periodicVerts, periodicEdges;
2028 
2029  // Construct two maps which determine how vertices and edges of
2030  // faces connect given a specific face orientation. The key of the
2031  // map is the number of vertices in the face, used to determine
2032  // difference between tris and quads.
2033  map<int, map<StdRegions::Orientation, vector<int>>> vmap;
2034  map<int, map<StdRegions::Orientation, vector<int>>> emap;
2035 
2036  map<StdRegions::Orientation, vector<int>> quadVertMap;
2037  quadVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
2038  quadVertMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {3, 2, 1, 0};
2039  quadVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 3, 2};
2040  quadVertMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
2041  quadVertMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {0, 3, 2, 1};
2042  quadVertMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
2043  quadVertMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
2044  quadVertMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {2, 1, 0, 3};
2045 
2046  map<StdRegions::Orientation, vector<int>> quadEdgeMap;
2047  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2, 3};
2048  quadEdgeMap[StdRegions::eDir1FwdDir1_Dir2BwdDir2] = {2, 1, 0, 3};
2049  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 3, 2, 1};
2050  quadEdgeMap[StdRegions::eDir1BwdDir1_Dir2BwdDir2] = {2, 3, 0, 1};
2051  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2FwdDir1] = {3, 2, 1, 0};
2052  quadEdgeMap[StdRegions::eDir1FwdDir2_Dir2BwdDir1] = {1, 2, 3, 0};
2053  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2FwdDir1] = {3, 0, 1, 2};
2054  quadEdgeMap[StdRegions::eDir1BwdDir2_Dir2BwdDir1] = {1, 0, 3, 2};
2055 
2056  map<StdRegions::Orientation, vector<int>> triVertMap;
2057  triVertMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
2058  triVertMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {1, 0, 2};
2059 
2060  map<StdRegions::Orientation, vector<int>> triEdgeMap;
2061  triEdgeMap[StdRegions::eDir1FwdDir1_Dir2FwdDir2] = {0, 1, 2};
2062  triEdgeMap[StdRegions::eDir1BwdDir1_Dir2FwdDir2] = {0, 2, 1};
2063 
2064  vmap[3] = triVertMap;
2065  vmap[4] = quadVertMap;
2066  emap[3] = triEdgeMap;
2067  emap[4] = quadEdgeMap;
2068 
2069  map<int, int> allCompPairs;
2070 
2071  // Collect composite id's of each periodic face for use if rotation
2072  // is required
2073  map<int, int> fIdToCompId;
2074 
2075  // Finally we have enough information to populate the periodic
2076  // vertex, edge and face maps. Begin by looping over all pairs of
2077  // periodic composites to determine pairs of periodic faces.
2078  for (auto &cIt : perComps)
2079  {
2081  const int id1 = cIt.first;
2082  const int id2 = cIt.second;
2083  std::string id1s = boost::lexical_cast<string>(id1);
2084  std::string id2s = boost::lexical_cast<string>(id2);
2085 
2086  if (compMap.count(id1) > 0)
2087  {
2088  c[0] = compMap[id1];
2089  }
2090 
2091  if (compMap.count(id2) > 0)
2092  {
2093  c[1] = compMap[id2];
2094  }
2095 
2096  // Loop over composite ordering to construct list of all
2097  // periodic faces, regardless of whether they are on this
2098  // process.
2099  map<int, int> compPairs;
2100 
2101  ASSERTL0(compOrder.count(id1) > 0,
2102  "Unable to find composite " + id1s + " in order map.");
2103  ASSERTL0(compOrder.count(id2) > 0,
2104  "Unable to find composite " + id2s + " in order map.");
2105  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
2106  "Periodic composites " + id1s + " and " + id2s +
2107  " should have the same number of elements.");
2108  ASSERTL0(compOrder[id1].size() > 0, "Periodic composites " +
2109  id1s + " and " + id2s +
2110  " are empty!");
2111 
2112  // Look up composite ordering to determine pairs.
2113  for (i = 0; i < compOrder[id1].size(); ++i)
2114  {
2115  int eId1 = compOrder[id1][i];
2116  int eId2 = compOrder[id2][i];
2117 
2118  ASSERTL0(compPairs.count(eId1) == 0, "Already paired.");
2119 
2120  // Sanity check that the faces are mutually periodic.
2121  if (compPairs.count(eId2) != 0)
2122  {
2123  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
2124  }
2125  compPairs[eId1] = eId2;
2126 
2127  // store a map of face ids to composite ids
2128  fIdToCompId[eId1] = id1;
2129  fIdToCompId[eId2] = id2;
2130  }
2131 
2132  // Now that we have all pairs of periodic faces, loop over the
2133  // ones local on this process and populate face/edge/vertex
2134  // maps.
2135  for (auto &pIt : compPairs)
2136  {
2137  int ids[2] = {pIt.first, pIt.second};
2138  bool local[2] = {locFaces.count(pIt.first) > 0,
2139  locFaces.count(pIt.second) > 0};
2140 
2141  ASSERTL0(coordMap.count(ids[0]) > 0 &&
2142  coordMap.count(ids[1]) > 0,
2143  "Unable to find face in coordinate map");
2144 
2145  allCompPairs[pIt.first] = pIt.second;
2146  allCompPairs[pIt.second] = pIt.first;
2147 
2148  // Loop up coordinates of the faces, check they have the
2149  // same number of vertices.
2150  SpatialDomains::PointGeomVector tmpVec[2] = {
2151  coordMap[ids[0]], coordMap[ids[1]]};
2152 
2153  ASSERTL0(tmpVec[0].size() == tmpVec[1].size(),
2154  "Two periodic faces have different number "
2155  "of vertices!");
2156 
2157  // o will store relative orientation of faces. Note that in
2158  // some transpose cases (Dir1FwdDir2_Dir2BwdDir1 and
2159  // Dir1BwdDir1_Dir2FwdDir1) it seems orientation will be
2160  // different going from face1->face2 instead of face2->face1
2161  // (check this).
2163  bool rotbnd = false;
2164  int dir = 0;
2165  NekDouble angle = 0.0;
2166  NekDouble sign = 0.0;
2167  NekDouble tol = 1e-8;
2168 
2169  // check to see if perioid boundary is rotated
2170  if (rotComp.count(fIdToCompId[pIt.first]))
2171  {
2172  rotbnd = true;
2173  dir = rotComp[fIdToCompId[pIt.first]].m_dir;
2174  angle = rotComp[fIdToCompId[pIt.first]].m_angle;
2175  tol = rotComp[fIdToCompId[pIt.first]].m_tol;
2176  }
2177 
2178  // Record periodic faces.
2179  for (i = 0; i < 2; ++i)
2180  {
2181  if (!local[i])
2182  {
2183  continue;
2184  }
2185 
2186  // Reference to the other face.
2187  int other = (i + 1) % 2;
2188 
2189  // angle is set up for i = 0 to i = 1
2190  sign = (i == 0) ? 1.0 : -1.0;
2191 
2192  // Calculate relative face orientation.
2193  if (tmpVec[0].size() == 3)
2194  {
2196  tmpVec[i], tmpVec[other], rotbnd, dir,
2197  sign * angle, tol);
2198  }
2199  else
2200  {
2202  tmpVec[i], tmpVec[other], rotbnd, dir,
2203  sign * angle, tol);
2204  }
2205 
2206  // Record face ID, orientation and whether other face is
2207  // local.
2208  PeriodicEntity ent(ids[other], o, local[other]);
2209  m_periodicFaces[ids[i]].push_back(ent);
2210  }
2211 
2212  int nFaceVerts = vertMap[ids[0]].size();
2213 
2214  // Determine periodic vertices.
2215  for (i = 0; i < 2; ++i)
2216  {
2217  int other = (i + 1) % 2;
2218 
2219  // angle is set up for i = 0 to i = 1
2220  sign = (i == 0) ? 1.0 : -1.0;
2221 
2222  // Calculate relative face orientation.
2223  if (tmpVec[0].size() == 3)
2224  {
2226  tmpVec[i], tmpVec[other], rotbnd, dir,
2227  sign * angle, tol);
2228  }
2229  else
2230  {
2232  tmpVec[i], tmpVec[other], rotbnd, dir,
2233  sign * angle, tol);
2234  }
2235 
2236  if (nFaceVerts == 3)
2237  {
2238  ASSERTL0(
2241  "Unsupported face orientation for face " +
2242  boost::lexical_cast<string>(ids[i]));
2243  }
2244 
2245  // Look up vertices for this face.
2246  vector<int> per1 = vertMap[ids[i]];
2247  vector<int> per2 = vertMap[ids[other]];
2248 
2249  // tmpMap will hold the pairs of vertices which are
2250  // periodic.
2251  map<int, pair<int, bool>> tmpMap;
2252 
2253  // Use vmap to determine which vertices connect given
2254  // the orientation o.
2255  for (j = 0; j < nFaceVerts; ++j)
2256  {
2257  int v = vmap[nFaceVerts][o][j];
2258  tmpMap[per1[j]] =
2259  make_pair(per2[v], locVerts.count(per2[v]) > 0);
2260  }
2261 
2262  // Now loop over tmpMap to associate periodic vertices.
2263  for (auto &mIt : tmpMap)
2264  {
2265  PeriodicEntity ent2(mIt.second.first,
2267  mIt.second.second);
2268 
2269  // See if this vertex has been recorded already.
2270  auto perIt = periodicVerts.find(mIt.first);
2271 
2272  if (perIt == periodicVerts.end())
2273  {
2274  // Vertex is new - just record this entity as
2275  // usual.
2276  periodicVerts[mIt.first].push_back(ent2);
2277  perIt = periodicVerts.find(mIt.first);
2278  }
2279  else
2280  {
2281  // Vertex is known - loop over the vertices
2282  // inside the record and potentially add vertex
2283  // mIt.second to the list.
2284  for (k = 0; k < perIt->second.size(); ++k)
2285  {
2286  if (perIt->second[k].id == mIt.second.first)
2287  {
2288  break;
2289  }
2290  }
2291 
2292  if (k == perIt->second.size())
2293  {
2294  perIt->second.push_back(ent2);
2295  }
2296  }
2297  }
2298  }
2299 
2300  // Determine periodic edges. Logic is the same as above,
2301  // and perhaps should be condensed to avoid replication.
2302  for (i = 0; i < 2; ++i)
2303  {
2304  int other = (i + 1) % 2;
2305 
2306  // angle is set up for i = 0 to i = 1
2307  sign = (i == 0) ? 1.0 : -1.0;
2308 
2309  if (tmpVec[0].size() == 3)
2310  {
2312  tmpVec[i], tmpVec[other], rotbnd, dir,
2313  sign * angle, tol);
2314  }
2315  else
2316  {
2318  tmpVec[i], tmpVec[other], rotbnd, dir,
2319  sign * angle, tol);
2320  }
2321 
2322  vector<int> per1 = edgeMap[ids[i]];
2323  vector<int> per2 = edgeMap[ids[other]];
2324 
2325  map<int, pair<int, bool>> tmpMap;
2326 
2327  for (j = 0; j < nFaceVerts; ++j)
2328  {
2329  int e = emap[nFaceVerts][o][j];
2330  tmpMap[per1[j]] =
2331  make_pair(per2[e], locEdges.count(per2[e]) > 0);
2332  }
2333 
2334  for (auto &mIt : tmpMap)
2335  {
2336  // Note we assume orientation of edges is forwards -
2337  // this may be reversed later.
2338  PeriodicEntity ent2(mIt.second.first,
2340  mIt.second.second);
2341  auto perIt = periodicEdges.find(mIt.first);
2342 
2343  if (perIt == periodicEdges.end())
2344  {
2345  periodicEdges[mIt.first].push_back(ent2);
2346  perIt = periodicEdges.find(mIt.first);
2347  }
2348  else
2349  {
2350  for (k = 0; k < perIt->second.size(); ++k)
2351  {
2352  if (perIt->second[k].id == mIt.second.first)
2353  {
2354  break;
2355  }
2356  }
2357 
2358  if (k == perIt->second.size())
2359  {
2360  perIt->second.push_back(ent2);
2361  }
2362  }
2363  }
2364  }
2365  }
2366  }
2367 
2368  // Make sure that the nubmer of face pairs and the
2369  // face Id to composite Id map match in size
2370  ASSERTL1(allCompPairs.size() == fIdToCompId.size(),
2371  "At this point the size of allCompPairs "
2372  "should have been the same as fIdToCompId");
2373 
2374  // also will need an edge id to composite id at
2375  // end of routine
2376  map<int, int> eIdToCompId;
2377 
2378  // Search for periodic vertices and edges which are not
2379  // in a periodic composite but lie in this process. First,
2380  // loop over all information we have from other
2381  // processors.
2382  for (cnt = i = 0; i < totFaces; ++i)
2383  {
2384  bool rotbnd = false;
2385  int dir = 0;
2386  NekDouble angle = 0.0;
2387  NekDouble tol = 1e-8;
2388 
2389  int faceId = faceIds[i];
2390 
2391  ASSERTL0(allCompPairs.count(faceId) > 0,
2392  "Unable to find matching periodic face. faceId = " +
2393  boost::lexical_cast<string>(faceId));
2394 
2395  int perFaceId = allCompPairs[faceId];
2396 
2397  // check to see if periodic boundary is rotated
2398  ASSERTL1((rotComp.size() == 0) || fIdToCompId.count(faceId) > 0,
2399  "Face " + boost::lexical_cast<string>(faceId) +
2400  " not found in fIdtoCompId map");
2401  if (rotComp.count(fIdToCompId[faceId]))
2402  {
2403  rotbnd = true;
2404  dir = rotComp[fIdToCompId[faceId]].m_dir;
2405  angle = rotComp[fIdToCompId[faceId]].m_angle;
2406  tol = rotComp[fIdToCompId[faceId]].m_tol;
2407  }
2408 
2409  for (j = 0; j < faceVerts[i]; ++j, ++cnt)
2410  {
2411  int vId = vertIds[cnt];
2412 
2413  auto perId = periodicVerts.find(vId);
2414 
2415  if (perId == periodicVerts.end())
2416  {
2417 
2418  // This vertex is not included in the
2419  // map. Figure out which vertex it is supposed
2420  // to be periodic with. perFaceId is the face
2421  // ID which is periodic with faceId. The logic
2422  // is much the same as the loop above.
2423  SpatialDomains::PointGeomVector tmpVec[2] = {
2424  coordMap[faceId], coordMap[perFaceId]};
2425 
2426  int nFaceVerts = tmpVec[0].size();
2428  nFaceVerts == 3
2430  tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2431  tol)
2433  tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2434  tol);
2435 
2436  // Use vmap to determine which vertex of the other face
2437  // should be periodic with this one.
2438  int perVertexId =
2439  vertMap[perFaceId][vmap[nFaceVerts][o][j]];
2440 
2441  PeriodicEntity ent(perVertexId,
2443  locVerts.count(perVertexId) > 0);
2444 
2445  periodicVerts[vId].push_back(ent);
2446  }
2447 
2448  int eId = edgeIds[cnt];
2449 
2450  perId = periodicEdges.find(eId);
2451 
2452  // this map is required at very end to determine rotation of
2453  // edges.
2454  if (rotbnd)
2455  {
2456  eIdToCompId[eId] = fIdToCompId[faceId];
2457  }
2458 
2459  if (perId == periodicEdges.end())
2460  {
2461  // This edge is not included in the map. Figure
2462  // out which edge it is supposed to be periodic
2463  // with. perFaceId is the face ID which is
2464  // periodic with faceId. The logic is much the
2465  // same as the loop above.
2466  SpatialDomains::PointGeomVector tmpVec[2] = {
2467  coordMap[faceId], coordMap[perFaceId]};
2468 
2469  int nFaceEdges = tmpVec[0].size();
2471  nFaceEdges == 3
2473  tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2474  tol)
2476  tmpVec[0], tmpVec[1], rotbnd, dir, angle,
2477  tol);
2478 
2479  // Use emap to determine which edge of the other
2480  // face should be periodic with this one.
2481  int perEdgeId =
2482  edgeMap[perFaceId][emap[nFaceEdges][o][j]];
2483 
2484  PeriodicEntity ent(perEdgeId, StdRegions::eForwards,
2485  locEdges.count(perEdgeId) > 0);
2486 
2487  periodicEdges[eId].push_back(ent);
2488 
2489  // this map is required at very end to
2490  // determine rotation of edges.
2491  if (rotbnd)
2492  {
2493  eIdToCompId[perEdgeId] = fIdToCompId[perFaceId];
2494  }
2495  }
2496  }
2497  }
2498 
2499  // Finally, we must loop over the periodicVerts and periodicEdges
2500  // map to complete connectivity information.
2501  for (auto &perIt : periodicVerts)
2502  {
2503  // For each vertex that is periodic with this one...
2504  for (i = 0; i < perIt.second.size(); ++i)
2505  {
2506  // Find the vertex in the periodicVerts map...
2507  auto perIt2 = periodicVerts.find(perIt.second[i].id);
2508  ASSERTL0(perIt2 != periodicVerts.end(),
2509  "Couldn't find periodic vertex.");
2510 
2511  // Now search through this vertex's list and make sure that
2512  // we have a record of any vertices which aren't in the
2513  // original list.
2514  for (j = 0; j < perIt2->second.size(); ++j)
2515  {
2516  if (perIt2->second[j].id == perIt.first)
2517  {
2518  continue;
2519  }
2520 
2521  for (k = 0; k < perIt.second.size(); ++k)
2522  {
2523  if (perIt2->second[j].id == perIt.second[k].id)
2524  {
2525  break;
2526  }
2527  }
2528 
2529  if (k == perIt.second.size())
2530  {
2531  perIt.second.push_back(perIt2->second[j]);
2532  }
2533  }
2534  }
2535  }
2536 
2537  for (auto &perIt : periodicEdges)
2538  {
2539  for (i = 0; i < perIt.second.size(); ++i)
2540  {
2541  auto perIt2 = periodicEdges.find(perIt.second[i].id);
2542  ASSERTL0(perIt2 != periodicEdges.end(),
2543  "Couldn't find periodic edge.");
2544 
2545  for (j = 0; j < perIt2->second.size(); ++j)
2546  {
2547  if (perIt2->second[j].id == perIt.first)
2548  {
2549  continue;
2550  }
2551 
2552  for (k = 0; k < perIt.second.size(); ++k)
2553  {
2554  if (perIt2->second[j].id == perIt.second[k].id)
2555  {
2556  break;
2557  }
2558  }
2559 
2560  if (k == perIt.second.size())
2561  {
2562  perIt.second.push_back(perIt2->second[j]);
2563  }
2564  }
2565  }
2566  }
2567 
2568  // Loop over periodic edges to determine relative edge orientations.
2569  for (auto &perIt : periodicEdges)
2570  {
2571  bool rotbnd = false;
2572  int dir = 0;
2573  NekDouble angle = 0.0;
2574  NekDouble tol = 1e-8;
2575 
2576  // Find edge coordinates
2577  auto eIt = eIdMap.find(perIt.first);
2578  SpatialDomains::PointGeom v[2] = {*vCoMap[eIt->second.first],
2579  *vCoMap[eIt->second.second]};
2580 
2581  // check to see if perioid boundary is rotated
2582  if (rotComp.count(eIdToCompId[perIt.first]))
2583  {
2584  rotbnd = true;
2585  dir = rotComp[eIdToCompId[perIt.first]].m_dir;
2586  angle = rotComp[eIdToCompId[perIt.first]].m_angle;
2587  tol = rotComp[eIdToCompId[perIt.first]].m_tol;
2588  }
2589 
2590  // Loop over each edge, and construct a vector that takes us
2591  // from one vertex to another. Use this to figure out which
2592  // vertex maps to which.
2593  for (i = 0; i < perIt.second.size(); ++i)
2594  {
2595  eIt = eIdMap.find(perIt.second[i].id);
2596 
2597  SpatialDomains::PointGeom w[2] = {
2598  *vCoMap[eIt->second.first],
2599  *vCoMap[eIt->second.second]};
2600 
2601  int vMap[2] = {-1, -1};
2602  if (rotbnd)
2603  {
2604 
2606 
2607  r.Rotate(v[0], dir, angle);
2608 
2609  if (r.dist(w[0]) < tol)
2610  {
2611  vMap[0] = 0;
2612  }
2613  else
2614  {
2615  r.Rotate(v[1], dir, angle);
2616  if (r.dist(w[0]) < tol)
2617  {
2618  vMap[0] = 1;
2619  }
2620  else
2621  {
2623  "Unable to align rotationally "
2624  "periodic edge vertex");
2625  }
2626  }
2627  }
2628  else // translation test
2629  {
2630  NekDouble cx =
2631  0.5 * (w[0](0) - v[0](0) + w[1](0) - v[1](0));
2632  NekDouble cy =
2633  0.5 * (w[0](1) - v[0](1) + w[1](1) - v[1](1));
2634  NekDouble cz =
2635  0.5 * (w[0](2) - v[0](2) + w[1](2) - v[1](2));
2636 
2637  for (j = 0; j < 2; ++j)
2638  {
2639  NekDouble x = v[j](0);
2640  NekDouble y = v[j](1);
2641  NekDouble z = v[j](2);
2642  for (k = 0; k < 2; ++k)
2643  {
2644  NekDouble x1 = w[k](0) - cx;
2645  NekDouble y1 = w[k](1) - cy;
2646  NekDouble z1 = w[k](2) - cz;
2647 
2648  if (sqrt((x1 - x) * (x1 - x) +
2649  (y1 - y) * (y1 - y) +
2650  (z1 - z) * (z1 - z)) < 1e-8)
2651  {
2652  vMap[k] = j;
2653  break;
2654  }
2655  }
2656  }
2657 
2658  // Sanity check the map.
2659  ASSERTL0(vMap[0] >= 0 && vMap[1] >= 0,
2660  "Unable to align periodic edge vertex.");
2661  ASSERTL0((vMap[0] == 0 || vMap[0] == 1) &&
2662  (vMap[1] == 0 || vMap[1] == 1) &&
2663  (vMap[0] != vMap[1]),
2664  "Unable to align periodic edge vertex.");
2665  }
2666 
2667  // If 0 -> 0 then edges are aligned already; otherwise
2668  // reverse the orientation.
2669  if (vMap[0] != 0)
2670  {
2671  perIt.second[i].orient = StdRegions::eBackwards;
2672  }
2673  }
2674  }
2675 
2676  // Do one final loop over periodic vertices/edges to remove
2677  // non-local vertices/edges from map.
2678  for (auto &perIt : periodicVerts)
2679  {
2680  if (locVerts.count(perIt.first) > 0)
2681  {
2682  m_periodicVerts.insert(perIt);
2683  }
2684  }
2685 
2686  for (auto &perIt : periodicEdges)
2687  {
2688  if (locEdges.count(perIt.first) > 0)
2689  {
2690  m_periodicEdges.insert(perIt);
2691  }
2692  }
2693  }
2694  break;
2695  default:
2696  ASSERTL1(false, "Not setup for this expansion");
2697  break;
2698  }
2699 }
2700 
2701 /**
2702  *
2703  */
2705  const GlobalLinSysKey &mkey)
2706 {
2708  "Routine currently only tested for HybridDGHelmholtz");
2709 
2711  "Full matrix global systems are not supported for HDG "
2712  "expansions");
2713 
2714  ASSERTL1(mkey.GetGlobalSysSolnType() == m_traceMap->GetGlobalSysSolnType(),
2715  "The local to global map is not set up for the requested "
2716  "solution type");
2717 
2718  GlobalLinSysSharedPtr glo_matrix;
2719  auto matrixIter = m_globalBndMat->find(mkey);
2720 
2721  if (matrixIter == m_globalBndMat->end())
2722  {
2723  glo_matrix = GenGlobalBndLinSys(mkey, m_traceMap);
2724  (*m_globalBndMat)[mkey] = glo_matrix;
2725  }
2726  else
2727  {
2728  glo_matrix = matrixIter->second;
2729  }
2730 
2731  return glo_matrix;
2732 }
2733 
2734 /**
2735  * For each boundary region, checks that the types and number of
2736  * boundary expansions in that region match.
2737  * @param In Field to compare with.
2738  * @return True if boundary conditions match.
2739  */
2741 {
2742  int i;
2743  bool returnval = true;
2744 
2745  for (i = 0; i < m_bndConditions.size(); ++i)
2746  {
2747  // check to see if boundary condition type is the same
2748  // and there are the same number of boundary
2749  // conditions in the boundary definition.
2750  if ((m_bndConditions[i]->GetBoundaryConditionType() !=
2751  In.m_bndConditions[i]->GetBoundaryConditionType()) ||
2752  (m_bndCondExpansions[i]->GetExpSize() !=
2753  In.m_bndCondExpansions[i]->GetExpSize()))
2754  {
2755  returnval = false;
2756  break;
2757  }
2758  }
2759 
2760  // Compare with all other processes. Return true only if all
2761  // processes report having the same boundary conditions.
2762  int vSame = (returnval ? 1 : 0);
2763  m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
2764 
2765  return (vSame == 1);
2766 }
2767 
2769 {
2770  if (m_negatedFluxNormal.size() == 0)
2771  {
2773  &elmtToTrace = m_traceMap->GetElmtToTrace();
2774 
2775  m_negatedFluxNormal.resize(2 * GetExpSize());
2776 
2777  for (int i = 0; i < GetExpSize(); ++i)
2778  {
2779 
2780  for (int v = 0; v < 2; ++v)
2781  {
2783  elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
2784 
2785  if (vertExp->GetLeftAdjacentElementExp()
2786  ->GetGeom()
2787  ->GetGlobalID() !=
2788  (*m_exp)[i]->GetGeom()->GetGlobalID())
2789  {
2790  m_negatedFluxNormal[2 * i + v] = true;
2791  }
2792  else
2793  {
2794  m_negatedFluxNormal[2 * i + v] = false;
2795  }
2796  }
2797  }
2798  }
2799 
2800  return m_negatedFluxNormal;
2801 }
2802 
2803 /**
2804  * \brief This method extracts the "forward" and "backward" trace
2805  * data from the array \a field and puts the data into output
2806  * vectors \a Fwd and \a Bwd.
2807  *
2808  * We first define the convention which defines "forwards" and
2809  * "backwards". First an association is made between the vertex/edge/face of
2810  * each element and its corresponding vertex/edge/face in the trace space
2811  * using the mapping #m_traceMap. The element can either be
2812  * left-adjacent or right-adjacent to this trace face (see
2813  * Expansion2D::GetLeftAdjacentElementExp). Boundary faces are
2814  * always left-adjacent since left-adjacency is populated first.
2815  *
2816  * If the element is left-adjacent we extract the trace data
2817  * from \a field into the forward trace space \a Fwd; otherwise,
2818  * we place it in the backwards trace space \a Bwd. In this way,
2819  * we form a unique set of trace normals since these are always
2820  * extracted from left-adjacent elements.
2821  *
2822  * \param field is a NekDouble array which contains the fielddata
2823  * from which we wish to extract the backward and forward
2824  * orientated trace/face arrays.
2825  *
2826  * \return Updates a NekDouble array \a Fwd and \a Bwd
2827  */
2831  bool FillBnd, // these should be template params so that compiler can remove
2832  // them
2833  bool PutFwdInBwdOnBCs, bool DoExchange)
2834 {
2835  // Is this zeroing necessary?
2836  // Zero forward/backward vectors.
2837  Vmath::Zero(Fwd.size(), Fwd, 1);
2838  Vmath::Zero(Bwd.size(), Bwd, 1);
2839 
2840  // Basis definition on each element
2841  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
2842  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
2843  {
2844  // blocked routine
2845  Array<OneD, NekDouble> tracevals(
2846  m_locTraceToTraceMap->GetNLocTracePts());
2847 
2848  m_locTraceToTraceMap->LocTracesFromField(field, tracevals);
2849  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, Fwd);
2850 
2851  Array<OneD, NekDouble> invals =
2852  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
2853  m_locTraceToTraceMap->InterpLocTracesToTrace(1, invals, Bwd);
2854  }
2855  else
2856  {
2857  // Loop over elements and collect forward expansion
2858  auto nexp = GetExpSize();
2859  Array<OneD, NekDouble> e_tmp;
2861 
2863  &elmtToTrace = m_traceMap->GetElmtToTrace();
2864 
2865  for (int n = 0, cnt = 0; n < nexp; ++n)
2866  {
2867  exp = (*m_exp)[n];
2868  auto phys_offset = GetPhys_Offset(n);
2869 
2870  for (int e = 0; e < exp->GetNtraces(); ++e, ++cnt)
2871  {
2872  auto offset =
2873  m_trace->GetPhys_Offset(elmtToTrace[n][e]->GetElmtId());
2874 
2875  e_tmp =
2876  (m_leftAdjacentTraces[cnt]) ? Fwd + offset : Bwd + offset;
2877 
2878  exp->GetTracePhysVals(e, elmtToTrace[n][e], field + phys_offset,
2879  e_tmp);
2880  }
2881  }
2882  }
2883 
2885 
2886  if (FillBnd)
2887  {
2888  v_FillBwdWithBoundCond(Fwd, Bwd, PutFwdInBwdOnBCs);
2889  }
2890 
2891  if (DoExchange)
2892  {
2893  // Do parallel exchange for forwards/backwards spaces.
2894  m_traceMap->GetAssemblyCommDG()->PerformExchange(Fwd, Bwd);
2895  }
2896 }
2897 
2900  bool PutFwdInBwdOnBCs)
2901 {
2902  // Fill boundary conditions into missing elements
2903  if (PutFwdInBwdOnBCs) // just set Bwd value to be Fwd value on BCs
2904  {
2905  // Fill boundary conditions into missing elements
2906  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2907  {
2908  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2910  {
2911  auto ne = m_bndCondExpansions[n]->GetExpSize();
2912  for (int e = 0; e < ne; ++e)
2913  {
2914  auto npts =
2915  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2916  auto id2 = m_trace->GetPhys_Offset(
2917  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2918  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2919  }
2920 
2921  cnt += ne;
2922  }
2923  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2925  m_bndConditions[n]->GetBoundaryConditionType() ==
2927  {
2928  auto ne = m_bndCondExpansions[n]->GetExpSize();
2929  for (int e = 0; e < ne; ++e)
2930  {
2931  auto npts =
2932  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2933  auto id2 = m_trace->GetPhys_Offset(
2934  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2935  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2936  }
2937  cnt += ne;
2938  }
2939  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
2941  {
2942  ASSERTL0(false, "Method not set up for this "
2943  "boundary condition.");
2944  }
2945  }
2946  }
2947  else
2948  {
2949  for (int n = 0, cnt = 0; n < m_bndCondExpansions.size(); ++n)
2950  {
2951  if (m_bndConditions[n]->GetBoundaryConditionType() ==
2953  {
2954  auto ne = m_bndCondExpansions[n]->GetExpSize();
2955  for (int e = 0; e < ne; ++e)
2956  {
2957  auto npts =
2958  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2959  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
2960  auto id2 = m_trace->GetPhys_Offset(
2961  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2962  Vmath::Vcopy(npts,
2963  &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
2964  &Bwd[id2], 1);
2965  }
2966  cnt += ne;
2967  }
2968  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2970  m_bndConditions[n]->GetBoundaryConditionType() ==
2972  {
2973  auto ne = m_bndCondExpansions[n]->GetExpSize();
2974  for (int e = 0; e < ne; ++e)
2975  {
2976  auto npts =
2977  m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
2978  auto id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
2979  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[id1] == 0.0,
2980  "Method not set up for non-zero "
2981  "Neumann boundary condition");
2982  auto id2 = m_trace->GetPhys_Offset(
2983  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
2984  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
2985  }
2986 
2987  cnt += ne;
2988  }
2989  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
2991  {
2992  // Do nothing
2993  }
2994  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
2996  {
2998  "Method not set up for this boundary "
2999  "condition.");
3000  }
3001  }
3002  }
3003 }
3004 
3006  const Array<OneD, const NekDouble> &Fwd,
3008 {
3009  // Basis definition on each element
3010  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3011  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
3012  {
3013  Array<OneD, NekDouble> tracevals(
3014  m_locTraceToTraceMap->GetNLocTracePts(), 0.0);
3015 
3016  Array<OneD, NekDouble> invals =
3017  tracevals + m_locTraceToTraceMap->GetNFwdLocTracePts();
3018 
3019  m_locTraceToTraceMap->InterpTraceToLocTrace(1, Bwd, invals);
3020 
3021  m_locTraceToTraceMap->InterpTraceToLocTrace(0, Fwd, tracevals);
3022 
3023  m_locTraceToTraceMap->AddLocTracesToField(tracevals, field);
3024  }
3025  else
3026  {
3028  "v_AddTraceQuadPhysToField not coded for eGauss_Lagrange");
3029  }
3030 }
3031 
3033 {
3034  ASSERTL1(m_physState == true, "local physical space is not true ");
3035  v_ExtractTracePhys(m_phys, outarray);
3036 }
3037 
3038 /**
3039  * @brief This method extracts the trace (verts in 1D) from
3040  * the field @a inarray and puts the values in @a outarray.
3041  *
3042  * It assumes the field is C0 continuous so that it can
3043  * overwrite the edge data when visited by the two adjacent
3044  * elements.
3045  *
3046  * @param inarray An array containing the 1D data from which we wish
3047  * to extract the edge data.
3048  * @param outarray The resulting edge information.
3049  *
3050  * This will not work for non-boundary expansions
3051  */
3053  const Array<OneD, const NekDouble> &inarray,
3054  Array<OneD, NekDouble> &outarray)
3055 {
3056  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3057  if ((basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3058  {
3059  Vmath::Zero(outarray.size(), outarray, 1);
3060 
3061  Array<OneD, NekDouble> tracevals(
3062  m_locTraceToTraceMap->GetNFwdLocTracePts());
3063  m_locTraceToTraceMap->FwdLocTracesFromField(inarray, tracevals);
3064  m_locTraceToTraceMap->InterpLocTracesToTrace(0, tracevals, outarray);
3065  m_traceMap->GetAssemblyCommDG()->PerformExchange(outarray, outarray);
3066  }
3067  else
3068  {
3069 
3070  // Loop over elemente and collect forward expansion
3071  int nexp = GetExpSize();
3072  int n, p, offset, phys_offset;
3073  Array<OneD, NekDouble> t_tmp;
3074 
3076  &elmtToTrace = m_traceMap->GetElmtToTrace();
3077 
3078  ASSERTL1(outarray.size() >= m_trace->GetNpoints(),
3079  "input array is of insufficient length");
3080 
3081  for (n = 0; n < nexp; ++n)
3082  {
3083  phys_offset = GetPhys_Offset(n);
3084 
3085  for (p = 0; p < (*m_exp)[n]->GetNtraces(); ++p)
3086  {
3087  offset =
3088  m_trace->GetPhys_Offset(elmtToTrace[n][p]->GetElmtId());
3089  (*m_exp)[n]->GetTracePhysVals(p, elmtToTrace[n][p],
3090  inarray + phys_offset,
3091  t_tmp = outarray + offset);
3092  }
3093  }
3094  }
3095 }
3096 
3097 /**
3098  * @brief Add trace contributions into elemental coefficient spaces.
3099  *
3100  * Given some quantity \f$ \vec{Fn} \f$, which conatins this
3101  * routine calculates the integral
3102  *
3103  * \f[
3104  * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
3105  * \f]
3106  *
3107  * and adds this to the coefficient space provided by outarray.
3108  *
3109  * @see Expansion3D::AddFaceNormBoundaryInt
3110  *
3111  * @param Fn The trace quantities.
3112  * @param outarray Resulting 3D coefficient space.
3113  *
3114  */
3116  Array<OneD, NekDouble> &outarray)
3117 {
3118  if (m_expType == e1D)
3119  {
3120  int n, offset, t_offset;
3121 
3123  &elmtToTrace = m_traceMap->GetElmtToTrace();
3124 
3125  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
3126 
3127  for (n = 0; n < GetExpSize(); ++n)
3128  {
3129  // Number of coefficients on each element
3130  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3131 
3132  offset = GetCoeff_Offset(n);
3133 
3134  // Implementation for every points except Gauss points
3135  if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
3137  {
3138  t_offset =
3139  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3140  if (negatedFluxNormal[2 * n])
3141  {
3142  outarray[offset] -= Fn[t_offset];
3143  }
3144  else
3145  {
3146  outarray[offset] += Fn[t_offset];
3147  }
3148 
3149  t_offset =
3150  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3151 
3152  if (negatedFluxNormal[2 * n + 1])
3153  {
3154  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] -=
3155  Fn[t_offset];
3156  }
3157  else
3158  {
3159  outarray[offset + (*m_exp)[n]->GetVertexMap(1)] +=
3160  Fn[t_offset];
3161  }
3162  }
3163  else
3164  {
3165  int j;
3166  static DNekMatSharedPtr m_Ixm, m_Ixp;
3167  static int sav_ncoeffs = 0;
3168  if (!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
3169  {
3171  const LibUtilities::PointsKey BS_p(
3173  const LibUtilities::BasisKey BS_k(
3174  LibUtilities::eGauss_Lagrange, e_ncoeffs, BS_p);
3175 
3176  BASE = LibUtilities::BasisManager()[BS_k];
3177 
3178  Array<OneD, NekDouble> coords(1, 0.0);
3179 
3180  coords[0] = -1.0;
3181  m_Ixm = BASE->GetI(coords);
3182 
3183  coords[0] = 1.0;
3184  m_Ixp = BASE->GetI(coords);
3185 
3186  sav_ncoeffs = e_ncoeffs;
3187  }
3188 
3189  t_offset =
3190  GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
3191 
3192  if (negatedFluxNormal[2 * n])
3193  {
3194  for (j = 0; j < e_ncoeffs; j++)
3195  {
3196  outarray[offset + j] -=
3197  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3198  }
3199  }
3200  else
3201  {
3202  for (j = 0; j < e_ncoeffs; j++)
3203  {
3204  outarray[offset + j] +=
3205  (m_Ixm->GetPtr())[j] * Fn[t_offset];
3206  }
3207  }
3208 
3209  t_offset =
3210  GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
3211 
3212  if (negatedFluxNormal[2 * n + 1])
3213  {
3214  for (j = 0; j < e_ncoeffs; j++)
3215  {
3216  outarray[offset + j] -=
3217  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3218  }
3219  }
3220  else
3221  {
3222  for (j = 0; j < e_ncoeffs; j++)
3223  {
3224  outarray[offset + j] +=
3225  (m_Ixp->GetPtr())[j] * Fn[t_offset];
3226  }
3227  }
3228  }
3229  }
3230  }
3231  else // other dimensions
3232  {
3233  // Basis definition on each element
3234  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
3235  if ((m_expType != e1D) &&
3236  (basis->GetBasisType() != LibUtilities::eGauss_Lagrange))
3237  {
3238  Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
3239  m_trace->IProductWRTBase(Fn, Fcoeffs);
3240 
3241  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
3242  outarray);
3243  }
3244  else
3245  {
3246  int e, n, offset, t_offset;
3247  Array<OneD, NekDouble> e_outarray;
3249  &elmtToTrace = m_traceMap->GetElmtToTrace();
3250 
3251  if (m_expType == e2D)
3252  {
3253  for (n = 0; n < GetExpSize(); ++n)
3254  {
3255  offset = GetCoeff_Offset(n);
3256  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3257  {
3258  t_offset = GetTrace()->GetPhys_Offset(
3259  elmtToTrace[n][e]->GetElmtId());
3260  (*m_exp)[n]->AddEdgeNormBoundaryInt(
3261  e, elmtToTrace[n][e], Fn + t_offset,
3262  e_outarray = outarray + offset);
3263  }
3264  }
3265  }
3266  else if (m_expType == e3D)
3267  {
3268  for (n = 0; n < GetExpSize(); ++n)
3269  {
3270  offset = GetCoeff_Offset(n);
3271  for (e = 0; e < (*m_exp)[n]->GetNtraces(); ++e)
3272  {
3273  t_offset = GetTrace()->GetPhys_Offset(
3274  elmtToTrace[n][e]->GetElmtId());
3275  (*m_exp)[n]->AddFaceNormBoundaryInt(
3276  e, elmtToTrace[n][e], Fn + t_offset,
3277  e_outarray = outarray + offset);
3278  }
3279  }
3280  }
3281  }
3282  }
3283 }
3284 
3285 /**
3286  * @brief Add trace contributions into elemental coefficient spaces.
3287  *
3288  * Given some quantity \f$ \vec{q} \f$, calculate the elemental integral
3289  *
3290  * \f[
3291  * \int_{\Omega^e} \vec{q}, \mathrm{d}S
3292  * \f]
3293  *
3294  * and adds this to the coefficient space provided by
3295  * outarray. The value of q is determined from the routine
3296  * IsLeftAdjacentTrace() which if true we use Fwd else we use
3297  * Bwd
3298  *
3299  * @param Fwd The trace quantities associated with left (fwd)
3300  * adjancent elmt.
3301  * @param Bwd The trace quantities associated with right (bwd)
3302  * adjacent elet.
3303  * @param outarray Resulting coefficient space.
3304  */
3306  const Array<OneD, const NekDouble> &Fwd,
3308 {
3309 
3310  ASSERTL0(m_expType != e1D, "This method is not setup or "
3311  "tested for 1D expansion");
3312 
3313  Array<OneD, NekDouble> Coeffs(m_trace->GetNcoeffs());
3314 
3315  m_trace->IProductWRTBase(Fwd, Coeffs);
3316  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, Coeffs, outarray);
3317  m_trace->IProductWRTBase(Bwd, Coeffs);
3318  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, Coeffs, outarray);
3319 }
3320 
3322  Array<OneD, NekDouble> &outarray,
3323  const StdRegions::ConstFactorMap &factors,
3324  const StdRegions::VarCoeffMap &varcoeff,
3325  const MultiRegions::VarFactorsMap &varfactors,
3326  const Array<OneD, const NekDouble> &dirForcing,
3327  const bool PhysSpaceForcing)
3328 {
3329  boost::ignore_unused(varfactors, dirForcing);
3330  int i, n, cnt, nbndry;
3331  int nexp = GetExpSize();
3332 
3334  DNekVec F(m_ncoeffs, f, eWrapper);
3335  Array<OneD, NekDouble> e_f, e_l;
3336 
3337  //----------------------------------
3338  // Setup RHS Inner product if required
3339  //----------------------------------
3340  if (PhysSpaceForcing)
3341  {
3342  IProductWRTBase(inarray, f);
3343  Vmath::Neg(m_ncoeffs, f, 1);
3344  }
3345  else
3346  {
3347  Vmath::Smul(m_ncoeffs, -1.0, inarray, 1, f, 1);
3348  }
3349 
3350  //----------------------------------
3351  // Solve continuous Boundary System
3352  //----------------------------------
3353  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
3354  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
3355  int e_ncoeffs;
3356 
3358  NullAssemblyMapSharedPtr, factors, varcoeff);
3359  const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(HDGLamToUKey);
3360 
3361  // Retrieve number of local trace space coefficients N_{\lambda},
3362  // and set up local elemental trace solution \lambda^e.
3363  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
3364  Array<OneD, NekDouble> bndrhs(LocBndCoeffs, 0.0);
3365  Array<OneD, NekDouble> loclambda(LocBndCoeffs, 0.0);
3366  DNekVec LocLambda(LocBndCoeffs, loclambda, eWrapper);
3367 
3368  //----------------------------------
3369  // Evaluate Trace Forcing
3370  // Kirby et al, 2010, P23, Step 5.
3371  //----------------------------------
3372  // Determing <u_lam,f> terms using HDGLamToU matrix
3373  for (cnt = n = 0; n < nexp; ++n)
3374  {
3375  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
3376 
3377  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
3378  e_f = f + m_coeff_offset[n];
3379  e_l = bndrhs + cnt;
3380 
3381  // use outarray as tmp space
3382  DNekVec Floc(nbndry, e_l, eWrapper);
3383  DNekVec ElmtFce(e_ncoeffs, e_f, eWrapper);
3384  Floc = Transpose(*(HDGLamToU->GetBlock(n, n))) * ElmtFce;
3385 
3386  cnt += nbndry;
3387  }
3388 
3389  Array<OneD, const int> bndCondMap =
3390  m_traceMap->GetBndCondCoeffsToLocalTraceMap();
3391  Array<OneD, const NekDouble> Sign = m_traceMap->GetLocalToGlobalBndSign();
3392 
3393  // Copy Dirichlet boundary conditions and weak forcing
3394  // into trace space
3395  int locid;
3396  cnt = 0;
3397  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3398  {
3399  Array<OneD, const NekDouble> bndcoeffs =
3400  m_bndCondExpansions[i]->GetCoeffs();
3401 
3402  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3404  {
3405  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3406  {
3407  locid = bndCondMap[cnt + j];
3408  loclambda[locid] = Sign[locid] * bndcoeffs[j];
3409  }
3410  }
3411  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3413  m_bndConditions[i]->GetBoundaryConditionType() ==
3415  {
3416  // Add weak boundary condition to trace forcing
3417  for (int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
3418  {
3419  locid = bndCondMap[cnt + j];
3420  bndrhs[locid] += Sign[locid] * bndcoeffs[j];
3421  }
3422  }
3423 
3424  cnt += (m_bndCondExpansions[i])->GetNcoeffs();
3425  }
3426 
3427  //----------------------------------
3428  // Solve trace problem: \Lambda = K^{-1} F
3429  // K is the HybridDGHelmBndLam matrix.
3430  //----------------------------------
3431  if (GloBndDofs - NumDirBCs > 0)
3432  {
3434  factors, varcoeff);
3436  LinSys->Solve(bndrhs, loclambda, m_traceMap);
3437 
3438  // For consistency with previous version put global
3439  // solution into m_trace->m_coeffs
3440  m_traceMap->LocalToGlobal(loclambda, m_trace->UpdateCoeffs());
3441  }
3442 
3443  //----------------------------------
3444  // Internal element solves
3445  //----------------------------------
3447  NullAssemblyMapSharedPtr, factors, varcoeff);
3448 
3449  const DNekScalBlkMatSharedPtr &InvHDGHelm = GetBlockMatrix(invHDGhelmkey);
3450  DNekVec out(m_ncoeffs, outarray, eWrapper);
3451  Vmath::Zero(m_ncoeffs, outarray, 1);
3452 
3453  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
3454  out = (*InvHDGHelm) * F + (*HDGLamToU) * LocLambda;
3455 }
3456 
3457 /* \brief This function evaluates the boundary conditions at a certain
3458  * time-level.
3459  *
3460  * Based on the boundary condition \f$g(\boldsymbol{x},t)\f$ evaluated
3461  * at a given time-level \a t, this function transforms the boundary
3462  * conditions onto the coefficients of the (one-dimensional) boundary
3463  * expansion. Depending on the type of boundary conditions, these
3464  * expansion coefficients are calculated in different ways:
3465  * - <b>Dirichlet boundary conditions</b><BR>
3466  * In order to ensure global \f$C^0\f$ continuity of the spectral/hp
3467  * approximation, the Dirichlet boundary conditions are projected onto
3468  * the boundary expansion by means of a modified \f$C^0\f$ continuous
3469  * Galerkin projection. This projection can be viewed as a collocation
3470  * projection at the vertices, followed by an \f$L^2\f$ projection on
3471  * the interior modes of the edges. The resulting coefficients
3472  * \f$\boldsymbol{\hat{u}}^{\mathcal{D}}\f$ will be stored for the
3473  * boundary expansion.
3474  * - <b>Neumann boundary conditions</b>
3475  * In the discrete Galerkin formulation of the problem to be solved,
3476  * the Neumann boundary conditions appear as the set of surface
3477  * integrals: \f[\boldsymbol{\hat{g}}=\int_{\Gamma}
3478  * \phi^e_n(\boldsymbol{x})g(\boldsymbol{x})d(\boldsymbol{x})\quad
3479  * \forall n \f]
3480  * As a result, it are the coefficients \f$\boldsymbol{\hat{g}}\f$
3481  * that will be stored in the boundary expansion
3482  *
3483  * @param time The time at which the boundary conditions
3484  * should be evaluated.
3485  * @param bndCondExpansions List of boundary conditions.
3486  * @param bndConditions Information about the boundary conditions.
3487  *
3488  * This will only be undertaken for time dependent
3489  * boundary conditions unless time == 0.0 which is the
3490  * case when the method is called from the constructor.
3491  */
3493  const std::string varName,
3494  const NekDouble x2_in,
3495  const NekDouble x3_in)
3496 {
3497  int i;
3498  int npoints;
3499 
3500  MultiRegions::ExpListSharedPtr locExpList;
3501 
3502  for (i = 0; i < m_bndCondExpansions.size(); ++i)
3503  {
3504  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
3505  {
3506  m_bndCondBndWeight[i] = 1.0;
3507  locExpList = m_bndCondExpansions[i];
3508 
3509  npoints = locExpList->GetNpoints();
3510  Array<OneD, NekDouble> x0(npoints, 0.0);
3511  Array<OneD, NekDouble> x1(npoints, 0.0);
3512  Array<OneD, NekDouble> x2(npoints, 0.0);
3513 
3514  locExpList->GetCoords(x0, x1, x2);
3515 
3516  if (x2_in != NekConstants::kNekUnsetDouble &&
3518  {
3519  Vmath::Fill(npoints, x2_in, x1, 1);
3520  Vmath::Fill(npoints, x3_in, x2, 1);
3521  }
3522  else if (x2_in != NekConstants::kNekUnsetDouble)
3523  {
3524  Vmath::Fill(npoints, x2_in, x2, 1);
3525  }
3526 
3527  // treat 1D expansions separately since we only
3528  // require an evaluation at a point rather than
3529  // any projections or inner products that are not
3530  // available in a PointExp
3531  if (m_expType == e1D)
3532  {
3533  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3535  {
3536  m_bndCondExpansions[i]->SetCoeff(
3537  0, (std::static_pointer_cast<
3538  SpatialDomains ::DirichletBoundaryCondition>(
3539  m_bndConditions[i])
3540  ->m_dirichletCondition)
3541  .Evaluate(x0[0], x1[0], x2[0], time));
3542  m_bndCondExpansions[i]->SetPhys(
3543  0, m_bndCondExpansions[i]->GetCoeff(0));
3544  }
3545  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3547  {
3548  m_bndCondExpansions[i]->SetCoeff(
3549  0, (std::static_pointer_cast<
3550  SpatialDomains ::NeumannBoundaryCondition>(
3551  m_bndConditions[i])
3552  ->m_neumannCondition)
3553  .Evaluate(x0[0], x1[0], x2[0], time));
3554  }
3555  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3557  {
3558  m_bndCondExpansions[i]->SetCoeff(
3559  0, (std::static_pointer_cast<
3560  SpatialDomains ::RobinBoundaryCondition>(
3561  m_bndConditions[i])
3562  ->m_robinFunction)
3563  .Evaluate(x0[0], x1[0], x2[0], time));
3564  }
3565  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3567  {
3568  continue;
3569  }
3570  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3572  {
3573  }
3574  else
3575  {
3577  "This type of BC not implemented yet");
3578  }
3579  }
3580  else // 2D and 3D versions
3581  {
3582  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3584  {
3586  std::static_pointer_cast<
3588  m_bndConditions[i]);
3589 
3590  Array<OneD, NekDouble> valuesFile(npoints, 1.0),
3591  valuesExp(npoints, 1.0);
3592 
3593  string filebcs = bcPtr->m_filename;
3594  string exprbcs = bcPtr->m_expr;
3595 
3596  if (filebcs != "")
3597  {
3598  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3599  locExpList);
3600  valuesFile = locExpList->GetPhys();
3601  }
3602 
3603  if (exprbcs != "")
3604  {
3605  LibUtilities::Equation condition =
3606  std::static_pointer_cast<
3608  m_bndConditions[i])
3609  ->m_dirichletCondition;
3610 
3611  condition.Evaluate(x0, x1, x2, time, valuesExp);
3612  }
3613 
3614  Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
3615  locExpList->UpdatePhys(), 1);
3616 
3617  locExpList->FwdTransBndConstrained(
3618  locExpList->GetPhys(), locExpList->UpdateCoeffs());
3619  }
3620  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3622  {
3624  std::static_pointer_cast<
3626  m_bndConditions[i]);
3627  string filebcs = bcPtr->m_filename;
3628  if (filebcs != "")
3629  {
3630  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3631  locExpList);
3632  }
3633  else
3634  {
3635  LibUtilities::Equation condition =
3636  std::static_pointer_cast<
3638  m_bndConditions[i])
3639  ->m_neumannCondition;
3640  condition.Evaluate(x0, x1, x2, time,
3641  locExpList->UpdatePhys());
3642  }
3643 
3644  locExpList->IProductWRTBase(locExpList->GetPhys(),
3645  locExpList->UpdateCoeffs());
3646  }
3647  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3649  {
3651  std::static_pointer_cast<
3653  m_bndConditions[i]);
3654 
3655  string filebcs = bcPtr->m_filename;
3656 
3657  if (filebcs != "")
3658  {
3659  ExtractFileBCs(filebcs, bcPtr->GetComm(), varName,
3660  locExpList);
3661  }
3662  else
3663  {
3664  LibUtilities::Equation condition =
3665  std::static_pointer_cast<
3667  m_bndConditions[i])
3668  ->m_robinFunction;
3669  condition.Evaluate(x0, x1, x2, time,
3670  locExpList->UpdatePhys());
3671  }
3672 
3673  locExpList->IProductWRTBase(locExpList->GetPhys(),
3674  locExpList->UpdateCoeffs());
3675  }
3676  else if (m_bndConditions[i]->GetBoundaryConditionType() ==
3678  {
3679  continue;
3680  }
3681  else
3682  {
3684  "This type of BC not implemented yet");
3685  }
3686  }
3687  }
3688  }
3689 }
3690 
3691 /**
3692  * @brief Fill the weight with m_bndCondBndWeight.
3693  */
3695  Array<OneD, NekDouble> &weightjmp)
3696 {
3697  int cnt;
3698  int npts;
3699  int e = 0;
3700 
3701  // Fill boundary conditions into missing elements
3702  int id2 = 0;
3703 
3704  for (int n = cnt = 0; n < m_bndCondExpansions.size(); ++n)
3705  {
3706 
3707  if (m_bndConditions[n]->GetBoundaryConditionType() ==
3709  {
3710  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3711  {
3712  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3713  id2 = m_trace->GetPhys_Offset(
3714  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3715  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3716  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3717  }
3718 
3719  cnt += e;
3720  }
3721  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
3723  m_bndConditions[n]->GetBoundaryConditionType() ==
3725  {
3726  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
3727  {
3728  npts = m_bndCondExpansions[n]->GetExp(e)->GetTotPoints();
3729  id2 = m_trace->GetPhys_Offset(
3730  m_traceMap->GetBndCondIDToGlobalTraceID(cnt + e));
3731  Vmath::Fill(npts, m_bndCondBndWeight[n], &weightave[id2], 1);
3732  Vmath::Fill(npts, 0.0, &weightjmp[id2], 1);
3733  }
3734 
3735  cnt += e;
3736  }
3737  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
3739  {
3741  "Method not set up for this boundary condition.");
3742  }
3743  }
3744 }
3745 
3746 // Set up a list of element ids and trace ids that link to the
3747 // boundary conditions
3749  Array<OneD, int> &TraceID)
3750 {
3751 
3752  if (m_BCtoElmMap.size() == 0)
3753  {
3754  switch (m_expType)
3755  {
3756  case e1D:
3757  {
3758  map<int, int> VertGID;
3759  int i, n, id;
3760  int bid, cnt, Vid;
3761  int nbcs = m_bndConditions.size();
3762 
3763  // make sure arrays are of sufficient length
3764  m_BCtoElmMap = Array<OneD, int>(nbcs, -1);
3766 
3767  // setup map of all global ids along boundary
3768  cnt = 0;
3769  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3770  {
3771  Vid = m_bndCondExpansions[n]
3772  ->GetExp(0)
3773  ->GetGeom()
3774  ->GetVertex(0)
3775  ->GetVid();
3776  VertGID[Vid] = cnt++;
3777  }
3778 
3779  // Loop over elements and find verts that match;
3781  for (cnt = n = 0; n < GetExpSize(); ++n)
3782  {
3783  exp = (*m_exp)[n];
3784  for (i = 0; i < exp->GetNverts(); ++i)
3785  {
3786  id = exp->GetGeom()->GetVid(i);
3787 
3788  if (VertGID.count(id) > 0)
3789  {
3790  bid = VertGID.find(id)->second;
3791  ASSERTL1(m_BCtoElmMap[bid] == -1,
3792  "Edge already set");
3793  m_BCtoElmMap[bid] = n;
3794  m_BCtoTraceMap[bid] = i;
3795  cnt++;
3796  }
3797  }
3798  }
3799  ASSERTL1(cnt == nbcs,
3800  "Failed to visit all boundary condtiions");
3801  }
3802  break;
3803  case e2D:
3804  {
3805  map<int, int> globalIdMap;
3806  int i, n;
3807  int cnt;
3808  int nbcs = 0;
3809 
3810  // Populate global ID map (takes global geometry
3811  // ID to local expansion list ID).
3812  for (i = 0; i < GetExpSize(); ++i)
3813  {
3814  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3815  }
3816 
3817  // Determine number of boundary condition expansions.
3818  for (i = 0; i < m_bndConditions.size(); ++i)
3819  {
3820  nbcs += m_bndCondExpansions[i]->GetExpSize();
3821  }
3822 
3823  // Initialize arrays
3826 
3828  cnt = 0;
3829  for (n = 0; n < m_bndCondExpansions.size(); ++n)
3830  {
3831  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3832  ++i, ++cnt)
3833  {
3834  exp1d = m_bndCondExpansions[n]
3835  ->GetExp(i)
3836  ->as<LocalRegions::Expansion1D>();
3837 
3838  // Use edge to element map from MeshGraph.
3840  m_graph->GetElementsFromEdge(exp1d->GetGeom1D());
3841 
3842  m_BCtoElmMap[cnt] =
3843  globalIdMap[(*tmp)[0].first->GetGlobalID()];
3844  m_BCtoTraceMap[cnt] = (*tmp)[0].second;
3845  }
3846  }
3847  }
3848  break;
3849  case e3D:
3850  {
3851  map<int, int> globalIdMap;
3852  int i, n;
3853  int cnt;
3854  int nbcs = 0;
3855 
3856  // Populate global ID map (takes global geometry ID to local
3857  // expansion list ID).
3859  for (i = 0; i < GetExpSize(); ++i)
3860  {
3861  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3862  }
3863 
3864  // Determine number of boundary condition expansions.
3865  for (i = 0; i < m_bndConditions.size(); ++i)
3866  {
3867  nbcs += m_bndCondExpansions[i]->GetExpSize();
3868  }
3869 
3870  // Initialize arrays
3873 
3875  for (cnt = n = 0; n < m_bndCondExpansions.size(); ++n)
3876  {
3877  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
3878  ++i, ++cnt)
3879  {
3880  exp2d = m_bndCondExpansions[n]
3881  ->GetExp(i)
3882  ->as<LocalRegions::Expansion2D>();
3883 
3885  m_graph->GetElementsFromFace(exp2d->GetGeom2D());
3886  m_BCtoElmMap[cnt] =
3887  globalIdMap[tmp->at(0).first->GetGlobalID()];
3888  m_BCtoTraceMap[cnt] = tmp->at(0).second;
3889  }
3890  }
3891  }
3892  break;
3893  default:
3894  ASSERTL1(false, "Not setup for this expansion");
3895  break;
3896  }
3897  }
3898 
3899  ElmtID = m_BCtoElmMap;
3900  TraceID = m_BCtoTraceMap;
3901 }
3902 
3904  std::shared_ptr<ExpList> &result,
3905  const bool DeclareCoeffPhysArrays)
3906 {
3907  int n, cnt, nq;
3908  int offsetOld, offsetNew;
3909  std::vector<unsigned int> eIDs;
3910 
3911  Array<OneD, int> ElmtID, TraceID;
3912  GetBoundaryToElmtMap(ElmtID, TraceID);
3913 
3914  // Skip other boundary regions
3915  for (cnt = n = 0; n < i; ++n)
3916  {
3917  cnt += m_bndCondExpansions[n]->GetExpSize();
3918  }
3919 
3920  // Populate eIDs with information from BoundaryToElmtMap
3921  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
3922  {
3923  eIDs.push_back(ElmtID[cnt + n]);
3924  }
3925 
3926  // Create expansion list
3928  *this, eIDs, DeclareCoeffPhysArrays, Collections::eNoCollection);
3929 
3930  // Copy phys and coeffs to new explist
3931  if (DeclareCoeffPhysArrays)
3932  {
3933  Array<OneD, NekDouble> tmp1, tmp2;
3934  for (n = 0; n < result->GetExpSize(); ++n)
3935  {
3936  nq = GetExp(ElmtID[cnt + n])->GetTotPoints();
3937  offsetOld = GetPhys_Offset(ElmtID[cnt + n]);
3938  offsetNew = result->GetPhys_Offset(n);
3939  Vmath::Vcopy(nq, tmp1 = GetPhys() + offsetOld, 1,
3940  tmp2 = result->UpdatePhys() + offsetNew, 1);
3941 
3942  nq = GetExp(ElmtID[cnt + n])->GetNcoeffs();
3943  offsetOld = GetCoeff_Offset(ElmtID[cnt + n]);
3944  offsetNew = result->GetCoeff_Offset(n);
3945  Vmath::Vcopy(nq, tmp1 = GetCoeffs() + offsetOld, 1,
3946  tmp2 = result->UpdateCoeffs() + offsetNew, 1);
3947  }
3948  }
3949 }
3950 
3951 /**
3952  * @brief Reset this field, so that geometry information can be updated.
3953  */
3955 {
3956  ExpList::v_Reset();
3957 
3958  // Reset boundary condition expansions.
3959  for (int n = 0; n < m_bndCondExpansions.size(); ++n)
3960  {
3961  m_bndCondExpansions[n]->Reset();
3962  }
3963 }
3964 
3965 /**
3966  * Search through the edge expansions and identify which ones
3967  * have Robin/Mixed type boundary conditions. If find a Robin
3968  * boundary then store the edge id of the boundary condition
3969  * and the array of points of the physical space boundary
3970  * condition which are hold the boundary condition primitive
3971  * variable coefficient at the quatrature points
3972  *
3973  * \return std map containing the robin boundary condition
3974  * info using a key of the element id
3975  *
3976  * There is a next member to allow for more than one Robin
3977  * boundary condition per element
3978  */
3979 map<int, RobinBCInfoSharedPtr> DisContField::v_GetRobinBCInfo(void)
3980 {
3981  int i, cnt;
3982  map<int, RobinBCInfoSharedPtr> returnval;
3983  Array<OneD, int> ElmtID, TraceID;
3984  GetBoundaryToElmtMap(ElmtID, TraceID);
3985 
3986  for (cnt = i = 0; i < m_bndCondExpansions.size(); ++i)
3987  {
3988  MultiRegions::ExpListSharedPtr locExpList;
3989 
3990  if (m_bndConditions[i]->GetBoundaryConditionType() ==
3992  {
3993  int e, elmtid;
3994  Array<OneD, NekDouble> Array_tmp;
3995 
3996  locExpList = m_bndCondExpansions[i];
3997 
3998  int npoints = locExpList->GetNpoints();
3999  Array<OneD, NekDouble> x0(npoints, 0.0);
4000  Array<OneD, NekDouble> x1(npoints, 0.0);
4001  Array<OneD, NekDouble> x2(npoints, 0.0);
4002  Array<OneD, NekDouble> coeffphys(npoints);
4003 
4004  locExpList->GetCoords(x0, x1, x2);
4005 
4006  LibUtilities::Equation coeffeqn =
4007  std::static_pointer_cast<
4009  ->m_robinPrimitiveCoeff;
4010 
4011  // evalaute coefficient
4012  coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
4013 
4014  for (e = 0; e < locExpList->GetExpSize(); ++e)
4015  {
4016  RobinBCInfoSharedPtr rInfo =
4018  TraceID[cnt + e],
4019  Array_tmp = coeffphys + locExpList->GetPhys_Offset(e));
4020 
4021  elmtid = ElmtID[cnt + e];
4022  // make link list if necessary
4023  if (returnval.count(elmtid) != 0)
4024  {
4025  rInfo->next = returnval.find(elmtid)->second;
4026  }
4027  returnval[elmtid] = rInfo;
4028  }
4029  }
4030  cnt += m_bndCondExpansions[i]->GetExpSize();
4031  }
4032 
4033  return returnval;
4034 }
4035 
4036 /**
4037  * @brief Calculate the \f$ L^2 \f$ error of the \f$ Q_{\rm dir} \f$
4038  * derivative using the consistent DG evaluation of \f$ Q_{\rm dir} \f$.
4039  *
4040  * The solution provided is of the primative variation at the quadrature
4041  * points and the derivative is compared to the discrete derivative at
4042  * these points, which is likely to be undesirable unless using a much
4043  * higher number of quadrature points than the polynomial order used to
4044  * evaluate \f$ Q_{\rm dir} \f$.
4045  */
4047  const Array<OneD, const NekDouble> &soln)
4048 {
4049  int i, e, ncoeff_edge;
4050  Array<OneD, const NekDouble> tmp_coeffs;
4051  Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
4052 
4054  m_traceMap->GetElmtToTrace();
4055 
4056  StdRegions::Orientation edgedir;
4057 
4058  int cnt;
4059  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4060  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
4061 
4062  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4063 
4064  edge_lambda = loc_lambda;
4065 
4066  // Calculate Q using standard DG formulation.
4067  for (i = cnt = 0; i < GetExpSize(); ++i)
4068  {
4069  // Probably a better way of setting up lambda than this.
4070  // Note cannot use PutCoeffsInToElmts since lambda space
4071  // is mapped during the solve.
4072  int nEdges = (*m_exp)[i]->GetGeom()->GetNumEdges();
4073  Array<OneD, Array<OneD, NekDouble>> edgeCoeffs(nEdges);
4074 
4075  for (e = 0; e < nEdges; ++e)
4076  {
4077  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4078  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
4079  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
4080  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
4081  elmtToTrace[i][e]->SetCoeffsToOrientation(edgedir, edgeCoeffs[e],
4082  edgeCoeffs[e]);
4083  edge_lambda = edge_lambda + ncoeff_edge;
4084  }
4085 
4086  (*m_exp)[i]->DGDeriv(dir, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4087  elmtToTrace[i], edgeCoeffs, out_tmp = out_d + cnt);
4088  cnt += (*m_exp)[i]->GetNcoeffs();
4089  }
4090 
4091  BwdTrans(out_d, m_phys);
4092  Vmath::Vsub(m_npoints, m_phys, 1, soln, 1, m_phys, 1);
4093  return L2(m_phys);
4094 }
4095 
4096 /**
4097  * @brief Evaluate HDG post-processing to increase polynomial order of
4098  * solution.
4099  *
4100  * This function takes the solution (assumed to be one order lower) in
4101  * physical space, and postprocesses at the current polynomial order by
4102  * solving the system:
4103  *
4104  * \f[
4105  * \begin{aligned}
4106  * (\nabla w, \nabla u^*) &= (\nabla w, u), \\
4107  * \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
4108  * \end{aligned}
4109  * \f]
4110  *
4111  * where \f$ u \f$ corresponds with the current solution as stored
4112  * inside #m_coeffs.
4113  *
4114  * @param outarray The resulting field \f$ u^* \f$.
4115  */
4117 {
4118  int i, cnt, e, ncoeff_trace;
4119  Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
4121  m_traceMap->GetElmtToTrace();
4122 
4123  StdRegions::Orientation edgedir;
4124 
4125  int nq_elmt, nm_elmt;
4126  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
4127  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), trace_lambda;
4128  Array<OneD, NekDouble> tmp_coeffs;
4129  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(), loc_lambda);
4130 
4131  trace_lambda = loc_lambda;
4132 
4133  int dim = (m_expType == e2D) ? 2 : 3;
4134 
4135  int num_points[] = {0, 0, 0};
4136  int num_modes[] = {0, 0, 0};
4137 
4138  // Calculate Q using standard DG formulation.
4139  for (i = cnt = 0; i < GetExpSize(); ++i)
4140  {
4141  nq_elmt = (*m_exp)[i]->GetTotPoints();
4142  nm_elmt = (*m_exp)[i]->GetNcoeffs();
4143  qrhs = Array<OneD, NekDouble>(nq_elmt);
4144  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
4145  force = Array<OneD, NekDouble>(2 * nm_elmt);
4146  out_tmp = force + nm_elmt;
4148 
4149  for (int j = 0; j < dim; ++j)
4150  {
4151  num_points[j] = (*m_exp)[i]->GetBasis(j)->GetNumPoints();
4152  num_modes[j] = (*m_exp)[i]->GetBasis(j)->GetNumModes();
4153  }
4154 
4155  // Probably a better way of setting up lambda than this. Note
4156  // cannot use PutCoeffsInToElmts since lambda space is mapped
4157  // during the solve.
4158  int nTraces = (*m_exp)[i]->GetNtraces();
4159  Array<OneD, Array<OneD, NekDouble>> traceCoeffs(nTraces);
4160 
4161  for (e = 0; e < (*m_exp)[i]->GetNtraces(); ++e)
4162  {
4163  edgedir = (*m_exp)[i]->GetTraceOrient(e);
4164  ncoeff_trace = elmtToTrace[i][e]->GetNcoeffs();
4165  traceCoeffs[e] = Array<OneD, NekDouble>(ncoeff_trace);
4166  Vmath::Vcopy(ncoeff_trace, trace_lambda, 1, traceCoeffs[e], 1);
4167  if (dim == 2)
4168  {
4169  elmtToTrace[i][e]->SetCoeffsToOrientation(
4170  edgedir, traceCoeffs[e], traceCoeffs[e]);
4171  }
4172  else
4173  {
4174  (*m_exp)[i]
4176  ->SetFaceToGeomOrientation(e, traceCoeffs[e]);
4177  }
4178  trace_lambda = trace_lambda + ncoeff_trace;
4179  }
4180 
4181  // creating orthogonal expansion (checking if we have quads or
4182  // triangles)
4183  LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
4184  switch (shape)
4185  {
4187  {
4188  const LibUtilities::PointsKey PkeyQ1(
4189  num_points[0], LibUtilities::eGaussLobattoLegendre);
4190  const LibUtilities::PointsKey PkeyQ2(
4191  num_points[1], LibUtilities::eGaussLobattoLegendre);
4193  num_modes[0], PkeyQ1);
4195  num_modes[1], PkeyQ2);
4197  std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
4198  (*m_exp)[i]->GetGeom());
4200  BkeyQ1, BkeyQ2, qGeom);
4201  }
4202  break;
4204  {
4205  const LibUtilities::PointsKey PkeyT1(
4206  num_points[0], LibUtilities::eGaussLobattoLegendre);
4207  const LibUtilities::PointsKey PkeyT2(
4208  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4210  num_modes[0], PkeyT1);
4212  num_modes[1], PkeyT2);
4214  std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
4215  (*m_exp)[i]->GetGeom());
4217  BkeyT1, BkeyT2, tGeom);
4218  }
4219  break;
4221  {
4222  const LibUtilities::PointsKey PkeyH1(
4223  num_points[0], LibUtilities::eGaussLobattoLegendre);
4224  const LibUtilities::PointsKey PkeyH2(
4225  num_points[1], LibUtilities::eGaussLobattoLegendre);
4226  const LibUtilities::PointsKey PkeyH3(
4227  num_points[2], LibUtilities::eGaussLobattoLegendre);
4229  num_modes[0], PkeyH1);
4231  num_modes[1], PkeyH2);
4233  num_modes[2], PkeyH3);
4235  std::dynamic_pointer_cast<SpatialDomains::HexGeom>(
4236  (*m_exp)[i]->GetGeom());
4238  BkeyH1, BkeyH2, BkeyH3, hGeom);
4239  }
4240  break;
4242  {
4243  const LibUtilities::PointsKey PkeyT1(
4244  num_points[0], LibUtilities::eGaussLobattoLegendre);
4245  const LibUtilities::PointsKey PkeyT2(
4246  num_points[1], LibUtilities::eGaussRadauMAlpha1Beta0);
4247  const LibUtilities::PointsKey PkeyT3(
4248  num_points[2], LibUtilities::eGaussRadauMAlpha2Beta0);
4250  num_modes[0], PkeyT1);
4252  num_modes[1], PkeyT2);
4254  num_modes[2], PkeyT3);
4256  std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
4257  (*m_exp)[i]->GetGeom());
4259  BkeyT1, BkeyT2, BkeyT3, tGeom);
4260  }
4261  break;
4262  default:
4264  "Wrong shape type, HDG postprocessing is not "
4265  "implemented");
4266  };
4267 
4268  // DGDeriv
4269  // (d/dx w, d/dx q_0)
4270  (*m_exp)[i]->DGDeriv(0, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4271  elmtToTrace[i], traceCoeffs, out_tmp);
4272  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4273  ppExp->IProductWRTDerivBase(0, qrhs, force);
4274 
4275  // + (d/dy w, d/dy q_1)
4276  (*m_exp)[i]->DGDeriv(1, tmp_coeffs = m_coeffs + m_coeff_offset[i],
4277  elmtToTrace[i], traceCoeffs, out_tmp);
4278 
4279  (*m_exp)[i]->BwdTrans(out_tmp, qrhs);
4280  ppExp->IProductWRTDerivBase(1, qrhs, out_tmp);
4281 
4282  Vmath::Vadd(nm_elmt, force, 1, out_tmp, 1, force, 1);
4283 
4284  // determine force[0] = (1,u)
4285  (*m_exp)[i]->BwdTrans(tmp_coeffs = m_coeffs + m_coeff_offset[i], qrhs);
4286  force[0] = (*m_exp)[i]->Integral(qrhs);
4287 
4288  // multiply by inverse Laplacian matrix
4289  // get matrix inverse
4291  ppExp->DetShapeType(), *ppExp);
4292  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
4293 
4294  NekVector<NekDouble> in(nm_elmt, force, eWrapper);
4295  NekVector<NekDouble> out(nm_elmt);
4296 
4297  out = (*lapsys) * in;
4298 
4299  // Transforming back to modified basis
4300  Array<OneD, NekDouble> work(nq_elmt);
4301  ppExp->BwdTrans(out.GetPtr(), work);
4302  (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
4303  }
4304 }
4305 
4307  const Array<OneD, const NekDouble> &Fwd,
4308  const Array<OneD, const NekDouble> &Bwd,
4309  Array<OneD, NekDouble> &locTraceFwd, Array<OneD, NekDouble> &locTraceBwd)
4310 {
4311  if (locTraceFwd != NullNekDouble1DArray)
4312  {
4313  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(0, Fwd,
4314  locTraceFwd);
4315  }
4316 
4317  if (locTraceBwd != NullNekDouble1DArray)
4318  {
4319  m_locTraceToTraceMap->InterpLocTracesToTraceTranspose(1, Bwd,
4320  locTraceBwd);
4321  }
4322 }
4323 
4325  const Array<OneD, const NekDouble> &FwdFlux,
4326  const Array<OneD, const NekDouble> &BwdFlux,
4327  Array<OneD, NekDouble> &outarray)
4328 {
4329  Array<OneD, NekDouble> FCoeffs(m_trace->GetNcoeffs());
4330 
4331  m_trace->IProductWRTBase(FwdFlux, FCoeffs);
4332  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(1, FCoeffs, outarray);
4333  m_trace->IProductWRTBase(BwdFlux, FCoeffs);
4334  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(0, FCoeffs, outarray);
4335 }
4336 } // namespace MultiRegions
4337 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Describes the specification for a Basis.
Definition: Basis.h:50
NekDouble Evaluate() const
Definition: Equation.cpp:95
Interpreter class for the evaluation of mathematical expressions.
Definition: Interpreter.h:78
int DefineFunction(const std::string &vlist, const std::string &expr)
Defines a function for the purposes of evaluation.
NekDouble Evaluate(const int id)
Evaluate a function which depends only on constants and/or parameters.
Defines a specification for a set of points.
Definition: Points.h:59
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
This class is the abstractio n of a global discontinuous two- dimensional spectral/hp element expansi...
Definition: DisContField.h:56
std::vector< int > m_periodicBwdCopy
Definition: DisContField.h:186
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
Definition: DisContField.h:173
PeriodicMap m_periodicFaces
A map which identifies pairs of periodic faces.
Definition: DisContField.h:178
std::set< int > m_boundaryTraces
A set storing the global IDs of any boundary Verts.
Definition: DisContField.h:163
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Definition: DisContField.h:198
virtual ~DisContField()
Destructor.
std::vector< bool > m_negatedFluxNormal
Definition: DisContField.h:360
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Definition: DisContField.h:158
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
Array< OneD, int > m_BCtoTraceMap
Definition: DisContField.h:59
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:400
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
Add trace contributions into elemental coefficient spaces.
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: DisContField.h:384
void FindPeriodicTraces(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
Array< OneD, NekDouble > m_bndCondBndWeight
Definition: DisContField.h:145
bool SameTypeOfBoundaryConditions(const DisContField &In)
Check to see if expansion has the same BCs as In.
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
Definition: DisContField.h:149
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.
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill the weight with m_bndCondBndWeight.
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
Definition: DisContField.h:152
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
Definition: DisContField.h:168
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
ExpListSharedPtr m_trace
Trace space storage for points between elements.
Definition: DisContField.h:155
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)
Evaluate all boundary conditions at a given time..
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &TraceID)
std::vector< bool > & GetNegatedFluxNormal(void)
DisContField()
Default constructor.
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
Definition: DisContField.h:143
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
Definition: DisContField.h:185
bool IsLeftAdjacentTrace(const int n, const int e)
This routine determines if an element is to the "left" of the adjacent trace, which arises from the i...
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
std::vector< bool > m_leftAdjacentTraces
Definition: DisContField.h:192
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Solve the Helmholtz equation.
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, const bool DeclareCoeffPhysArrays=true)
Discretises the boundary conditions.
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2670
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1158
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1787
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2319
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1641
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:2010
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2450
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:5038
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:2100
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:2195
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:2232
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2204
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1210
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:3304
std::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:2431
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2173
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1184
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1126
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1196
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1136
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2437
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2223
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1132
NekDouble L2(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error with respect to soln of the global spectral/hp element approximat...
Definition: ExpList.h:551
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2658
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
Definition: ExpList.h:1849
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:2240
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1119
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1175
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:217
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:244
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:234
void Rotate(PointGeom &a, int dir, NekDouble angle)
_this = rotation of a by angle 'angle' around axis dir
Definition: PointGeom.cpp:142
NekDouble dist(PointGeom &a)
return distance between this and input a
Definition: PointGeom.cpp:180
static StdRegions::Orientation GetFaceOrientation(const QuadGeom &face1, const QuadGeom &face2, bool doRot=false, int dir=0, NekDouble angle=0.0, NekDouble tol=1e-8)
Get the orientation of face1.
Definition: QuadGeom.cpp:137
static StdRegions::Orientation GetFaceOrientation(const TriGeom &face1, const TriGeom &face2, bool doRot, int dir, NekDouble angle, NekDouble tol)
Definition: TriGeom.cpp:117
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:59
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:48
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:48
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:50
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:53
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1633
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
static const NekDouble kNekUnsetDouble
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:210
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshGraph.h:109
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
std::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:208
std::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:289
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:212
std::shared_ptr< std::vector< std::pair< GeometrySharedPtr, int > > > GeometryLinkSharedPtr
Definition: MeshGraph.h:168
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshGraph.h:110
std::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:213
std::vector< PointGeomSharedPtr > PointGeomVector
Definition: Geometry2D.h:64
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:138
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:219
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:84
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:172
std::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:214
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:218
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:86
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
std::shared_ptr< Composite > CompositeSharedPtr
Definition: MeshGraph.h:137
std::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:215
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:240
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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:209
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:518
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:895
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:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291
int id
Geometry ID of entity.
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
bool isLocal
Flag specifying if this entity is local to this partition.
NekDouble m_tol
Tolerance to rotation is considered identical.
int m_dir
Axis of rotation. 0 = 'x', 1 = 'y', 2 = 'z'.
NekDouble m_angle
Angle of rotation in radians.