Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DisContField2D.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File DisContField2D.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 2D domain with boundary conditions using
33 // LDG flux.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <LocalRegions/MatrixKey.h>
40 #include <LocalRegions/Expansion.h>
41 #include <LocalRegions/QuadExp.h>
42 #include <LocalRegions/TriExp.h>
46 
47 using namespace std;
48 
49 namespace Nektar
50 {
51  namespace MultiRegions
52  {
53  /**
54  * @class DisContField2D
55  * Abstraction of a global discontinuous two-dimensional spectral/hp
56  * element expansion which approximates the solution of a set of
57  * partial differential equations.
58  */
59 
60  /**
61  * @brief Default constructor.
62  */
63  DisContField2D::DisContField2D(void)
64  : ExpList2D (),
65  m_bndCondExpansions(),
66  m_bndConditions (),
67  m_trace (NullExpListSharedPtr)
68  {
69  }
70 
72  const DisContField2D &In,
73  const bool DeclareCoeffPhysArrays)
74  : ExpList2D (In,DeclareCoeffPhysArrays),
75  m_bndCondExpansions (In.m_bndCondExpansions),
76  m_bndConditions (In.m_bndConditions),
77  m_globalBndMat (In.m_globalBndMat),
78  m_traceMap (In.m_traceMap),
79  m_locTraceToTraceMap (In.m_locTraceToTraceMap),
80  m_boundaryEdges (In.m_boundaryEdges),
81  m_periodicVerts (In.m_periodicVerts),
82  m_periodicEdges (In.m_periodicEdges),
83  m_periodicFwdCopy (In.m_periodicFwdCopy),
84  m_periodicBwdCopy (In.m_periodicBwdCopy),
85  m_leftAdjacentEdges (In.m_leftAdjacentEdges)
86  {
87  if (In.m_trace)
88  {
90  *boost::dynamic_pointer_cast<ExpList1D>(In.m_trace),
91  DeclareCoeffPhysArrays);
92  }
93  }
94 
95  /**
96  * @brief Constructs a global discontinuous field based on an input
97  * mesh with boundary conditions.
98  */
100  const LibUtilities::SessionReaderSharedPtr &pSession,
101  const SpatialDomains::MeshGraphSharedPtr &graph2D,
102  const std::string &variable,
103  const bool SetUpJustDG,
104  const bool DeclareCoeffPhysArrays,
105  const Collections::ImplementationType ImpType)
106  : ExpList2D(pSession, graph2D, DeclareCoeffPhysArrays, variable,
107  ImpType),
108  m_bndCondExpansions(),
109  m_bndConditions(),
110  m_trace(NullExpListSharedPtr),
111  m_periodicVerts(),
112  m_periodicEdges(),
113  m_periodicFwdCopy(),
114  m_periodicBwdCopy()
115  {
116 
117  if (variable.compare("DefaultVar") != 0) // do not set up BCs if default variable
118  {
120  GenerateBoundaryConditionExpansion(graph2D, bcs, variable,
121  DeclareCoeffPhysArrays);
122 
123  if (DeclareCoeffPhysArrays)
124  {
125  EvaluateBoundaryConditions(0.0, variable);
126  }
127 
128  // Find periodic edges for this variable.
129  FindPeriodicEdges(bcs, variable);
130  }
131 
132  if (SetUpJustDG)
133  {
134  SetUpDG(variable);
135  }
136  else
137  {
138  // Set element edges to point to Robin BC edges if required.
139  int i, cnt;
140  Array<OneD, int> ElmtID, EdgeID;
141  GetBoundaryToElmtMap(ElmtID, EdgeID);
142 
143  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
144  {
146  int e;
147  locExpList = m_bndCondExpansions[i];
148 
149  for(e = 0; e < locExpList->GetExpSize(); ++e)
150  {
152  (*m_exp)[ElmtID[cnt+e]]->
153  as<LocalRegions::Expansion2D>();
155  locExpList->GetExp(e)->
156  as<LocalRegions::Expansion1D>();
158  locExpList->GetExp(e)->
159  as<LocalRegions::Expansion> ();
160 
161  exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
162  exp1d->SetAdjacentElementExp(EdgeID[cnt+e], exp2d);
163  }
164  cnt += m_bndCondExpansions[i]->GetExpSize();
165  }
166 
167  if(m_session->DefinesSolverInfo("PROJECTION"))
168  {
169  std::string ProjectStr =
170  m_session->GetSolverInfo("PROJECTION");
171  if((ProjectStr == "MixedCGDG") ||
172  (ProjectStr == "Mixed_CG_Discontinuous"))
173  {
174  SetUpDG();
175  }
176  else
177  {
179  }
180  }
181  else
182  {
184  }
185  }
186  }
187 
188 
189  /*
190  * @brief Copy type constructor which declares new boundary conditions
191  * and re-uses mapping info and trace space if possible
192  */
194  const DisContField2D &In,
195  const SpatialDomains::MeshGraphSharedPtr &graph2D,
196  const std::string &variable,
197  const bool SetUpJustDG,
198  const bool DeclareCoeffPhysArrays):
199  ExpList2D(In,DeclareCoeffPhysArrays),
200  m_trace(NullExpListSharedPtr)
201  {
202  // Set up boundary conditions for this variable.
203  // Do not set up BCs if default variable
204  if(variable.compare("DefaultVar") != 0)
205  {
207  GenerateBoundaryConditionExpansion(graph2D, bcs, variable);
208 
209  if (DeclareCoeffPhysArrays)
210  {
211  EvaluateBoundaryConditions(0.0, variable);
212  }
213 
215  {
216  // Find periodic edges for this variable.
217  FindPeriodicEdges(bcs, variable);
218 
219  if(SetUpJustDG)
220  {
221  SetUpDG();
222  }
223  else
224  {
225  // set elmt edges to point to robin bc edges if required
226  int i, cnt = 0;
227  Array<OneD, int> ElmtID,EdgeID;
228  GetBoundaryToElmtMap(ElmtID,EdgeID);
229 
230  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
231  {
233 
234  int e;
235  locExpList = m_bndCondExpansions[i];
236 
237  for(e = 0; e < locExpList->GetExpSize(); ++e)
238  {
240  = (*m_exp)[ElmtID[cnt+e]]->
241  as<LocalRegions::Expansion2D>();
243  = locExpList->GetExp(e)->
244  as<LocalRegions::Expansion1D>();
246  = locExpList->GetExp(e)->
247  as<LocalRegions::Expansion> ();
248 
249  exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
250  exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
251  exp2d);
252  }
253 
254  cnt += m_bndCondExpansions[i]->GetExpSize();
255  }
256 
257 
258  if (m_session->DefinesSolverInfo("PROJECTION"))
259  {
260  std::string ProjectStr =
261  m_session->GetSolverInfo("PROJECTION");
262 
263  if ((ProjectStr == "MixedCGDG") ||
264  (ProjectStr == "Mixed_CG_Discontinuous"))
265  {
266  SetUpDG();
267  }
268  else
269  {
271  }
272  }
273  else
274  {
276  }
277  }
278  }
279  else
280  {
281  if (SetUpJustDG)
282  {
284  m_trace = In.m_trace;
285  m_traceMap = In.m_traceMap;
293  }
294  else
295  {
297  m_trace = In.m_trace;
298  m_traceMap = In.m_traceMap;
306 
307  // set elmt edges to point to robin bc edges if required
308  int i, cnt = 0;
309  Array<OneD, int> ElmtID, EdgeID;
310  GetBoundaryToElmtMap(ElmtID, EdgeID);
311 
312  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
313  {
315 
316  int e;
317  locExpList = m_bndCondExpansions[i];
318 
319  for (e = 0; e < locExpList->GetExpSize(); ++e)
320  {
322  = (*m_exp)[ElmtID[cnt+e]]->
323  as<LocalRegions::Expansion2D>();
325  = locExpList->GetExp(e)->
326  as<LocalRegions::Expansion1D>();
328  = locExpList->GetExp(e)->
329  as<LocalRegions::Expansion> ();
330 
331  exp2d->SetEdgeExp(EdgeID[cnt+e], exp1d);
332  exp1d->SetAdjacentElementExp(EdgeID[cnt+e],
333  exp2d);
334  }
335  cnt += m_bndCondExpansions[i]->GetExpSize();
336  }
337 
339  }
340  }
341  }
342  }
343 
344  /**
345  * @brief Default destructor.
346  */
348  {
349  }
350 
352  const GlobalLinSysKey &mkey)
353  {
355  "Routine currently only tested for HybridDGHelmholtz");
356  ASSERTL1(mkey.GetGlobalSysSolnType() ==
357  m_traceMap->GetGlobalSysSolnType(),
358  "The local to global map is not set up for the requested "
359  "solution type");
360 
361  GlobalLinSysSharedPtr glo_matrix;
362  GlobalLinSysMap::iterator matrixIter = m_globalBndMat->find(mkey);
363 
364  if(matrixIter == m_globalBndMat->end())
365  {
366  glo_matrix = GenGlobalBndLinSys(mkey,m_traceMap);
367  (*m_globalBndMat)[mkey] = glo_matrix;
368  }
369  else
370  {
371  glo_matrix = matrixIter->second;
372  }
373 
374  return glo_matrix;
375  }
376 
377  /**
378  * @brief Set up all DG member variables and maps.
379  */
380  void DisContField2D::SetUpDG(const std::string variable)
381  {
382  // Check for multiple calls
384  {
385  return;
386  }
387 
388  ExpList1DSharedPtr trace;
390  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
391  m_graph);
392 
393  // Set up matrix map
396 
397  // Set up trace space
400  graph2D, m_periodicEdges);
401 
402  m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
404  AllocateSharedPtr(m_session, graph2D, trace, *this,
407  variable);
408 
409  if (m_session->DefinesCmdLineArgument("verbose"))
410  {
411  m_traceMap->PrintStats(std::cout, variable);
412  }
413 
415  &elmtToTrace = m_traceMap->GetElmtToTrace();
416 
417  // Scatter trace segments to 2D elements. For each element,
418  // we find the trace segment associated to each edge. The
419  // element then retains a pointer to the trace space segments,
420  // to ensure uniqueness of normals when retrieving from two
421  //adjoining elements which do not lie in a plane.
422  for (int i = 0; i < m_exp->size(); ++i)
423  {
424  for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
425  {
427  (*m_exp)[i]->as<LocalRegions::Expansion2D>();
429  elmtToTrace[i][j]->as<LocalRegions::Expansion1D>();
430  LocalRegions::ExpansionSharedPtr exp = elmtToTrace[i][j];;
431  exp2d->SetEdgeExp (j, exp1d);
432  exp1d->SetAdjacentElementExp(j, exp2d);
433  }
434  }
435 
436  // Set up physical normals
438 
439  // Set up information for parallel and periodic problems.
440  for (int i = 0; i < m_trace->GetExpSize(); ++i)
441  {
443  m_trace->GetExp(i)->as<LocalRegions::Expansion1D>();
444 
445  int offset = m_trace->GetPhys_Offset(i);
446  int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
447  PeriodicMap::iterator pIt = m_periodicEdges.find(traceGeomId);
448 
449  if (pIt != m_periodicEdges.end() && !pIt->second[0].isLocal)
450  {
451  if (traceGeomId != min(pIt->second[0].id, traceGeomId))
452  {
453  traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
454  traceEl->GetLeftAdjacentElementEdge());
455  }
456  }
457  else if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
458  {
459  traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
460  traceEl->GetLeftAdjacentElementEdge());
461  }
462  }
463 
464  int cnt, n, e;
465 
466  // Identify boundary edges
467  for (cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
468  {
469  if (m_bndConditions[n]->GetBoundaryConditionType() !=
471  {
472  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
473  {
474  m_boundaryEdges.insert(
475  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
476  }
477  }
478  cnt += m_bndCondExpansions[n]->GetExpSize();
479  }
480 
481  // Set up information for periodic boundary conditions.
482  boost::unordered_map<int,pair<int,int> > perEdgeToExpMap;
483  boost::unordered_map<int,pair<int,int> >::iterator it2;
484  for (cnt = n = 0; n < m_exp->size(); ++n)
485  {
486  for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
487  {
489  (*m_exp)[n]->GetGeom()->GetEid(e));
490 
491  if (it != m_periodicEdges.end())
492  {
493  perEdgeToExpMap[it->first] = make_pair(n, e);
494  }
495  }
496  }
497 
498  // Set up left-adjacent edge list.
499  m_leftAdjacentEdges.resize(cnt);
500  cnt = 0;
501  for (int i = 0; i < m_exp->size(); ++i)
502  {
503  for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j, ++cnt)
504  {
506  }
507  }
508 
509  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
510  cnt = 0;
511  for (int n = 0; n < m_exp->size(); ++n)
512  {
513  for (int e = 0; e < (*m_exp)[n]->GetNedges(); ++e, ++cnt)
514  {
515  int edgeGeomId = (*m_exp)[n]->GetGeom()->GetEid(e);
516  int offset = m_trace->GetPhys_Offset(
517  elmtToTrace[n][e]->GetElmtId());
518 
519  // Check to see if this face is periodic.
520  PeriodicMap::iterator it = m_periodicEdges.find(edgeGeomId);
521 
522  if (it != m_periodicEdges.end())
523  {
524  const PeriodicEntity &ent = it->second[0];
525  it2 = perEdgeToExpMap.find(ent.id);
526 
527  if (it2 == perEdgeToExpMap.end())
528  {
529  if (m_session->GetComm()->
530  GetRowComm()->GetSize() > 1 && !ent.isLocal)
531  {
532  continue;
533  }
534  else
535  {
536  ASSERTL1(false, "Periodic edge not found!");
537  }
538  }
539 
541  "Periodic edge in non-forward space?");
542 
543  int offset2 = m_trace->GetPhys_Offset(
544  elmtToTrace[it2->second.first][it2->second.second]->
545  GetElmtId());
546 
547  // Calculate relative orientations between edges to
548  // calculate copying map.
549  int nquad = elmtToTrace[n][e]->GetNumPoints(0);
550 
551  vector<int> tmpBwd(nquad);
552  vector<int> tmpFwd(nquad);
553 
554  if (ent.orient == StdRegions::eForwards)
555  {
556  for (int i = 0; i < nquad; ++i)
557  {
558  tmpBwd[i] = offset2 + i;
559  tmpFwd[i] = offset + i;
560  }
561  }
562  else
563  {
564  for (int i = 0; i < nquad; ++i)
565  {
566  tmpBwd[i] = offset2 + i;
567  tmpFwd[i] = offset + nquad - i - 1;
568  }
569  }
570 
571  for (int i = 0; i < nquad; ++i)
572  {
573  m_periodicFwdCopy.push_back(tmpFwd[i]);
574  m_periodicBwdCopy.push_back(tmpBwd[i]);
575  }
576  }
577  }
578  }
579 
581  AllocateSharedPtr(*this, m_trace, elmtToTrace,
583  }
584 
585  /**
586  * For each boundary region, checks that the types and number of
587  * boundary expansions in that region match.
588  * @param In ContField2D to compare with.
589  * @return True if boundary conditions match.
590  */
592  const DisContField2D &In)
593  {
594  int i;
595  bool returnval = true;
596 
597  for(i = 0; i < m_bndConditions.num_elements(); ++i)
598  {
599 
600  // check to see if boundary condition type is the same
601  // and there are the same number of boundary
602  // conditions in the boundary definition.
603  if((m_bndConditions[i]->GetBoundaryConditionType()
604  != In.m_bndConditions[i]->GetBoundaryConditionType())||
606  != In.m_bndCondExpansions[i]->GetExpSize()))
607  {
608  returnval = false;
609  break;
610  }
611  }
612 
613  // Compare with all other processes. Return true only if all
614  // processes report having the same boundary conditions.
615  int vSame = (returnval?1:0);
616  m_comm->GetRowComm()->AllReduce(vSame, LibUtilities::ReduceMin);
617 
618  return (vSame == 1);
619  }
620 
621  /**
622  * \brief This function discretises the boundary conditions by setting
623  * up a list of one-dimensional boundary expansions.
624  *
625  * According to their boundary region, the separate segmental boundary
626  * expansions are bundled together in an object of the class
627  * MultiRegions#ExpList1D.
628  *
629  * \param graph2D A mesh, containing information about the domain and
630  * the spectral/hp element expansion.
631  * \param bcs An entity containing information about the
632  * boundaries and boundary conditions.
633  * \param variable An optional parameter to indicate for which variable
634  * the boundary conditions should be discretised.
635  */
637  const SpatialDomains::MeshGraphSharedPtr &graph2D,
639  const std::string &variable,
640  const bool DeclareCoeffPhysArrays)
641  {
642  int cnt = 0;
646  bcs.GetBoundaryRegions();
647  const SpatialDomains::BoundaryConditionCollection &bconditions =
648  bcs.GetBoundaryConditions();
649  SpatialDomains::BoundaryRegionCollection::const_iterator it;
650 
651  // count the number of non-periodic boundary regions
652  for (it = bregions.begin(); it != bregions.end(); ++it)
653  {
654  bc = GetBoundaryCondition(bconditions, it->first, variable);
655 
656  if (bc->GetBoundaryConditionType() != SpatialDomains::ePeriodic)
657  {
658  cnt++;
659  }
660  }
661 
664  m_bndConditions =
666 
667  cnt = 0;
668 
669  // list non-periodic boundaries
670  for (it = bregions.begin(); it != bregions.end(); ++it)
671  {
672  bc = GetBoundaryCondition(bconditions, it->first, variable);
673 
674  if (bc->GetBoundaryConditionType() != SpatialDomains::ePeriodic)
675  {
677  ::AllocateSharedPtr(m_session, *(it->second), graph2D,
678  DeclareCoeffPhysArrays, variable);
679 
680  m_bndCondExpansions[cnt] = locExpList;
681  m_bndConditions[cnt] = bc;
682 
683 
684  std::string type = m_bndConditions[cnt]->GetUserDefined();
685 
686  // Set up normals on non-Dirichlet boundary
687  // conditions. Second two conditions ideally
688  // should be in local solver setup (when made into factory)
689  if((bc->GetBoundaryConditionType() !=
691  boost::iequals(type,"I") ||
692  boost::iequals(type,"CalcBC"))
693  {
695  }
696 
697  cnt++;
698  }
699  }
700  }
701 
702  /**
703  * @brief Determine the periodic edges and vertices for the given graph.
704  *
705  * Note that much of this routine is the same as the three-dimensional
706  * version, which therefore has much better documentation.
707  *
708  * @param bcs Information about the boundary conditions.
709  * @param variable Specifies the field.
710  *
711  * @see DisContField3D::FindPeriodicFaces
712  */
715  const std::string &variable)
716  {
718  = bcs.GetBoundaryRegions();
720  = bcs.GetBoundaryConditions();
722  = boost::dynamic_pointer_cast<
724  SpatialDomains::BoundaryRegionCollection::const_iterator it;
725 
727  m_session->GetComm()->GetRowComm();
729  m_session->GetCompositeOrdering();
730  LibUtilities::BndRegionOrdering bndRegOrder =
731  m_session->GetBndRegionOrdering();
733  m_graph->GetComposites();
734 
735  // Unique collection of pairs of periodic composites (i.e. if
736  // composites 1 and 2 are periodic then this map will contain either
737  // the pair (1,2) or (2,1) but not both).
738  map<int,int> perComps;
739  map<int,vector<int> > allVerts;
740  set<int> locVerts;
741  map<int,StdRegions::Orientation> allEdges;
742 
743  int region1ID, region2ID, i, j, k, cnt;
745 
746  // Set up a set of all local verts and edges.
747  for(i = 0; i < (*m_exp).size(); ++i)
748  {
749  for(j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
750  {
751  int id = (*m_exp)[i]->GetGeom()->GetVid(j);
752  locVerts.insert(id);
753  }
754  }
755 
756  // Construct list of all periodic pairs local to this process.
757  for (it = bregions.begin(); it != bregions.end(); ++it)
758  {
759  locBCond = GetBoundaryCondition(
760  bconditions, it->first, variable);
761 
762  if (locBCond->GetBoundaryConditionType()
764  {
765  continue;
766  }
767 
768  // Identify periodic boundary region IDs.
769  region1ID = it->first;
770  region2ID = boost::static_pointer_cast<
772  locBCond)->m_connectedBoundaryRegion;
773 
774  // From this identify composites. Note that in serial this will
775  // be an empty map.
776  int cId1, cId2;
777  if (vComm->GetSize() == 1)
778  {
779  cId1 = it->second->begin()->first;
780  cId2 = bregions.find(region2ID)->second->begin()->first;
781  }
782  else
783  {
784  cId1 = bndRegOrder.find(region1ID)->second[0];
785  cId2 = bndRegOrder.find(region2ID)->second[0];
786  }
787 
788  ASSERTL0(it->second->size() == 1,
789  "Boundary region "+boost::lexical_cast<string>(
790  region1ID)+" should only contain 1 composite.");
791 
792  // Construct set containing all periodic edges on this process
793  SpatialDomains::Composite c = it->second->begin()->second;
794 
795  vector<unsigned int> tmpOrder;
796 
797  for (i = 0; i < c->size(); ++i)
798  {
800  boost::dynamic_pointer_cast<
801  SpatialDomains::SegGeom>((*c)[i]);
802  ASSERTL0(segGeom, "Unable to cast to shared ptr");
803 
805  graph2D->GetElementsFromEdge(segGeom);
806  ASSERTL0(elmt->size() == 1,
807  "The periodic boundaries belong to "
808  "more than one element of the mesh");
809 
811  boost::dynamic_pointer_cast<SpatialDomains::Geometry2D>(
812  (*elmt)[0]->m_Element);
813 
814  allEdges[(*c)[i]->GetGlobalID()] =
815  geom->GetEorient((*elmt)[0]->m_EdgeIndx);
816 
817  // In serial mesh partitioning will not have occurred so
818  // need to fill composite ordering map manually.
819  if (vComm->GetSize() == 1)
820  {
821  tmpOrder.push_back((*c)[i]->GetGlobalID());
822  }
823 
824  vector<int> vertList(2);
825  vertList[0] = segGeom->GetVid(0);
826  vertList[1] = segGeom->GetVid(1);
827  allVerts[(*c)[i]->GetGlobalID()] = vertList;
828  }
829 
830  if (vComm->GetSize() == 1)
831  {
832  compOrder[it->second->begin()->first] = tmpOrder;
833  }
834 
835  // See if we already have either region1 or region2 stored in
836  // perComps map.
837  if (perComps.count(cId1) == 0)
838  {
839  if (perComps.count(cId2) == 0)
840  {
841  perComps[cId1] = cId2;
842  }
843  else
844  {
845  std::stringstream ss;
846  ss << "Boundary region " << cId2 << " should be "
847  << "periodic with " << perComps[cId2] << " but "
848  << "found " << cId1 << " instead!";
849  ASSERTL0(perComps[cId2] == cId1, ss.str());
850  }
851  }
852  else
853  {
854  std::stringstream ss;
855  ss << "Boundary region " << cId1 << " should be "
856  << "periodic with " << perComps[cId1] << " but "
857  << "found " << cId2 << " instead!";
858  ASSERTL0(perComps[cId1] == cId1, ss.str());
859  }
860  }
861 
862  // Process local edge list to obtain relative edge orientations.
863  int n = vComm->GetSize();
864  int p = vComm->GetRank();
865  int totEdges;
866  Array<OneD, int> edgecounts(n,0);
867  Array<OneD, int> edgeoffset(n,0);
868  Array<OneD, int> vertoffset(n,0);
869 
870  edgecounts[p] = allEdges.size();
871  vComm->AllReduce(edgecounts, LibUtilities::ReduceSum);
872 
873  edgeoffset[0] = 0;
874  for (i = 1; i < n; ++i)
875  {
876  edgeoffset[i] = edgeoffset[i-1] + edgecounts[i-1];
877  }
878 
879  totEdges = Vmath::Vsum(n, edgecounts, 1);
880  Array<OneD, int> edgeIds (totEdges, 0);
881  Array<OneD, int> edgeOrient(totEdges, 0);
882  Array<OneD, int> edgeVerts (totEdges, 0);
883 
885 
886  for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
887  {
888  edgeIds [edgeoffset[p] + i ] = sIt->first;
889  edgeOrient[edgeoffset[p] + i ] = sIt->second;
890  edgeVerts [edgeoffset[p] + i++] = allVerts[sIt->first].size();
891  }
892 
893  vComm->AllReduce(edgeIds, LibUtilities::ReduceSum);
894  vComm->AllReduce(edgeOrient, LibUtilities::ReduceSum);
895  vComm->AllReduce(edgeVerts, LibUtilities::ReduceSum);
896 
897  // Calculate number of vertices on each processor.
898  Array<OneD, int> procVerts(n,0);
899  int nTotVerts;
900 
901  // Note if there are no periodic edges at all calling Vsum will
902  // cause a segfault.
903  if (totEdges > 0)
904  {
905  nTotVerts = Vmath::Vsum(totEdges, edgeVerts, 1);
906  }
907  else
908  {
909  nTotVerts = 0;
910  }
911 
912  for (i = 0; i < n; ++i)
913  {
914  if (edgecounts[i] > 0)
915  {
916  procVerts[i] = Vmath::Vsum(
917  edgecounts[i], edgeVerts + edgeoffset[i], 1);
918  }
919  else
920  {
921  procVerts[i] = 0;
922  }
923  }
924  vertoffset[0] = 0;
925 
926  for (i = 1; i < n; ++i)
927  {
928  vertoffset[i] = vertoffset[i-1] + procVerts[i-1];
929  }
930 
931  Array<OneD, int> vertIds(nTotVerts, 0);
932  for (i = 0, sIt = allEdges.begin(); sIt != allEdges.end(); ++sIt)
933  {
934  for (j = 0; j < allVerts[sIt->first].size(); ++j)
935  {
936  vertIds[vertoffset[p] + i++] = allVerts[sIt->first][j];
937  }
938  }
939 
940  vComm->AllReduce(vertIds, LibUtilities::ReduceSum);
941 
942  // For simplicity's sake create a map of edge id -> orientation.
943  map<int, StdRegions::Orientation> orientMap;
944  map<int, vector<int> > vertMap;
945 
946  for (cnt = i = 0; i < totEdges; ++i)
947  {
948  ASSERTL0(orientMap.count(edgeIds[i]) == 0,
949  "Already found edge in orientation map!");
950  orientMap[edgeIds[i]] = (StdRegions::Orientation)edgeOrient[i];
951 
952  vector<int> verts(edgeVerts[i]);
953 
954  for (j = 0; j < edgeVerts[i]; ++j)
955  {
956  verts[j] = vertIds[cnt++];
957  }
958  vertMap[edgeIds[i]] = verts;
959  }
960 
961  // Go through list of composites and figure out which edges are
962  // parallel from original ordering in session file. This includes
963  // composites which are not necessarily on this process.
964  map<int,int>::iterator cIt, pIt;
965  map<int,int>::const_iterator oIt;
966 
967  map<int,int> allCompPairs;
968 
969  // Store temporary map of periodic vertices which will hold all
970  // periodic vertices on the entire mesh so that doubly periodic
971  // vertices can be counted properly across partitions. Local
972  // vertices are copied into m_periodicVerts at the end of the
973  // function.
974  PeriodicMap periodicVerts;
975 
976  for (cIt = perComps.begin(); cIt != perComps.end(); ++cIt)
977  {
979  const int id1 = cIt->first;
980  const int id2 = cIt->second;
981  std::string id1s = boost::lexical_cast<string>(id1);
982  std::string id2s = boost::lexical_cast<string>(id2);
983 
984  if (compMap.count(id1) > 0)
985  {
986  c[0] = compMap[id1];
987  }
988 
989  if (compMap.count(id2) > 0)
990  {
991  c[1] = compMap[id2];
992  }
993 
994  ASSERTL0(c[0] || c[1],
995  "Both composites not found on this process!");
996 
997  // Loop over composite ordering to construct list of all
998  // periodic edges regardless of whether they are on this
999  // process.
1000  map<int,int> compPairs;
1001 
1002  ASSERTL0(compOrder.count(id1) > 0,
1003  "Unable to find composite "+id1s+" in order map.");
1004  ASSERTL0(compOrder.count(id2) > 0,
1005  "Unable to find composite "+id2s+" in order map.");
1006  ASSERTL0(compOrder[id1].size() == compOrder[id2].size(),
1007  "Periodic composites "+id1s+" and "+id2s+
1008  " should have the same number of elements.");
1009  ASSERTL0(compOrder[id1].size() > 0,
1010  "Periodic composites "+id1s+" and "+id2s+
1011  " are empty!");
1012 
1013  // TODO: Add more checks.
1014  for (i = 0; i < compOrder[id1].size(); ++i)
1015  {
1016  int eId1 = compOrder[id1][i];
1017  int eId2 = compOrder[id2][i];
1018 
1019  ASSERTL0(compPairs.count(eId1) == 0,
1020  "Already paired.");
1021 
1022  if (compPairs.count(eId2) != 0)
1023  {
1024  ASSERTL0(compPairs[eId2] == eId1, "Pairing incorrect");
1025  }
1026  compPairs[eId1] = eId2;
1027  }
1028 
1029  // Construct set of all edges that we have locally on this
1030  // processor.
1031  set<int> locEdges;
1032  set<int>::iterator sIt;
1033  for (i = 0; i < 2; ++i)
1034  {
1035  if (!c[i])
1036  {
1037  continue;
1038  }
1039 
1040  if (c[i]->size() > 0)
1041  {
1042  for (j = 0; j < c[i]->size(); ++j)
1043  {
1044  locEdges.insert(c[i]->at(j)->GetGlobalID());
1045  }
1046  }
1047  }
1048 
1049  // Loop over all edges in the geometry composite.
1050  for (pIt = compPairs.begin(); pIt != compPairs.end(); ++pIt)
1051  {
1052  int ids [2] = {pIt->first, pIt->second};
1053  bool local[2] = {locEdges.count(pIt->first) > 0,
1054  locEdges.count(pIt->second) > 0};
1055 
1056  ASSERTL0(orientMap.count(ids[0]) > 0 &&
1057  orientMap.count(ids[1]) > 0,
1058  "Unable to find edge in orientation map.");
1059 
1060  allCompPairs[pIt->first ] = pIt->second;
1061  allCompPairs[pIt->second] = pIt->first;
1062 
1063  for (i = 0; i < 2; ++i)
1064  {
1065  if (!local[i])
1066  {
1067  continue;
1068  }
1069 
1070  int other = (i+1) % 2;
1071 
1073  orientMap[ids[i]] == orientMap[ids[other]] ?
1076 
1077  PeriodicEntity ent(ids [other], o,
1078  local[other]);
1079  m_periodicEdges[ids[i]].push_back(ent);
1080  }
1081 
1082  for (i = 0; i < 2; ++i)
1083  {
1084  int other = (i+1) % 2;
1085 
1087  orientMap[ids[i]] == orientMap[ids[other]] ?
1090 
1091  // Determine periodic vertices.
1092  vector<int> perVerts1 = vertMap[ids[i]];
1093  vector<int> perVerts2 = vertMap[ids[other]];
1094 
1095  map<int, pair<int, bool> > tmpMap;
1096  map<int, pair<int, bool> >::iterator mIt;
1097  if (o == StdRegions::eForwards)
1098  {
1099  tmpMap[perVerts1[0]] = make_pair(
1100  perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1101  tmpMap[perVerts1[1]] = make_pair(
1102  perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1103  }
1104  else
1105  {
1106  tmpMap[perVerts1[0]] = make_pair(
1107  perVerts2[1], locVerts.count(perVerts2[1]) > 0);
1108  tmpMap[perVerts1[1]] = make_pair(
1109  perVerts2[0], locVerts.count(perVerts2[0]) > 0);
1110  }
1111 
1112  for (mIt = tmpMap.begin(); mIt != tmpMap.end(); ++mIt)
1113  {
1114  // See if this vertex has been recorded already.
1115  PeriodicEntity ent2(mIt->second.first,
1117  mIt->second.second);
1118  PeriodicMap::iterator perIt = periodicVerts.find(
1119  mIt->first);
1120 
1121  if (perIt == periodicVerts.end())
1122  {
1123  periodicVerts[mIt->first].push_back(ent2);
1124  perIt = periodicVerts.find(mIt->first);
1125  }
1126  else
1127  {
1128  bool doAdd = true;
1129  for (j = 0; j < perIt->second.size(); ++j)
1130  {
1131  if (perIt->second[j].id == mIt->second.first)
1132  {
1133  doAdd = false;
1134  break;
1135  }
1136  }
1137 
1138  if (doAdd)
1139  {
1140  perIt->second.push_back(ent2);
1141  }
1142  }
1143  }
1144  }
1145  }
1146  }
1147 
1148  Array<OneD, int> pairSizes(n, 0);
1149  pairSizes[p] = allCompPairs.size();
1150  vComm->AllReduce(pairSizes, LibUtilities::ReduceSum);
1151 
1152  int totPairSizes = Vmath::Vsum(n, pairSizes, 1);
1153 
1154  Array<OneD, int> pairOffsets(n, 0);
1155  pairOffsets[0] = 0;
1156 
1157  for (i = 1; i < n; ++i)
1158  {
1159  pairOffsets[i] = pairOffsets[i-1] + pairSizes[i-1];
1160  }
1161 
1162  Array<OneD, int> first (totPairSizes, 0);
1163  Array<OneD, int> second(totPairSizes, 0);
1164 
1165  cnt = pairOffsets[p];
1166 
1167  for (pIt = allCompPairs.begin(); pIt != allCompPairs.end(); ++pIt)
1168  {
1169  first [cnt ] = pIt->first;
1170  second[cnt++] = pIt->second;
1171  }
1172 
1173  vComm->AllReduce(first, LibUtilities::ReduceSum);
1174  vComm->AllReduce(second, LibUtilities::ReduceSum);
1175 
1176  allCompPairs.clear();
1177 
1178  for(cnt = 0; cnt < totPairSizes; ++cnt)
1179  {
1180  allCompPairs[first[cnt]] = second[cnt];
1181  }
1182 
1183  // Search for periodic vertices and edges which are not in
1184  // a periodic composite but lie in this process. First, loop
1185  // over all information we have from other processors.
1186  for (cnt = i = 0; i < totEdges; ++i)
1187  {
1188  int edgeId = edgeIds[i];
1189 
1190  ASSERTL0(allCompPairs.count(edgeId) > 0,
1191  "Unable to find matching periodic edge.");
1192 
1193  int perEdgeId = allCompPairs[edgeId];
1194 
1195  for (j = 0; j < edgeVerts[i]; ++j, ++cnt)
1196  {
1197  int vId = vertIds[cnt];
1198 
1199  PeriodicMap::iterator perId = periodicVerts.find(vId);
1200 
1201  if (perId == periodicVerts.end())
1202  {
1203  // This vertex is not included in the map. Figure
1204  // out which vertex it is supposed to be periodic
1205  // with. perEdgeId is the edge ID which is periodic with
1206  // edgeId. The logic is much the same as the loop above.
1207  int perVertexId =
1208  orientMap[edgeId] == orientMap[perEdgeId] ?
1209  vertMap[perEdgeId][(j+1)%2] : vertMap[perEdgeId][j];
1210 
1211  PeriodicEntity ent(perVertexId,
1213  locVerts.count(perVertexId) > 0);
1214 
1215  periodicVerts[vId].push_back(ent);
1216  }
1217  }
1218  }
1219 
1220  // Loop over all periodic vertices to complete connectivity
1221  // information.
1222  PeriodicMap::iterator perIt, perIt2;
1223  for (perIt = periodicVerts.begin();
1224  perIt != periodicVerts.end(); ++perIt)
1225  {
1226  // Loop over associated vertices.
1227  for (i = 0; i < perIt->second.size(); ++i)
1228  {
1229  perIt2 = periodicVerts.find(perIt->second[i].id);
1230  ASSERTL0(perIt2 != periodicVerts.end(),
1231  "Couldn't find periodic vertex.");
1232 
1233  for (j = 0; j < perIt2->second.size(); ++j)
1234  {
1235  if (perIt2->second[j].id == perIt->first)
1236  {
1237  continue;
1238  }
1239 
1240  bool doAdd = true;
1241  for (k = 0; k < perIt->second.size(); ++k)
1242  {
1243  if (perIt2->second[j].id == perIt->second[k].id)
1244  {
1245  doAdd = false;
1246  break;
1247  }
1248  }
1249 
1250  if (doAdd)
1251  {
1252  perIt->second.push_back(perIt2->second[j]);
1253  }
1254  }
1255  }
1256  }
1257 
1258  // Do one final loop over periodic vertices to remove non-local
1259  // vertices from map.
1260  for (perIt = periodicVerts.begin();
1261  perIt != periodicVerts.end(); ++perIt)
1262  {
1263  if (locVerts.count(perIt->first) > 0)
1264  {
1265  m_periodicVerts.insert(*perIt);
1266  }
1267  }
1268  }
1269 
1270  bool DisContField2D::IsLeftAdjacentEdge(const int n, const int e)
1271  {
1272  set<int>::iterator it;
1274  m_traceMap->GetElmtToTrace()[n][e]->
1275  as<LocalRegions::Expansion1D>();
1276 
1277 
1278  bool fwd = true;
1279  if (traceEl->GetLeftAdjacentElementEdge () == -1 ||
1280  traceEl->GetRightAdjacentElementEdge() == -1)
1281  {
1282  // Boundary edge (1 connected element). Do nothing in
1283  // serial.
1284  it = m_boundaryEdges.find(traceEl->GetElmtId());
1285 
1286  // If the edge does not have a boundary condition set on
1287  // it, then assume it is a partition edge.
1288  if (it == m_boundaryEdges.end())
1289  {
1290  int traceGeomId = traceEl->GetGeom1D()->GetGlobalID();
1292  traceGeomId);
1293 
1294  if (pIt != m_periodicEdges.end() && !pIt->second[0].isLocal)
1295  {
1296  fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
1297  }
1298  else
1299  {
1300  int offset = m_trace->GetPhys_Offset(traceEl->GetElmtId());
1301 
1302  fwd = m_traceMap->
1303  GetTraceToUniversalMapUnique(offset) >= 0;
1304  }
1305  }
1306  }
1307  else if (traceEl->GetLeftAdjacentElementEdge () != -1 &&
1308  traceEl->GetRightAdjacentElementEdge() != -1)
1309  {
1310  // Non-boundary edge (2 connected elements).
1311  fwd = dynamic_cast<Nektar::StdRegions::StdExpansion*>
1312  (traceEl->GetLeftAdjacentElementExp().get()) ==
1313  (*m_exp)[n].get();
1314  }
1315  else
1316  {
1317  ASSERTL2(false, "Unconnected trace element!");
1318  }
1319 
1320  return fwd;
1321  }
1322 
1323  // Construct the two trace vectors of the inner and outer
1324  // trace solution from the field contained in m_phys, where
1325  // the Weak dirichlet boundary conditions are listed in the
1326  // outer part of the vecotr
1330  {
1331  v_GetFwdBwdTracePhys(m_phys, Fwd, Bwd);
1332  }
1333 
1334  /**
1335  * @brief This method extracts the "forward" and "backward" trace data
1336  * from the array @a field and puts the data into output vectors @a Fwd
1337  * and @a Bwd.
1338  *
1339  * We first define the convention which defines "forwards" and
1340  * "backwards". First an association is made between the edge of each
1341  * element and its corresponding edge in the trace space using the
1342  * mapping #m_traceMap. The element can either be left-adjacent or
1343  * right-adjacent to this trace edge (see
1344  * Expansion1D::GetLeftAdjacentElementExp). Boundary edges are always
1345  * left-adjacent since left-adjacency is populated first.
1346  *
1347  * If the element is left-adjacent we extract the edge trace data from
1348  * @a field into the forward trace space @a Fwd; otherwise, we place it
1349  * in the backwards trace space @a Bwd. In this way, we form a unique
1350  * set of trace normals since these are always extracted from
1351  * left-adjacent elements.
1352  *
1353  * @param field is a NekDouble array which contains the 2D data
1354  * from which we wish to extract the backward and
1355  * forward orientated trace/edge arrays.
1356  * @param Fwd The resulting forwards space.
1357  * @param Bwd The resulting backwards space.
1358  */
1360  const Array<OneD, const NekDouble> &field,
1363  {
1364  int cnt, n, e, npts, phys_offset;
1365 
1366  // Zero forward/backward vectors.
1367  Vmath::Zero(Fwd.num_elements(), Fwd, 1);
1368  Vmath::Zero(Bwd.num_elements(), Bwd, 1);
1369 
1370  // Basis definition on each element
1371  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
1372  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
1373  {
1374 
1375  // blocked routine
1377  GetNLocTracePts());
1378 
1379  m_locTraceToTraceMap->LocTracesFromField(field, edgevals);
1380  m_locTraceToTraceMap->InterpLocEdgesToTrace(0, edgevals, Fwd);
1381 
1382  Array<OneD, NekDouble> invals = edgevals + m_locTraceToTraceMap->
1383  GetNFwdLocTracePts();
1384  m_locTraceToTraceMap->InterpLocEdgesToTrace(1, invals, Bwd);
1385  }
1386  else
1387  {
1388  // Loop over elements and collect forward expansion
1389  int nexp = GetExpSize();
1390  Array<OneD,NekDouble> e_tmp;
1392  boost::unordered_map<int,pair<int,int> >::iterator it3;
1394 
1396  &elmtToTrace = m_traceMap->GetElmtToTrace();
1397 
1398  for(cnt = n = 0; n < nexp; ++n)
1399  {
1400  exp2d = (*m_exp)[n]->as<LocalRegions::Expansion2D>();
1401  phys_offset = GetPhys_Offset(n);
1402 
1403  for(e = 0; e < exp2d->GetNedges(); ++e, ++cnt)
1404  {
1405  int offset = m_trace->GetPhys_Offset(
1406  elmtToTrace[n][e]->GetElmtId());
1407 
1408  if (m_leftAdjacentEdges[cnt])
1409  {
1410  exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1411  field + phys_offset,
1412  e_tmp = Fwd + offset);
1413  }
1414  else
1415  {
1416  exp2d->GetEdgePhysVals(e, elmtToTrace[n][e],
1417  field + phys_offset,
1418  e_tmp = Bwd + offset);
1419  }
1420 
1421  }
1422  }
1423  }
1424 
1425  // Fill boundary conditions into missing elements
1426  int id1, id2 = 0;
1427 
1428  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1429  {
1430  if (m_bndConditions[n]->GetBoundaryConditionType() ==
1432  {
1433  for (e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
1434  {
1435  npts = m_bndCondExpansions[n]->GetExp(e)->GetNumPoints(0);
1436  id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
1437  id2 = m_trace->GetPhys_Offset(m_traceMap->
1438  GetBndCondTraceToGlobalTraceMap(cnt+e));
1439  Vmath::Vcopy(npts,
1440  &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
1441  &Bwd[id2], 1);
1442  }
1443 
1444  cnt += e;
1445  }
1446  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
1448  m_bndConditions[n]->GetBoundaryConditionType() ==
1450  {
1451  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
1452  {
1453  npts = m_bndCondExpansions[n]->GetExp(e)->GetNumPoints(0);
1454  id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
1455  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[id1]==0.0,
1456  "Method not set up for non-zero Neumann "
1457  "boundary condition");
1458  id2 = m_trace->GetPhys_Offset(
1459  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
1460  Vmath::Vcopy(npts, &Fwd[id2], 1, &Bwd[id2], 1);
1461  }
1462 
1463  cnt += e;
1464  }
1465  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
1467  {
1468  ASSERTL0(false,
1469  "Method not set up for this boundary condition.");
1470  }
1471  }
1472 
1473  // Copy any periodic boundary conditions.
1474  for (n = 0; n < m_periodicFwdCopy.size(); ++n)
1475  {
1476  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
1477  }
1478 
1479  // Do parallel exchange for forwards/backwards spaces.
1480  m_traceMap->UniversalTraceAssemble(Fwd);
1481  m_traceMap->UniversalTraceAssemble(Bwd);
1482  }
1483 
1484 
1486  Array<OneD, NekDouble> &outarray)
1487  {
1488  ASSERTL1(m_physState == true,
1489  "Field must be in physical state to extract trace space.");
1490 
1491  v_ExtractTracePhys(m_phys, outarray);
1492  }
1493 
1494  /**
1495  * @brief This method extracts the trace (edges in 2D) from the field @a
1496  * inarray and puts the values in @a outarray.
1497  *
1498  * It assumes the field is C0 continuous so that it can overwrite the
1499  * edge data when visited by the two adjacent elements.
1500  *
1501  * @param inarray An array containing the 2D data from which we wish
1502  * to extract the edge data.
1503  * @param outarray The resulting edge information.
1504  */
1506  const Array<OneD, const NekDouble> &inarray,
1507  Array<OneD, NekDouble> &outarray)
1508  {
1509  // Loop over elemente and collect forward expansion
1510  int nexp = GetExpSize();
1511  int n, e, offset, phys_offset;
1512  Array<OneD,NekDouble> e_tmp;
1514  &elmtToTrace = m_traceMap->GetElmtToTrace();
1515 
1516  ASSERTL1(outarray.num_elements() >= m_trace->GetNpoints(),
1517  "input array is of insufficient length");
1518 
1519  // use m_trace tmp space in element to fill values
1520  for(n = 0; n < nexp; ++n)
1521  {
1522  phys_offset = GetPhys_Offset(n);
1523 
1524  for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1525  {
1526  offset = m_trace->GetPhys_Offset(
1527  elmtToTrace[n][e]->GetElmtId());
1528  (*m_exp)[n]->GetEdgePhysVals(e, elmtToTrace[n][e],
1529  inarray + phys_offset,
1530  e_tmp = outarray + offset);
1531  }
1532  }
1533  }
1534 
1536  const Array<OneD, const NekDouble> &Fx,
1537  const Array<OneD, const NekDouble> &Fy,
1538  Array<OneD, NekDouble> &outarray)
1539  {
1540  int e, n, offset, t_offset;
1541  Array<OneD, NekDouble> e_outarray;
1543  &elmtToTrace = m_traceMap->GetElmtToTrace();
1544 
1545  for(n = 0; n < GetExpSize(); ++n)
1546  {
1547  offset = GetCoeff_Offset(n);
1548  for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1549  {
1550  t_offset = GetTrace()->GetPhys_Offset(
1551  elmtToTrace[n][e]->GetElmtId());
1552 
1553  (*m_exp)[n]->AddEdgeNormBoundaryInt(
1554  e,elmtToTrace[n][e],
1555  Fx + t_offset,
1556  Fy + t_offset,
1557  e_outarray = outarray+offset);
1558  }
1559  }
1560  }
1561 
1562  /**
1563  * @brief Add trace contributions into elemental coefficient spaces.
1564  *
1565  * Given some quantity \f$ \vec{Fn} \f$, which conatins this
1566  * routine calculates the integral
1567  *
1568  * \f[
1569  * \int_{\Omega^e} \vec{Fn}, \mathrm{d}S
1570  * \f]
1571  *
1572  * and adds this to the coefficient space provided by outarray.
1573  *
1574  * @see Expansion2D::AddEdgeNormBoundaryInt
1575  *
1576  * @param Fn The trace quantities.
1577  * @param outarray Resulting 2D coefficient space.
1578  */
1580  const Array<OneD, const NekDouble> &Fn,
1581  Array<OneD, NekDouble> &outarray)
1582  {
1583  // Basis definition on each element
1584  LibUtilities::BasisSharedPtr basis = (*m_exp)[0]->GetBasis(0);
1585  if (basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
1586  {
1587  Array<OneD, NekDouble> Fcoeffs(m_trace->GetNcoeffs());
1588  m_trace->IProductWRTBase(Fn, Fcoeffs);
1589 
1590  m_locTraceToTraceMap->AddTraceCoeffsToFieldCoeffs(Fcoeffs,
1591  outarray);
1592  }
1593  else
1594  {
1595  int e, n, offset, t_offset;
1596  Array<OneD, NekDouble> e_outarray;
1598  &elmtToTrace = m_traceMap->GetElmtToTrace();
1599 
1600  for(n = 0; n < GetExpSize(); ++n)
1601  {
1602  offset = GetCoeff_Offset(n);
1603  for(e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1604  {
1605  t_offset = GetTrace()->GetPhys_Offset(
1606  elmtToTrace[n][e]->GetElmtId());
1607  (*m_exp)[n]->AddEdgeNormBoundaryInt(
1608  e, elmtToTrace[n][e],
1609  Fn+t_offset,
1610  e_outarray = outarray+offset);
1611  }
1612  }
1613  }
1614  }
1615 
1616 
1617  /**
1618  * @brief Add trace contributions into elemental coefficient spaces.
1619  *
1620  * Given some quantity \f$ \vec{q} \f$, calculate the elemental integral
1621  *
1622  * \f[
1623  * \int_{\Omega^e} \vec{q}, \mathrm{d}S
1624  * \f]
1625  *
1626  * and adds this to the coefficient space provided by
1627  * outarray. The value of q is determined from the routine
1628  * IsLeftAdjacentEdge() which if true we use Fwd else we use
1629  * Bwd
1630  *
1631  * @see Expansion2D::AddEdgeNormBoundaryInt
1632  *
1633  * @param Fwd The trace quantities associated with left (fwd)
1634  * adjancent elmt.
1635  * @param Bwd The trace quantities associated with right (bwd)
1636  * adjacent elet.
1637  * @param outarray Resulting 2D coefficient space.
1638  */
1640  const Array<OneD, const NekDouble> &Fwd,
1641  const Array<OneD, const NekDouble> &Bwd,
1642  Array<OneD, NekDouble> &outarray)
1643  {
1644  int e,n,offset, t_offset;
1645  Array<OneD, NekDouble> e_outarray;
1647  &elmtToTrace = m_traceMap->GetElmtToTrace();
1648 
1649  for (n = 0; n < GetExpSize(); ++n)
1650  {
1651  offset = GetCoeff_Offset(n);
1652  for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
1653  {
1654  t_offset = GetTrace()->GetPhys_Offset(
1655  elmtToTrace[n][e]->GetElmtId());
1656 
1657  // Evaluate upwind flux less local edge
1658  if (IsLeftAdjacentEdge(n, e))
1659  {
1660  (*m_exp)[n]->AddEdgeNormBoundaryInt(
1661  e, elmtToTrace[n][e], Fwd+t_offset,
1662  e_outarray = outarray+offset);
1663  }
1664  else
1665  {
1666  (*m_exp)[n]->AddEdgeNormBoundaryInt(
1667  e, elmtToTrace[n][e], Bwd+t_offset,
1668  e_outarray = outarray+offset);
1669  }
1670 
1671  }
1672  }
1673  }
1674 
1675  /**
1676  * @brief Set up a list of element IDs and edge IDs that link to the
1677  * boundary conditions.
1678  */
1680  Array<OneD, int> &ElmtID,
1681  Array<OneD, int> &EdgeID)
1682  {
1683  if (m_BCtoElmMap.num_elements() == 0)
1684  {
1685  map<int, int> globalIdMap;
1686  int i,n;
1687  int cnt;
1688  int nbcs = 0;
1689 
1691  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
1692  m_graph);
1693 
1694  // Populate global ID map (takes global geometry ID to local
1695  // expansion list ID).
1696  for (i = 0; i < GetExpSize(); ++i)
1697  {
1698  globalIdMap[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1699  }
1700 
1701  // Determine number of boundary condition expansions.
1702  for(i = 0; i < m_bndConditions.num_elements(); ++i)
1703  {
1704  nbcs += m_bndCondExpansions[i]->GetExpSize();
1705  }
1706 
1707  // Initialize arrays
1710 
1712  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1713  {
1714  for (i = 0; i < m_bndCondExpansions[n]->GetExpSize();
1715  ++i, ++cnt)
1716  {
1717  exp1d = m_bndCondExpansions[n]->GetExp(i)->
1718  as<LocalRegions::Expansion1D>();
1719  // Use edge to element map from MeshGraph2D.
1721  graph2D->GetElementsFromEdge(exp1d->GetGeom1D());
1722 
1723  m_BCtoElmMap[cnt] = globalIdMap[(*tmp)[0]->
1724  m_Element->GetGlobalID()];
1725  m_BCtoEdgMap[cnt] = (*tmp)[0]->m_EdgeIndx;
1726  }
1727  }
1728  }
1729  ElmtID = m_BCtoElmMap;
1730  EdgeID = m_BCtoEdgMap;
1731  }
1732 
1734  boost::shared_ptr<ExpList> &result,
1735  const bool DeclareCoeffPhysArrays)
1736  {
1737  int n, cnt, nq;
1738  int offsetOld, offsetNew;
1739  std::vector<unsigned int> eIDs;
1740 
1741  Array<OneD, int> ElmtID,EdgeID;
1742  GetBoundaryToElmtMap(ElmtID,EdgeID);
1743 
1744  // Skip other boundary regions
1745  for (cnt = n = 0; n < i; ++n)
1746  {
1747  cnt += m_bndCondExpansions[n]->GetExpSize();
1748  }
1749 
1750  // Populate eIDs with information from BoundaryToElmtMap
1751  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
1752  {
1753  eIDs.push_back(ElmtID[cnt+n]);
1754  }
1755 
1756  // Create expansion list
1757  result =
1759  (*this, eIDs, DeclareCoeffPhysArrays);
1760 
1761  // Copy phys and coeffs to new explist
1762  if( DeclareCoeffPhysArrays)
1763  {
1764  Array<OneD, NekDouble> tmp1, tmp2;
1765  for (n = 0; n < result->GetExpSize(); ++n)
1766  {
1767  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
1768  offsetOld = GetPhys_Offset(ElmtID[cnt+n]);
1769  offsetNew = result->GetPhys_Offset(n);
1770  Vmath::Vcopy(nq, tmp1 = GetPhys()+ offsetOld, 1,
1771  tmp2 = result->UpdatePhys()+ offsetNew, 1);
1772 
1773  nq = GetExp(ElmtID[cnt+n])->GetNcoeffs();
1774  offsetOld = GetCoeff_Offset(ElmtID[cnt+n]);
1775  offsetNew = result->GetCoeff_Offset(n);
1776  Vmath::Vcopy(nq, tmp1 = GetCoeffs()+ offsetOld, 1,
1777  tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
1778  }
1779  }
1780  }
1781 
1782  /**
1783  * @brief Reset this field, so that geometry information can be updated.
1784  */
1786  {
1787  ExpList::v_Reset();
1788 
1789  // Reset boundary condition expansions.
1790  for (int n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1791  {
1792  m_bndCondExpansions[n]->Reset();
1793  }
1794  }
1795 
1796  /**
1797  * @brief Calculate the \f$ L^2 \f$ error of the \f$ Q_{\rm dir} \f$
1798  * derivative using the consistent DG evaluation of \f$ Q_{\rm dir} \f$.
1799  *
1800  * The solution provided is of the primative variation at the quadrature
1801  * points and the derivative is compared to the discrete derivative at
1802  * these points, which is likely to be undesirable unless using a much
1803  * higher number of quadrature points than the polynomial order used to
1804  * evaluate \f$ Q_{\rm dir} \f$.
1805  */
1807  const int dir,
1808  const Array<OneD, const NekDouble> &soln)
1809  {
1810  int i,e,ncoeff_edge;
1811  Array<OneD, const NekDouble> tmp_coeffs;
1812  Array<OneD, NekDouble> out_d(m_ncoeffs), out_tmp;
1813 
1815  &elmtToTrace = m_traceMap->GetElmtToTrace();
1816 
1817  StdRegions::Orientation edgedir;
1818 
1819  int cnt;
1820  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
1821  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
1822  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
1823 
1824  edge_lambda = loc_lambda;
1825 
1826  // Calculate Q using standard DG formulation.
1827  for(i = cnt = 0; i < GetExpSize(); ++i)
1828  {
1829  // Probably a better way of setting up lambda than this.
1830  // Note cannot use PutCoeffsInToElmts since lambda space
1831  // is mapped during the solve.
1832  int nEdges = (*m_exp)[i]->GetNedges();
1833  Array<OneD, Array<OneD, NekDouble> > edgeCoeffs(nEdges);
1834 
1835  for(e = 0; e < nEdges; ++e)
1836  {
1837  edgedir = (*m_exp)[i]->GetEorient(e);
1838  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
1839  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
1840  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
1841  elmtToTrace[i][e]->SetCoeffsToOrientation(
1842  edgedir, edgeCoeffs[e], edgeCoeffs[e]);
1843  edge_lambda = edge_lambda + ncoeff_edge;
1844  }
1845 
1846  (*m_exp)[i]->DGDeriv(dir,
1847  tmp_coeffs=m_coeffs+m_coeff_offset[i],
1848  elmtToTrace[i],
1849  edgeCoeffs,
1850  out_tmp = out_d+cnt);
1851  cnt += (*m_exp)[i]->GetNcoeffs();
1852  }
1853 
1854  BwdTrans(out_d,m_phys);
1855  Vmath::Vsub(m_npoints,m_phys,1,soln,1,m_phys,1);
1856  return L2(m_phys);
1857  }
1858 
1860  const Array<OneD, const NekDouble> &inarray,
1861  Array<OneD, NekDouble> &outarray,
1862  const FlagList &flags,
1863  const StdRegions::ConstFactorMap &factors,
1864  const StdRegions::VarCoeffMap &varcoeff,
1865  const Array<OneD, const NekDouble> &dirForcing,
1866  const bool PhysSpaceForcing)
1867 
1868  {
1869  int i,j,n,cnt,cnt1,nbndry;
1870  int nexp = GetExpSize();
1872 
1874  DNekVec F(m_ncoeffs,f,eWrapper);
1875  Array<OneD,NekDouble> e_f, e_l;
1876 
1877  //----------------------------------
1878  // Setup RHS Inner product
1879  //----------------------------------
1880  if(PhysSpaceForcing)
1881  {
1882  IProductWRTBase(inarray,f);
1883  Vmath::Neg(m_ncoeffs,f,1);
1884  }
1885  else
1886  {
1887  Vmath::Smul(m_ncoeffs,-1.0,inarray,1,f,1);
1888  }
1889 
1890  //----------------------------------
1891  // Solve continuous flux System
1892  //----------------------------------
1893  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
1894  int NumDirichlet = m_traceMap->GetNumLocalDirBndCoeffs();
1895  int e_ncoeffs,id;
1896 
1897  // Retrieve block matrix of U^e
1899  NullAssemblyMapSharedPtr,factors,varcoeff);
1900  const DNekScalBlkMatSharedPtr &HDGLamToU = GetBlockMatrix(
1901  HDGLamToUKey);
1902 
1903  // Retrieve global trace space storage, \Lambda, from trace
1904  // expansion
1905  Array<OneD,NekDouble> BndSol = m_trace->UpdateCoeffs();
1906 
1907  // Create trace space forcing, F
1908  Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
1909 
1910  // Zero \Lambda
1911  Vmath::Zero(GloBndDofs,BndSol,1);
1912 
1913  // Retrieve number of local trace space coefficients N_{\lambda},
1914  // and set up local elemental trace solution \lambda^e.
1915  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
1916  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
1917  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
1918 
1919  //----------------------------------
1920  // Evaluate Trace Forcing vector F
1921  // Kirby et al, 2010, P23, Step 5.
1922  //----------------------------------
1923  // Loop over all expansions in the domain
1924  for(cnt = cnt1 = n = 0; n < nexp; ++n)
1925  {
1926  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1927 
1928  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1929  e_f = f + cnt;
1930  e_l = loc_lambda + cnt1;
1931 
1932  // Local trace space \lambda^e
1933  DNekVec Floc (nbndry, e_l, eWrapper);
1934  // Local forcing f^e
1935  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
1936  // Compute local (U^e)^{\top} f^e
1937  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1938 
1939  cnt += e_ncoeffs;
1940  cnt1 += nbndry;
1941  }
1942 
1943  // Assemble local \lambda_e into global \Lambda
1944  m_traceMap->AssembleBnd(loc_lambda,BndRhs);
1945 
1946  // Copy Dirichlet boundary conditions and weak forcing into trace
1947  // space
1948  cnt = 0;
1949  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1950  {
1951  if(m_bndConditions[i]->GetBoundaryConditionType() ==
1953  {
1954  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
1955  {
1956  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1957  BndSol[id] = m_bndCondExpansions[i]->GetCoeffs()[j];
1958  }
1959  }
1960  else
1961  {
1962  //Add weak boundary condition to trace forcing
1963  for(j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
1964  {
1965  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(cnt++);
1966  BndRhs[id] += m_bndCondExpansions[i]->GetCoeffs()[j];
1967  }
1968  }
1969  }
1970 
1971  //----------------------------------
1972  // Solve trace problem: \Lambda = K^{-1} F
1973  // K is the HybridDGHelmBndLam matrix.
1974  //----------------------------------
1975  if(GloBndDofs - NumDirichlet > 0)
1976  {
1978  m_traceMap,factors,varcoeff);
1980  LinSys->Solve(BndRhs,BndSol,m_traceMap);
1981  }
1982 
1983  //----------------------------------
1984  // Internal element solves
1985  //----------------------------------
1987  NullAssemblyMapSharedPtr,factors,varcoeff);
1988  const DNekScalBlkMatSharedPtr& InvHDGHelm = GetBlockMatrix(
1989  invHDGhelmkey);
1990  DNekVec out(m_ncoeffs,outarray,eWrapper);
1991  Vmath::Zero(m_ncoeffs,outarray,1);
1992 
1993  // get local trace solution from BndSol
1994  m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1995 
1996  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
1997  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
1998  }
1999 
2000 
2001  /**
2002  * @brief Calculates the result of the multiplication of a global matrix
2003  * of type specified by @a mkey with a vector given by @a inarray.
2004  *
2005  * @param mkey Key representing desired matrix multiplication.
2006  * @param inarray Input vector.
2007  * @param outarray Resulting multiplication.
2008  */
2010  const GlobalMatrixKey &gkey,
2011  const Array<OneD,const NekDouble> &inarray,
2012  Array<OneD, NekDouble> &outarray,
2013  CoeffState coeffstate)
2014  {
2015  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
2016  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
2017  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
2018  const DNekScalBlkMatSharedPtr& HDGHelm = GetBlockMatrix(gkey);
2019 
2020  m_traceMap->GlobalToLocalBnd(inarray, loc_lambda);
2021  LocLambda = (*HDGHelm) * LocLambda;
2022  m_traceMap->AssembleBnd(loc_lambda,outarray);
2023  }
2024 
2025  /**
2026  * @brief Search through the edge expansions and identify which ones
2027  * have Robin/Mixed type boundary conditions.
2028  *
2029  * If a Robin boundary is found then store the edge ID of the boundary
2030  * condition and the array of points of the physical space boundary
2031  * condition which are hold the boundary condition primitive variable
2032  * coefficient at the quatrature points
2033  *
2034  * @return A map containing the Robin boundary condition information
2035  * using a key of the element ID.
2036  */
2037  map<int, RobinBCInfoSharedPtr> DisContField2D::v_GetRobinBCInfo(void)
2038  {
2039  int i,cnt;
2040  map<int, RobinBCInfoSharedPtr> returnval;
2041  Array<OneD, int> ElmtID,EdgeID;
2042  GetBoundaryToElmtMap(ElmtID,EdgeID);
2043 
2044  for(cnt = i = 0; i < m_bndCondExpansions.num_elements(); ++i)
2045  {
2046  MultiRegions::ExpListSharedPtr locExpList;
2047 
2048  if(m_bndConditions[i]->GetBoundaryConditionType() ==
2050  {
2051  int e,elmtid;
2052  Array<OneD, NekDouble> Array_tmp;
2053 
2054  locExpList = m_bndCondExpansions[i];
2055 
2056  int npoints = locExpList->GetNpoints();
2057  Array<OneD, NekDouble> x0(npoints, 0.0);
2058  Array<OneD, NekDouble> x1(npoints, 0.0);
2059  Array<OneD, NekDouble> x2(npoints, 0.0);
2060  Array<OneD, NekDouble> coeffphys(npoints);
2061 
2062  locExpList->GetCoords(x0, x1, x2);
2063 
2064  LibUtilities::Equation coeffeqn =
2065  boost::static_pointer_cast<
2067  (m_bndConditions[i])->m_robinPrimitiveCoeff;
2068 
2069  // evalaute coefficient
2070  coeffeqn.Evaluate(x0, x1, x2, 0.0, coeffphys);
2071 
2072  for(e = 0; e < locExpList->GetExpSize(); ++e)
2073  {
2074  RobinBCInfoSharedPtr rInfo =
2077  EdgeID[cnt+e],
2078  Array_tmp = coeffphys +
2079  locExpList->GetPhys_Offset(e));
2080 
2081  elmtid = ElmtID[cnt+e];
2082  // make link list if necessary
2083  if(returnval.count(elmtid) != 0)
2084  {
2085  rInfo->next = returnval.find(elmtid)->second;
2086  }
2087  returnval[elmtid] = rInfo;
2088  }
2089  }
2090  cnt += m_bndCondExpansions[i]->GetExpSize();
2091  }
2092 
2093  return returnval;
2094  }
2095 
2096  /**
2097  * @brief Evaluate HDG post-processing to increase polynomial order of
2098  * solution.
2099  *
2100  * This function takes the solution (assumed to be one order lower) in
2101  * physical space, and postprocesses at the current polynomial order by
2102  * solving the system:
2103  *
2104  * \f[
2105  * \begin{aligned}
2106  * (\nabla w, \nabla u^*) &= (\nabla w, u), \\
2107  * \langle \nabla u^*, 1 \rangle &= \langle \nabla u, 1 \rangle
2108  * \end{aligned}
2109  * \f]
2110  *
2111  * where \f$ u \f$ corresponds with the current solution as stored
2112  * inside #m_coeffs.
2113  *
2114  * @param outarray The resulting field \f$ u^* \f$.
2115  */
2117  Array<OneD, NekDouble> &outarray)
2118  {
2119  int i,cnt,e,ncoeff_edge;
2120  Array<OneD, NekDouble> force, out_tmp, qrhs, qrhs1;
2122  &elmtToTrace = m_traceMap->GetElmtToTrace();
2123 
2124  StdRegions::Orientation edgedir;
2125 
2126  int nq_elmt, nm_elmt;
2127  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
2128  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs), edge_lambda;
2129  Array<OneD, NekDouble> tmp_coeffs;
2130  m_traceMap->GlobalToLocalBnd(m_trace->GetCoeffs(),loc_lambda);
2131 
2132  edge_lambda = loc_lambda;
2133 
2134  // Calculate Q using standard DG formulation.
2135  for(i = cnt = 0; i < GetExpSize(); ++i)
2136  {
2137  nq_elmt = (*m_exp)[i]->GetTotPoints();
2138  nm_elmt = (*m_exp)[i]->GetNcoeffs();
2139  qrhs = Array<OneD, NekDouble>(nq_elmt);
2140  qrhs1 = Array<OneD, NekDouble>(nq_elmt);
2141  force = Array<OneD, NekDouble>(2*nm_elmt);
2142  out_tmp = force + nm_elmt;
2144 
2145  int num_points0 = (*m_exp)[i]->GetBasis(0)->GetNumPoints();
2146  int num_points1 = (*m_exp)[i]->GetBasis(1)->GetNumPoints();
2147  int num_modes0 = (*m_exp)[i]->GetBasis(0)->GetNumModes();
2148  int num_modes1 = (*m_exp)[i]->GetBasis(1)->GetNumModes();
2149 
2150  // Probably a better way of setting up lambda than this. Note
2151  // cannot use PutCoeffsInToElmts since lambda space is mapped
2152  // during the solve.
2153  int nEdges = (*m_exp)[i]->GetNedges();
2154  Array<OneD, Array<OneD, NekDouble> > edgeCoeffs(nEdges);
2155 
2156  for(e = 0; e < (*m_exp)[i]->GetNedges(); ++e)
2157  {
2158  edgedir = (*m_exp)[i]->GetEorient(e);
2159  ncoeff_edge = elmtToTrace[i][e]->GetNcoeffs();
2160  edgeCoeffs[e] = Array<OneD, NekDouble>(ncoeff_edge);
2161  Vmath::Vcopy(ncoeff_edge, edge_lambda, 1, edgeCoeffs[e], 1);
2162  elmtToTrace[i][e]->SetCoeffsToOrientation(
2163  edgedir, edgeCoeffs[e], edgeCoeffs[e]);
2164  edge_lambda = edge_lambda + ncoeff_edge;
2165  }
2166 
2167  //creating orthogonal expansion (checking if we have quads or triangles)
2168  LibUtilities::ShapeType shape = (*m_exp)[i]->DetShapeType();
2169  switch(shape)
2170  {
2172  {
2175  LibUtilities::BasisKey BkeyQ1(LibUtilities::eOrtho_A, num_modes0, PkeyQ1);
2176  LibUtilities::BasisKey BkeyQ2(LibUtilities::eOrtho_A, num_modes1, PkeyQ2);
2177  SpatialDomains::QuadGeomSharedPtr qGeom = boost::dynamic_pointer_cast<SpatialDomains::QuadGeom>((*m_exp)[i]->GetGeom());
2178  ppExp = MemoryManager<LocalRegions::QuadExp>::AllocateSharedPtr(BkeyQ1, BkeyQ2, qGeom);
2179  }
2180  break;
2182  {
2185  LibUtilities::BasisKey BkeyT1(LibUtilities::eOrtho_A, num_modes0, PkeyT1);
2186  LibUtilities::BasisKey BkeyT2(LibUtilities::eOrtho_B, num_modes1, PkeyT2);
2187  SpatialDomains::TriGeomSharedPtr tGeom = boost::dynamic_pointer_cast<SpatialDomains::TriGeom>((*m_exp)[i]->GetGeom());
2188  ppExp = MemoryManager<LocalRegions::TriExp>::AllocateSharedPtr(BkeyT1, BkeyT2, tGeom);
2189  }
2190  break;
2191  default:
2192  ASSERTL0(false, "Wrong shape type, HDG postprocessing is not implemented");
2193  };
2194 
2195 
2196  //DGDeriv
2197  // (d/dx w, d/dx q_0)
2198  (*m_exp)[i]->DGDeriv(
2199  0,tmp_coeffs = m_coeffs + m_coeff_offset[i],
2200  elmtToTrace[i], edgeCoeffs, out_tmp);
2201  (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2202  //(*m_exp)[i]->IProductWRTDerivBase(0,qrhs,force);
2203  ppExp->IProductWRTDerivBase(0,qrhs,force);
2204 
2205 
2206  // + (d/dy w, d/dy q_1)
2207  (*m_exp)[i]->DGDeriv(
2208  1,tmp_coeffs = m_coeffs + m_coeff_offset[i],
2209  elmtToTrace[i], edgeCoeffs, out_tmp);
2210 
2211  (*m_exp)[i]->BwdTrans(out_tmp,qrhs);
2212  //(*m_exp)[i]->IProductWRTDerivBase(1,qrhs,out_tmp);
2213  ppExp->IProductWRTDerivBase(1,qrhs,out_tmp);
2214 
2215  Vmath::Vadd(nm_elmt,force,1,out_tmp,1,force,1);
2216 
2217  // determine force[0] = (1,u)
2218  (*m_exp)[i]->BwdTrans(
2219  tmp_coeffs = m_coeffs + m_coeff_offset[i],qrhs);
2220  force[0] = (*m_exp)[i]->Integral(qrhs);
2221 
2222  // multiply by inverse Laplacian matrix
2223  // get matrix inverse
2224  LocalRegions::MatrixKey lapkey(StdRegions::eInvLaplacianWithUnityMean, ppExp->DetShapeType(), *ppExp);
2225  DNekScalMatSharedPtr lapsys = ppExp->GetLocMatrix(lapkey);
2226 
2227  NekVector<NekDouble> in (nm_elmt,force,eWrapper);
2228  NekVector<NekDouble> out(nm_elmt);
2229 
2230  out = (*lapsys)*in;
2231 
2232  // Transforming back to modified basis
2233  Array<OneD, NekDouble> work(nq_elmt);
2234  ppExp->BwdTrans(out.GetPtr(), work);
2235  (*m_exp)[i]->FwdTrans(work, tmp_coeffs = outarray + m_coeff_offset[i]);
2236  }
2237  }
2238 
2239  /**
2240  * Evaluates the boundary condition expansions, \a bndCondExpansions,
2241  * given the information provided by \a bndConditions.
2242  * @param time The time at which the boundary conditions
2243  * should be evaluated.
2244  * @param bndCondExpansions List of boundary conditions.
2245  * @param bndConditions Information about the boundary conditions.
2246  *
2247  * This will only be undertaken for time dependent
2248  * boundary conditions unless time == 0.0 which is the
2249  * case when the method is called from the constructor.
2250  */
2252  const NekDouble time,
2253  const std::string varName,
2254  const NekDouble x2_in,
2255  const NekDouble x3_in)
2256  {
2257  int i;
2258  int npoints;
2259  int nbnd = m_bndCondExpansions.num_elements();
2260 
2261  MultiRegions::ExpListSharedPtr locExpList;
2262 
2263  for (i = 0; i < nbnd; ++i)
2264  {
2265  if (time == 0.0 ||
2266  m_bndConditions[i]->IsTimeDependent())
2267  {
2268  locExpList = m_bndCondExpansions[i];
2269  npoints = locExpList->GetNpoints();
2270  Array<OneD, NekDouble> x0(npoints, 0.0);
2271  Array<OneD, NekDouble> x1(npoints, 0.0);
2272  Array<OneD, NekDouble> x2(npoints, 0.0);
2273 
2274  // Homogeneous input case for x2.
2275  if (x2_in == NekConstants::kNekUnsetDouble)
2276  {
2277  locExpList->GetCoords(x0, x1, x2);
2278  }
2279  else
2280  {
2281  locExpList->GetCoords(x0, x1, x2);
2282  Vmath::Fill(npoints, x2_in, x2, 1);
2283  }
2284 
2285  if (m_bndConditions[i]->GetBoundaryConditionType()
2287  {
2289  = boost::static_pointer_cast<
2291  m_bndConditions[i]);
2292  string filebcs = bcPtr->m_filename;
2293 
2294  if (filebcs != "")
2295  {
2296  ExtractFileBCs(filebcs, bcPtr->m_comm, varName, locExpList);
2297  }
2298  else
2299  {
2300  LibUtilities::Equation condition =
2301  boost::static_pointer_cast<
2303  (m_bndConditions[i])->
2304  m_dirichletCondition;
2305 
2306  condition.Evaluate(x0, x1, x2, time,
2307  locExpList->UpdatePhys());
2308  }
2309 
2310  locExpList->FwdTrans_BndConstrained(
2311  locExpList->GetPhys(),
2312  locExpList->UpdateCoeffs());
2313  }
2314  else if (m_bndConditions[i]->GetBoundaryConditionType()
2316  {
2317  SpatialDomains::NeumannBCShPtr bcPtr = boost::static_pointer_cast<
2319  m_bndConditions[i]);
2320  string filebcs = bcPtr->m_filename;
2321  if (filebcs != "")
2322  {
2323  ExtractFileBCs(filebcs, bcPtr->m_comm, varName, locExpList);
2324  }
2325  else
2326  {
2327  LibUtilities::Equation condition =
2328  boost::static_pointer_cast<
2330  (m_bndConditions[i])->
2331  m_neumannCondition;
2332  condition.Evaluate(x0, x1, x2, time,
2333  locExpList->UpdatePhys());
2334  }
2335 
2336  locExpList->IProductWRTBase(
2337  locExpList->GetPhys(),
2338  locExpList->UpdateCoeffs());
2339  }
2340  else if (m_bndConditions[i]->GetBoundaryConditionType()
2342  {
2343  SpatialDomains::RobinBCShPtr bcPtr = boost::static_pointer_cast<
2345  (m_bndConditions[i]);
2346  string filebcs = bcPtr->m_filename;
2347 
2348  if (filebcs != "")
2349  {
2350  ExtractFileBCs(filebcs, bcPtr->m_comm, varName, locExpList);
2351  }
2352  else
2353  {
2354  LibUtilities::Equation condition =
2355  boost::static_pointer_cast<
2357  (m_bndConditions[i])->
2358  m_robinFunction;
2359  condition.Evaluate(x0, x1, x2, time,
2360  locExpList->UpdatePhys());
2361  }
2362 
2363  locExpList->IProductWRTBase(
2364  locExpList->GetPhys(),
2365  locExpList->UpdateCoeffs());
2366  }
2367  else
2368  {
2369  ASSERTL0(false, "This type of BC not implemented yet");
2370  }
2371  }
2372  }
2373  }
2374  } // end of namespace
2375 } //end of namespace
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:914
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:1938
virtual void v_ExtractTracePhys(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This method extracts the trace (edges in 2D) from the field inarray and puts the values in outarray...
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1537
void FindPeriodicEdges(const SpatialDomains::BoundaryConditions &bcs, const std::string &variable)
Determine the periodic edges and vertices for the given graph.
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const boost::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:2011
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:2076
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1477
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2233
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
bool isLocal
Flag specifying if this entity is local to this partition.
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
Definition: ExpList.h:2084
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:248
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:46
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
boost::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: QuadGeom.h:54
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2067
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
STL namespace.
PeriodicMap m_periodicEdges
A map which identifies pairs of periodic edges.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1691
std::map< int, std::vector< unsigned int > > BndRegionOrdering
Definition: MeshPartition.h:54
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2286
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1015
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.
std::map< int, std::vector< unsigned int > > CompositeOrdering
Definition: MeshPartition.h:53
boost::shared_ptr< NeumannBoundaryCondition > NeumannBCShPtr
Definition: Conditions.h:221
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:998
boost::shared_ptr< MeshGraph2D > MeshGraph2DSharedPtr
Definition: MeshGraph2D.h:238
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2046
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
Search through the edge expansions and identify which ones have Robin/Mixed type boundary conditions...
boost::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:1294
boost::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Gauss Radau pinned at x=-1, .
Definition: PointsType.h:59
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Principle Orthogonal Functions .
Definition: BasisType.h:47
boost::shared_ptr< DirichletBoundaryCondition > DirichletBCShPtr
Definition: Conditions.h:220
bool SameTypeOfBoundaryConditions(const DisContField2D &In)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1047
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
DisContField2D()
Default constructor.
The base class for all shapes.
Definition: StdExpansion.h:69
boost::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2143
LocTraceToTraceMapSharedPtr m_locTraceToTraceMap
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
An object which contains the discretised boundary conditions.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1036
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)
boost::shared_ptr< ExpList1D > ExpList1DSharedPtr
Shared pointer to an ExpList1D object.
Definition: ExpList1D.h:50
virtual ~DisContField2D()
Default destructor.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:227
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:60
std::set< int > m_boundaryEdges
A set storing the global IDs of any boundary edges.
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1024
NekDouble Evaluate() const
Definition: Equation.h:102
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:976
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:217
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
static std::string npts
Definition: InputFld.cpp:43
Principle Orthogonal Functions .
Definition: BasisType.h:46
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:972
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:969
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Defines a list of flags.
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:226
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
bool IsLeftAdjacentEdge(const int n, const int e)
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
static const NekDouble kNekUnsetDouble
Describe a linear system.
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:115
StdRegions::Orientation orient
Orientation of entity within higher dimensional entity.
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:2037
void SetUpDG(const std::string="DefaultVar")
Set up all DG member variables and maps.
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry2D.h:59
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:343
boost::shared_ptr< ElementEdgeVector > ElementEdgeVectorSharedPtr
Definition: MeshGraph.h:134
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:966
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:55
2D geometry information
Definition: Geometry2D.h:65
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1641
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
Definition: ExpList2D.h:60
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1485
virtual void v_GetFwdBwdTracePhys(const Array< OneD, const NekDouble > &field, Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
This method extracts the "forward" and "backward" trace data from the array field and puts the data i...
boost::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
void EvaluateHDGPostProcessing(Array< OneD, NekDouble > &outarray)
Evaluate HDG post-processing to increase polynomial order of solution.
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:3059
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:52
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:238
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Calculates the result of the multiplication of a global matrix of type specified by mkey with a vecto...
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
boost::shared_ptr< Basis > BasisSharedPtr
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:737
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph2D, const SpatialDomains::BoundaryConditions &bcs, const std::string &variable, const bool DeclareCoeffPhysArrays=true)
This function discretises the boundary conditions by setting up a list of one-dimensional boundary ex...
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Set up a list of element IDs and edge IDs that link to the boundary conditions.
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
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:528
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
NekDouble L2_DGDeriv(const int dir, const Array< OneD, const NekDouble > &soln)
Calculate the error of the derivative using the consistent DG evaluation of .
Describes the specification for a Basis.
Definition: Basis.h:50
Array< OneD, DataType > & GetPtr()
Definition: NekVector.cpp:230
boost::shared_ptr< RobinBoundaryCondition > RobinBCShPtr
Definition: Conditions.h:222
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:52
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:299
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:49