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