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