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