Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DisContField1D.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File DisContField1D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Field definition for 1D domain with boundary
33 // conditions using LDG-H
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
38 #include <StdRegions/StdSegExp.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  /**
49  * @class DisContField1D
50  * This class augments the list of local expansions inherited from
51  * ExpList1D with boundary conditions. Inter-element boundaries are
52  * handled using an discontinuous Galerkin scheme.
53  */
54 
55  /**
56  * Constructs an empty expansion list with no boundary conditions.
57  */
58  DisContField1D::DisContField1D():
59  ExpList1D(),
60  m_bndCondExpansions(),
61  m_bndConditions()
62  {
63  }
64 
65 
66  /**
67  * An expansion list for the boundary expansions is generated first for
68  * the field. These are subsequently evaluated for time zero. The trace
69  * map is then constructed.
70  * @param graph1D A mesh, containing information about the domain
71  * and the spectral/hp element expansions.
72  * @param bcs Information about the enforced boundary
73  * conditions.
74  * @param variable The session variable associated with the
75  * boundary conditions to enforce.
76  * @param solnType Type of global system to use.
77  */
81  const std::string &variable,
82  const bool SetUpJustDG,
83  const Collections::ImplementationType ImpType)
84  : ExpList1D(pSession,graph1D,true, ImpType),
85  m_bndCondExpansions(),
86  m_bndConditions()
87  {
88  if (variable.compare("DefaultVar") != 0)
89  {
91 
92  GenerateBoundaryConditionExpansion(graph1D,bcs,variable);
93  EvaluateBoundaryConditions(0.0, variable);
94  ApplyGeomInfo();
95  FindPeriodicVertices(bcs,variable);
96  }
97 
98  if(SetUpJustDG)
99  {
100  SetUpDG(variable);
101  }
102  else
103  {
104  int i;
105  Array<OneD, int> ElmtID, VertexID;
106  GetBoundaryToElmtMap(ElmtID, VertexID);
107 
108  for(i = 0; i < m_bndCondExpansions.num_elements(); ++i)
109  {
111  locExpList = m_bndCondExpansions[i];
112 
114  (*m_exp)[ElmtID[i]]->
115  as<LocalRegions::Expansion1D>();
117  locExpList->GetExp(0)->
118  as<LocalRegions::Expansion0D>();
119 
120  exp0d->SetAdjacentElementExp(VertexID[i], exp1d);
121  }
122 
124  }
125 
126  }
127 
128  void DisContField1D::SetUpDG(const std::string &variable)
129  {
130  // Check for multiple calls
132  {
133  return;
134  }
135 
137  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph1D>(
138  m_graph);
139 
141 
146  *m_exp,graph1D,
148 
149  m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
150 
152  AllocateSharedPtr(m_session, graph1D, trace, *this,
154  m_periodicVerts, variable);
155 
156  if (m_session->DefinesCmdLineArgument("verbose"))
157  {
158  m_traceMap->PrintStats(std::cout, variable);
159  }
160 
161  // Scatter trace points to 1D elements. For each element, we find
162  // the trace point associated to each vertex. The element then
163  // retains a pointer to the trace space points, to ensure
164  // uniqueness of normals when retrieving from two adjoining
165  // elements which do not lie in a plane.
166 
167  int ElmtPointGeom = 0;
168  int TracePointGeom = 0;
171  for (int i = 0; i < m_exp->size(); ++i)
172  {
173  exp1d = (*m_exp)[i]->as<LocalRegions::Expansion1D>();
174  for (int j = 0; j < exp1d->GetNverts(); ++j)
175  {
176  ElmtPointGeom = (exp1d->GetGeom1D())->GetVid(j);
177 
178  for (int k = 0; k < m_trace->GetExpSize(); ++k)
179  {
180  TracePointGeom = m_trace->GetExp(k)->GetGeom()->GetVid(0);
181 
182  if (TracePointGeom == ElmtPointGeom)
183  {
184  exp0d = m_trace->GetExp(k)->as<LocalRegions::Expansion0D>();
185  exp0d->SetAdjacentElementExp(j,exp1d);
186  break;
187  }
188  }
189  }
190  }
191 
193 
194  int cnt, n, e;
195 
196  // Identify boundary verts
197  for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
198  {
199  if (m_bndConditions[n]->GetBoundaryConditionType() !=
201  {
202  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
203  {
204  m_boundaryVerts.insert(
205  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
206  }
207  }
208  else
209  {
210  ASSERTL0(false,"Periodic verts need setting up");
211  }
212  cnt += m_bndCondExpansions[n]->GetExpSize();
213  }
214 
215  // Set up left-adjacent edge list.
216  m_leftAdjacentVerts.resize(2*((*m_exp).size()));
217 
218  // count size of trace
219  for (cnt = n = 0; n < m_exp->size(); ++n)
220  {
221  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v, ++cnt)
222  {
224  }
225  }
226 
227 
228  boost::unordered_map<int,pair<int,int> > perVertToExpMap;
229  boost::unordered_map<int,pair<int,int> >::iterator it2;
230  for (n = 0; n < m_exp->size(); ++n)
231  {
232  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
233  {
235  (*m_exp)[n]->GetGeom()->GetVid(v));
236 
237  if (it != m_periodicVerts.end())
238  {
239  perVertToExpMap[it->first] = make_pair(n,v);
240  }
241  }
242  }
243 
244 
245  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
246  for (n = 0; n < m_exp->size(); ++n)
247  {
248  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
249  {
250  int vertGeomId = (*m_exp)[n]->GetGeom()->GetVid(v);
251 
252  // Check to see if this face is periodic.
253  PeriodicMap::iterator it = m_periodicVerts.find(vertGeomId);
254 
255  if (it != m_periodicVerts.end())
256  {
257  const PeriodicEntity &ent = it->second[0];
258  it2 = perVertToExpMap.find(ent.id);
259 
260  if (it2 == perVertToExpMap.end())
261  {
262  if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
263  !ent.isLocal)
264  {
265  continue;
266  }
267  else
268  {
269  ASSERTL1(false, "Periodic vert not found!");
270  }
271  }
272 
273  int offset = m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
274  [n][v]->GetElmtId());
275  m_periodicFwdCopy.push_back(offset);
276 
277  int offset2 = m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
278  [it2->second.first]
279  [it2->second.second]->GetElmtId());
280  m_periodicBwdCopy.push_back(offset2);
281  }
282  }
283  }
284  }
285 
286 
287  bool DisContField1D::IsLeftAdjacentVertex(const int n, const int e)
288  {
291  m_traceMap->GetElmtToTrace()[n][e]->as<LocalRegions::Expansion0D>();
292 
293 
294  bool fwd = true;
295  if (traceEl->GetLeftAdjacentElementVertex () == -1 ||
296  traceEl->GetRightAdjacentElementVertex() == -1)
297  {
298  // Boundary edge (1 connected element). Do nothing in
299  // serial.
300  it = m_boundaryVerts.find(traceEl->GetElmtId());
301 
302  // If the edge does not have a boundary condition set on
303  // it, then assume it is a partition edge or periodic.
304  if (it == m_boundaryVerts.end())
305  {
306  int traceGeomId = traceEl->GetGeom0D()->GetGlobalID();
308  traceGeomId);
309 
310  if (pIt != m_periodicVerts.end() && !pIt->second[0].isLocal)
311  {
312  fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
313  }
314  else
315  {
316  int offset = m_trace->GetPhys_Offset(traceEl->GetElmtId());
317  fwd = m_traceMap->
318  GetTraceToUniversalMapUnique(offset) >= 0;
319  }
320  }
321  }
322  else if (traceEl->GetLeftAdjacentElementVertex () != -1 &&
323  traceEl->GetRightAdjacentElementVertex() != -1)
324  {
325  // Non-boundary edge (2 connected elements).
326  fwd = dynamic_cast<Nektar::StdRegions::StdExpansion*>
327  (traceEl->GetLeftAdjacentElementExp().get()) ==
328  (*m_exp)[n].get();
329  }
330  else
331  {
332  ASSERTL2(false, "Unconnected trace element!");
333  }
334 
335  return fwd;
336  }
337 
338 
339  // Given all boundary regions for the whole solution determine
340  // which ones (if any) are part of domain and ensure all other
341  // conditions are given as UserDefined Dirichlet.
343  const SpatialDomains::CompositeMap &domain,
345  const std::string &variable)
346  {
348 
350 
352  map<int,int> GeometryToRegionsMap;
353 
354  SpatialDomains::BoundaryRegionCollection::const_iterator it;
355 
357  = Allbcs.GetBoundaryRegions();
359  = Allbcs.GetBoundaryConditions();
360 
361  // Set up a map of all boundary regions
362  for(it = bregions.begin(); it != bregions.end(); ++it)
363  {
365  for (bregionIt = it->second->begin();
366  bregionIt != it->second->end(); bregionIt++)
367  {
368  // can assume that all regions only contain one point in 1D
369  // Really do not need loop above
370  int id = (*(bregionIt->second))[0]->GetGlobalID();
371  GeometryToRegionsMap[id] = it->first;
372  }
373  }
374 
376  map<int,SpatialDomains::GeometrySharedPtr> EndOfDomain;
377 
378  // Now find out which points in domain have only one vertex
379  for(domIt = domain.begin(); domIt != domain.end(); ++domIt)
380  {
381  SpatialDomains::Composite geomvector = domIt->second;
382  for(int i = 0; i < geomvector->size(); ++i)
383  {
384  for(int j = 0; j < 2; ++j)
385  {
386  int vid = (*geomvector)[i]->GetVid(j);
387  if(EndOfDomain.count(vid) == 0)
388  {
389  EndOfDomain[vid] = (*geomvector)[i]->GetVertex(j);
390  }
391  else
392  {
393  EndOfDomain.erase(vid);
394  }
395  }
396  }
397  }
398  ASSERTL1(EndOfDomain.size() == 2,"Did not find two ends of domain");
399 
401  int numNewBc = 1;
402  for(regIt = EndOfDomain.begin(); regIt != EndOfDomain.end(); ++regIt)
403  {
404  if(GeometryToRegionsMap.count(regIt->first) != 0) // Set up boundary condition up
405  {
406  map<int,int>::iterator iter = GeometryToRegionsMap.find(regIt->first);
407  ASSERTL1(iter != GeometryToRegionsMap.end(),"Failied to find GeometryToRegionMap");
408  int regionId = iter->second;
409  SpatialDomains::BoundaryRegionCollection::const_iterator bregionsIter = bregions.find(regionId);
410  ASSERTL1(bregionsIter != bregions.end(),"Failed to find boundary region");
411  SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
412  returnval->AddBoundaryRegions (regionId,breg);
413 
414  SpatialDomains::BoundaryConditionCollection::const_iterator bconditionsIter = bconditions.find(regionId);
415  ASSERTL1(bconditionsIter != bconditions.end(),"Failed to find boundary collection");
416  SpatialDomains::BoundaryConditionMapShPtr bcond = bconditionsIter->second;
417  returnval->AddBoundaryConditions(regionId,bcond);
418  }
419  else // Set up an undefined region.
420  {
422 
423  // Set up Composite (GemetryVector) to contain vertex and put into bRegion
425  gvec->push_back(regIt->second);
426  (*breg)[regIt->first] = gvec;
427 
428  returnval->AddBoundaryRegions(bregions.size()+numNewBc,breg);
429 
431 
432  // Set up just boundary condition for this variable.
434  (*bCondition)[variable] = notDefinedCondition;
435 
436  returnval->AddBoundaryConditions(bregions.size()+numNewBc,bCondition);
437  ++numNewBc;
438 
439  }
440  }
441 
442  return returnval;
443  }
444 
445  /**
446  * Constructor for use in multidomain computations where a
447  * domain list can be passed instead of graph1D
448  *
449  * @param domain Subdomain specified in the inputfile from
450  * which the DisContField1D is set up
451  */
453  const LibUtilities::SessionReaderSharedPtr &pSession,
454  const SpatialDomains::MeshGraphSharedPtr &graph1D,
455  const SpatialDomains::CompositeMap &domain,
456  const SpatialDomains::BoundaryConditions &Allbcs,
457  const std::string &variable,
458  bool SetToOneSpaceDimension,
459  const Collections::ImplementationType ImpType):
460  ExpList1D(pSession,graph1D,domain, true,variable,SetToOneSpaceDimension,ImpType),
461  m_bndCondExpansions(),
462  m_bndConditions()
463  {
464  if (variable.compare("DefaultVar") != 0)
465  {
466  SpatialDomains::BoundaryConditionsSharedPtr DomBCs = GetDomainBCs(domain,Allbcs,variable);
467 
468  GenerateBoundaryConditionExpansion(graph1D,*DomBCs,variable);
469  EvaluateBoundaryConditions(0.0, variable);
470  ApplyGeomInfo();
471  FindPeriodicVertices(*DomBCs,variable);
472  }
473 
474  SetUpDG(variable);
475  }
476 
477  /**
478  * Constructs a field as a copy of an existing field.
479  * @param In Existing DisContField1D object to copy.
480  */
482  ExpList1D(In),
483  m_bndCondExpansions(In.m_bndCondExpansions),
484  m_bndConditions(In.m_bndConditions),
485  m_globalBndMat(In.m_globalBndMat),
486  m_trace(In.m_trace),
487  m_traceMap(In.m_traceMap),
488  m_boundaryVerts(In.m_boundaryVerts),
489  m_periodicVerts(In.m_periodicVerts),
490  m_periodicFwdCopy(In.m_periodicFwdCopy),
491  m_periodicBwdCopy(In.m_periodicBwdCopy),
492  m_leftAdjacentVerts(In.m_leftAdjacentVerts)
493  {
494  }
495 
496 
497  /**
498  * Constructs a field as a copy of an existing explist1D field.
499  * @param In Existing ExpList1D object to copy.
500  */
502  ExpList1D(In)
503  {
504  }
505 
506  /**
507  *
508  */
510  {
511  }
512 
513 
514  /**
515  * Generate the boundary condition expansion list
516  * @param graph1D A mesh, containing information about the domain
517  * and the spectral/hp element expansions.
518  * @param bcs Information about the enforced boundary
519  * conditions.
520  * @param variable The session variable associated with the
521  * boundary conditions to enforce.
522  */
524  const SpatialDomains::MeshGraphSharedPtr &graph1D,
526  const std::string variable)
527  {
528  int cnt = 0;
529 
531  = bcs.GetBoundaryRegions();
533  = bcs.GetBoundaryConditions();
534  SpatialDomains::BoundaryRegionCollection::const_iterator it;
535 
536  // count the number of non-periodic boundary points
537  for (it = bregions.begin(); it != bregions.end(); ++it)
538  {
539  const SpatialDomains::BoundaryConditionShPtr boundaryCondition =
540  GetBoundaryCondition(bconditions, it->first, variable);
541  if (boundaryCondition->GetBoundaryConditionType() !=
543  {
545  for (bregionIt = it->second->begin();
546  bregionIt != it->second->end(); bregionIt++)
547  {
548  cnt += bregionIt->second->size();
549  }
550  }
551  }
552 
555 
558 
559  SetBoundaryConditionExpansion(graph1D,bcs,variable,
561  m_bndConditions);
562  }
563 
564 
565  /**
566  * @param graph1D A mesh containing information about the domain
567  * and the spectral/hp element expansion.
568  * @param bcs Information about the boundary conditions.
569  * @param variable Specifies the field.
570  * @param periodicVertices Map into which the list of periodic
571  * vertices is placed.
572  */
575  const std::string variable)
576  {
578  = bcs.GetBoundaryRegions();
580  = bcs.GetBoundaryConditions();
581 
583  = boost::dynamic_pointer_cast<
585  SpatialDomains::BoundaryRegionCollection::const_iterator it;
586 
588  m_session->GetComm()->GetRowComm();
589 
590  int i, region1ID, region2ID;
591 
593 
594  map<int,int> BregionToVertMap;
595 
596  // Construct list of all periodic Region and their global vertex on
597  // this process.
598  for (it = bregions.begin(); it != bregions.end(); ++it)
599  {
600  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
601 
602  if (locBCond->GetBoundaryConditionType()
604  {
605  continue;
606  }
607 
608  int id = (*(it->second->begin()->second))[0]->GetGlobalID();
609 
610  BregionToVertMap[it->first] = id;
611  }
612 
614  set<int> islocal;
615 
616  int n = vComm->GetSize();
617  int p = vComm->GetRank();
618 
619  Array<OneD, int> nregions(n, 0);
620  nregions[p] = BregionToVertMap.size();
621  vComm->AllReduce(nregions, LibUtilities::ReduceSum);
622 
623  int totRegions = Vmath::Vsum(n, nregions, 1);
624 
625  Array<OneD, int> regOffset(n, 0);
626 
627  for (i = 1; i < n; ++i)
628  {
629  regOffset[i] = regOffset[i-1] + nregions[i-1];
630  }
631 
632  Array<OneD, int> bregmap(totRegions, 0);
633  Array<OneD, int> bregid (totRegions, 0);
634  for(i = regOffset[p], iit = BregionToVertMap.begin();
635  iit != BregionToVertMap.end(); ++iit, ++i)
636  {
637  bregid [i] = iit->first;
638  bregmap[i] = iit->second;
639  islocal.insert(iit->first);
640  }
641 
642  vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
643  vComm->AllReduce(bregid, LibUtilities::ReduceSum);
644 
645  for (int i = 0; i < totRegions; ++i)
646  {
647  BregionToVertMap[bregid[i]] = bregmap[i];
648  }
649 
650  // Construct list of all periodic pairs local to this process.
651  for (it = bregions.begin(); it != bregions.end(); ++it)
652  {
653  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
654 
655  if (locBCond->GetBoundaryConditionType()
657  {
658  continue;
659  }
660 
661  // Identify periodic boundary region IDs.
662  region1ID = it->first;
663  region2ID = boost::static_pointer_cast<
665  locBCond)->m_connectedBoundaryRegion;
666 
667  ASSERTL0(BregionToVertMap.count(region1ID) != 0,
668  "Cannot determine vertex of region1ID");
669 
670  ASSERTL0(BregionToVertMap.count(region2ID) != 0,
671  "Cannot determine vertex of region2ID");
672 
673  PeriodicEntity ent(BregionToVertMap[region2ID],
675  islocal.count(region2ID) != 0);
676 
677  m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
678  }
679  }
680 
681 
682  /**
683  * @param graph1D A mesh containing information about the domain
684  * and the Spectral/hp element expansion.
685  * @param bcs Information about the boundary conditions.
686  * @param variable Specifies the field.
687  * @param bndCondExpansions Array of ExpList1D objects each
688  * containing a 1D spectral/hp element expansion
689  * on a single boundary region.
690  * @param bncConditions Array of BoundaryCondition objects which
691  * contain information about the boundary
692  * conditions on the different boundary regions.
693  */
695  const SpatialDomains::MeshGraphSharedPtr &graph1D,
697  const std::string variable,
699  &bndCondExpansions,
700  Array<OneD, SpatialDomains
701  ::BoundaryConditionShPtr> &bndConditions)
702  {
703  int k;
704  int cnt = 0;
705 
707  = bcs.GetBoundaryRegions();
709  = bcs.GetBoundaryConditions();
710  SpatialDomains::BoundaryRegionCollection::const_iterator it;
711 
715 
716  cnt = 0;
717  for (it = bregions.begin(); it != bregions.end(); ++it)
718  {
719  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
720 
721  if (locBCond->GetBoundaryConditionType() !=
723  {
725  for (bregionIt = it->second->begin();
726  bregionIt != it->second->end(); bregionIt++)
727  {
728  for (k = 0; k < bregionIt->second->size(); k++)
729  {
730  if((vert = boost::dynamic_pointer_cast
732  (*bregionIt->second)[k])))
733  {
734  locPointExp
737  bndCondExpansions[cnt] = locPointExp;
738  bndConditions[cnt++] = locBCond;
739  }
740  else
741  {
742  ASSERTL0(false,
743  "dynamic cast to a vertex failed");
744  }
745  }
746  }
747  }
748  }
749  }
750 
751  /**
752  *
753  */
755  const GlobalLinSysKey &mkey)
756  {
758  "Routine currently only tested for HybridDGHelmholtz");
759 
761  "Full matrix global systems are not supported for HDG "
762  "expansions");
763 
765  ==m_traceMap->GetGlobalSysSolnType(),
766  "The local to global map is not set up for the requested "
767  "solution type");
768 
769  GlobalLinSysSharedPtr glo_matrix;
770  GlobalLinSysMap::iterator matrixIter = m_globalBndMat->find(mkey);
771 
772  if (matrixIter == m_globalBndMat->end())
773  {
774  glo_matrix = GenGlobalBndLinSys(mkey,m_traceMap);
775  (*m_globalBndMat)[mkey] = glo_matrix;
776  }
777  else
778  {
779  glo_matrix = matrixIter->second;
780  }
781 
782  return glo_matrix;
783  }
784 
785 
787  {
788  if(m_negatedFluxNormal.size() == 0)
789  {
791  &elmtToTrace = m_traceMap->GetElmtToTrace();
792 
793  m_negatedFluxNormal.resize(2*GetExpSize());
794 
795  for(int i = 0; i < GetExpSize(); ++i)
796  {
797 
798  for(int v = 0; v < 2; ++v)
799  {
800 
802  elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
803 
804  if(vertExp->GetLeftAdjacentElementExp()->GetGeom()->GetGlobalID() != (*m_exp)[i]->GetGeom()->GetGlobalID())
805  {
806  m_negatedFluxNormal[2*i+v] = true;
807  }
808  else
809  {
810  m_negatedFluxNormal[2*i+v] = false;
811  }
812  }
813  }
814 
815  }
816 
817  return m_negatedFluxNormal;
818  }
819 
820  /**
821  * Generate the forward or backward state for each trace point.
822  * @param Fwd Forward state.
823  * @param Bwd Backward state.
824  */
827  {
828  v_GetFwdBwdTracePhys(m_phys,Fwd,Bwd);
829  }
830 
831 
832  /**
833  * @brief This method extracts the "forward" and "backward" trace data
834  * from the array @a field and puts the data into output vectors @a Fwd
835  * and @a Bwd.
836  *
837  * We first define the convention which defines "forwards" and
838  * "backwards". First an association is made between the vertex of each
839  * element and its corresponding vertex in the trace space using the
840  * mapping #m_traceMap. The element can either be left-adjacent or
841  * right-adjacent to this trace edge (see
842  * Expansion0D::GetLeftAdjacentElementExp). Boundary edges are never
843  * left-adjacent since elemental left-adjacency is populated first.
844  *
845  * If the element is left-adjacent we extract the vertex trace data from
846  * @a field into the forward trace space @a Fwd; otherwise, we place it
847  * in the backwards trace space @a Bwd. In this way, we form a unique
848  * set of trace normals since these are always extracted from
849  * left-adjacent elements.
850  *
851  * @param field is a NekDouble array which contains the 1D data
852  * from which we wish to extract the backward and forward
853  * orientated trace/edge arrays.
854  * @param Fwd The resulting forwards space.
855  * @param Bwd The resulting backwards space.
856  */
857 
859  const Array<OneD, const NekDouble> &field,
862  {
863  // Counter variables
864  int n, v;
865 
866  // Number of elements
867  int nElements = GetExpSize();
868 
869  // Initial index of each element
870  int phys_offset;
871 
873  &elmtToTrace = m_traceMap->GetElmtToTrace();
874 
875  // Set forward and backard state to zero
876  Vmath::Zero(Fwd.num_elements(), Fwd, 1);
877  Vmath::Zero(Bwd.num_elements(), Bwd, 1);
878 
879  int cnt;
880 
881  // Loop on the elements
882  for (cnt = n = 0; n < nElements; ++n)
883  {
884  // Set the offset of each element
885  phys_offset = GetPhys_Offset(n);
886 
887  for(v = 0; v < 2; ++v, ++cnt)
888  {
889  int offset = m_trace->GetPhys_Offset(elmtToTrace[n][v]->GetElmtId());
890 
891  if (m_leftAdjacentVerts[cnt])
892  {
893  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
894  Fwd[offset]);
895  }
896  else
897  {
898  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
899  Bwd[offset]);
900  }
901  }
902  }
903 
904  // Fill boundary conditions into missing elements.
905  int id = 0;
906 
907  for(cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
908  {
909  if (m_bndConditions[n]->GetBoundaryConditionType() ==
911  {
912  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
913  Bwd[id] = m_bndCondExpansions[n]->GetPhys()[0]; //this is not getting the correct value?
914  cnt++;
915  }
916  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
918  m_bndConditions[n]->GetBoundaryConditionType() ==
920  {
921  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[0]==0.0,
922  "Method not set up for non-zero Neumann "
923  "boundary condition");
924  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
925  Bwd[id] = Fwd[id];
926 
927  cnt++;
928  }
929  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
931  {
932  // Do nothing
933  }
934  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
936  {
937  ASSERTL0(false,
938  "Method not set up for this boundary condition.");
939  }
940  }
941 
942  // Copy any periodic boundary conditions.
943  for (n = 0; n < m_periodicFwdCopy.size(); ++n)
944  {
945  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
946  }
947 
948  // Do parallel exchange for forwards/backwards spaces.
949  m_traceMap->UniversalTraceAssemble(Fwd);
950  m_traceMap->UniversalTraceAssemble(Bwd);
951 
952  }
953 
954 
956  Array<OneD, NekDouble> &outarray)
957  {
958  ASSERTL1(m_physState == true,"local physical space is not true ");
959  v_ExtractTracePhys(m_phys, outarray);
960  }
961 
962  /**
963  * @brief This method extracts the trace (verts in 1D) from the field @a
964  * inarray and puts the values in @a outarray.
965  *
966  * It assumes the field is C0 continuous so that it can overwrite the
967  * edge data when visited by the two adjacent elements.
968  *
969  * @param inarray An array containing the 1D data from which we wish
970  * to extract the edge data.
971  * @param outarray The resulting edge information.
972  *
973  * This will not work for non-boundary expansions
974  */
976  const Array<OneD, const NekDouble> &inarray,
977  Array<OneD, NekDouble> &outarray)
978  {
979  // Loop over elemente and collect forward expansion
980  int nexp = GetExpSize();
981  int n,p,offset,phys_offset;
982 
983  ASSERTL1(outarray.num_elements() >= m_trace->GetExpSize(),
984  "input array is of insufficient length");
985 
986  for (n = 0; n < nexp; ++n)
987  {
988  phys_offset = GetPhys_Offset(n);
989 
990  for (p = 0; p < (*m_exp)[n]->GetNverts(); ++p)
991  {
992  offset = m_trace->GetPhys_Offset(
993  (m_traceMap->GetElmtToTrace())[n][p]->GetElmtId());
994  (*m_exp)[n]->GetVertexPhysVals(p, inarray + phys_offset,
995  outarray[offset]);
996  }
997  }
998  }
999 
1001  const Array<OneD, const NekDouble> &Fn,
1002  Array<OneD, NekDouble> &outarray)
1003  {
1004  int n,offset, t_offset;
1005 
1007  &elmtToTrace = m_traceMap->GetElmtToTrace();
1008 
1009  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
1010 
1011  for (n = 0; n < GetExpSize(); ++n)
1012  {
1013  // Number of coefficients on each element
1014  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1015 
1016  offset = GetCoeff_Offset(n);
1017 
1018  // Implementation for every points except Gauss points
1019  if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
1021  {
1022  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1023  if(negatedFluxNormal[2*n])
1024  {
1025  outarray[offset] -= Fn[t_offset];
1026  }
1027  else
1028  {
1029  outarray[offset] += Fn[t_offset];
1030  }
1031 
1032  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1033 
1034  if(negatedFluxNormal[2*n+1])
1035  {
1036  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] -= Fn[t_offset];
1037  }
1038  else
1039  {
1040  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] += Fn[t_offset];
1041  }
1042 
1043  }
1044  else
1045  {
1046 #if 0
1047  DNekMatSharedPtr m_Ixm;
1050  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1052  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1053 
1054  BASE = LibUtilities::BasisManager()[BS_k];
1055 
1056  Array<OneD, NekDouble> coords(3, 0.0);
1057 
1058  int j;
1059 
1060  for(p = 0; p < 2; ++p)
1061  {
1062  NekDouble vertnorm = 0.0;
1063  for (int i=0; i<((*m_exp)[n]->
1064  GetVertexNormal(p)).num_elements(); i++)
1065  {
1066  vertnorm += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
1067  coords[0] = vertnorm ;
1068  }
1069 
1070  t_offset = GetTrace()->GetPhys_Offset(n+p);
1071 
1072  if (vertnorm >= 0.0)
1073  {
1074  m_Ixm = BASE->GetI(coords);
1075 
1076 
1077  for (j = 0; j < e_ncoeffs; j++)
1078  {
1079  outarray[offset + j] +=
1080  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1081  }
1082  }
1083 
1084  if (vertnorm < 0.0)
1085  {
1086  m_Ixm = BASE->GetI(coords);
1087 
1088  for (j = 0; j < e_ncoeffs; j++)
1089  {
1090  outarray[offset + j] -=
1091  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1092  }
1093  }
1094  }
1095 #else
1096  int j;
1097  static DNekMatSharedPtr m_Ixm, m_Ixp;
1098  static int sav_ncoeffs = 0;
1099  if(!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
1100  {
1103  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1105  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1106 
1107  BASE = LibUtilities::BasisManager()[BS_k];
1108 
1109  Array<OneD, NekDouble> coords(1, 0.0);
1110 
1111  coords[0] = -1.0;
1112  m_Ixm = BASE->GetI(coords);
1113 
1114  coords[0] = 1.0;
1115  m_Ixp = BASE->GetI(coords);
1116 
1117  sav_ncoeffs = e_ncoeffs;
1118  }
1119 
1120  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1121  if(negatedFluxNormal[2*n])
1122  {
1123  for (j = 0; j < e_ncoeffs; j++)
1124  {
1125  outarray[offset + j] -=
1126  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1127  }
1128  }
1129  else
1130  {
1131  for (j = 0; j < e_ncoeffs; j++)
1132  {
1133  outarray[offset + j] +=
1134  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1135  }
1136  }
1137 
1138  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1139  if (negatedFluxNormal[2*n+1])
1140  {
1141  for (j = 0; j < e_ncoeffs; j++)
1142  {
1143  outarray[offset + j] -=
1144  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1145  }
1146  }
1147  else
1148  {
1149  for (j = 0; j < e_ncoeffs; j++)
1150  {
1151  outarray[offset + j] +=
1152  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1153  }
1154  }
1155 #endif
1156  }
1157  }
1158  }
1159 
1160 
1162  const Array<OneD, const NekDouble> &inarray,
1163  Array<OneD, NekDouble> &outarray,
1164  const FlagList &flags,
1165  const StdRegions::ConstFactorMap &factors,
1166  const StdRegions::VarCoeffMap &varcoeff,
1167  const Array<OneD, const NekDouble> &dirForcing,
1168  const bool PhysSpaceForcing)
1169  {
1170  int i,n,cnt,nbndry;
1171  int nexp = GetExpSize();
1173  DNekVec F(m_ncoeffs,f,eWrapper);
1174  Array<OneD,NekDouble> e_f, e_l;
1175 
1176  //----------------------------------
1177  // Setup RHS Inner product if required
1178  //----------------------------------
1179  if(PhysSpaceForcing)
1180  {
1181  IProductWRTBase(inarray,f);
1182  Vmath::Neg(m_ncoeffs,f,1);
1183  }
1184  else
1185  {
1186  Vmath::Smul(m_ncoeffs,-1.0,inarray,1,f,1);
1187  }
1188 
1189  //----------------------------------
1190  // Solve continuous Boundary System
1191  //----------------------------------
1192  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
1193  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
1194  int e_ncoeffs,id;
1195 
1196  GlobalMatrixKey HDGLamToUKey(
1199  factors,
1200  varcoeff);
1201 
1202  const DNekScalBlkMatSharedPtr &HDGLamToU =
1203  GetBlockMatrix(HDGLamToUKey);
1204 
1205  // Retrieve global trace space storage, \Lambda, from trace expansion
1207  (m_traceMap->GetNumLocalBndCoeffs());
1208 
1209 
1210  Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
1211  // Zero trace space
1212  Vmath::Zero(GloBndDofs,BndSol,1);
1213 
1214  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
1215  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
1216  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
1217 
1218  //----------------------------------
1219  // Evaluate Trace Forcing
1220  //----------------------------------
1221  // Determing <u_lam,f> terms using HDGLamToU matrix
1222  for (cnt = n = 0; n < nexp; ++n)
1223  {
1224  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1225 
1226  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1227  e_f = f+m_coeff_offset[n];
1228  e_l = loc_lambda + cnt;
1229 
1230  // use outarray as tmp space
1231  DNekVec Floc (nbndry, e_l, eWrapper);
1232  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
1233  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1234 
1235  cnt += nbndry;
1236  }
1237 
1238  // Assemble into global operator
1239  m_traceMap->AssembleBnd(loc_lambda,BndRhs);
1240 
1241  cnt = 0;
1242  // Copy Dirichlet boundary conditions into trace space
1243  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1244  {
1245  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1247  {
1248  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1249  BndSol[id] = m_bndCondExpansions[i]->GetCoeff(0);
1250  }
1251  else
1252  {
1253  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1254  BndRhs[id] += m_bndCondExpansions[i]->GetCoeff(0);
1255  }
1256  }
1257 
1258  //----------------------------------
1259  // Solve trace problem
1260  //----------------------------------
1261  if (GloBndDofs - NumDirBCs > 0)
1262  {
1264  m_traceMap,factors);
1266  LinSys->Solve(BndRhs,BndSol,m_traceMap);
1267  }
1268 
1269  //----------------------------------
1270  // Internal element solves
1271  //----------------------------------
1274  factors);
1275 
1276  const DNekScalBlkMatSharedPtr& InvHDGHelm =
1277  GetBlockMatrix(invHDGhelmkey);
1278  DNekVec out(m_ncoeffs,outarray,eWrapper);
1279  Vmath::Zero(m_ncoeffs,outarray,1);
1280 
1281  // get local trace solution from BndSol
1282  m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1283 
1284  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
1285  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
1286  }
1287 
1288  /**
1289  * Based on the expression \f$g(x,t)\f$ for the boundary conditions,
1290  * this function evaluates the boundary conditions for all boundaries
1291  * at time-level \a t.
1292  * @param time The time at which the boundary conditions
1293  * should be evaluated.
1294  * @param bndCondExpansions List of boundary expansions.
1295  * @param bndConditions Information about the boundary conditions.
1296  */
1298  const NekDouble time,
1299  const std::string varName,
1300  const NekDouble x2_in,
1301  const NekDouble x3_in)
1302  {
1303  int i;
1304 
1305  Array<OneD, NekDouble> x0(1);
1306  Array<OneD, NekDouble> x1(1);
1307  Array<OneD, NekDouble> x2(1);
1308 
1309  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1310  {
1311  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
1312  {
1313  m_bndCondExpansions[i]->GetCoords(x0, x1, x2);
1314 
1315  if (x2_in != NekConstants::kNekUnsetDouble && x3_in !=
1317  {
1318  x1[0] = x2_in;
1319  x2[0] = x3_in;
1320  }
1321 
1322  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1324  {
1325  m_bndCondExpansions[i]->SetCoeff(0,
1326  (boost::static_pointer_cast<SpatialDomains
1327  ::DirichletBoundaryCondition>(m_bndConditions[i])
1328  ->m_dirichletCondition).Evaluate(x0[0],x1[0],x2[0],time));
1329  m_bndCondExpansions[i]->SetPhys(0,m_bndCondExpansions[i]->GetCoeff(0));
1330  }
1331  else if (m_bndConditions[i]->GetBoundaryConditionType()
1333  {
1334  m_bndCondExpansions[i]->SetCoeff(0,
1335  (boost::static_pointer_cast<SpatialDomains
1336  ::NeumannBoundaryCondition>(m_bndConditions[i])
1337  ->m_neumannCondition).Evaluate(x0[0],x1[0],x2[0],time));
1338  }
1339  else if (m_bndConditions[i]->GetBoundaryConditionType()
1341  {
1342  m_bndCondExpansions[i]->SetCoeff(0,
1343  (boost::static_pointer_cast<SpatialDomains
1344  ::RobinBoundaryCondition>(m_bndConditions[i])
1345  ->m_robinFunction).Evaluate(x0[0],x1[0],x2[0],time));
1346 
1347  }
1348  else if (m_bndConditions[i]->GetBoundaryConditionType()
1350  {
1351  }
1352  else
1353  {
1354  ASSERTL0(false, "This type of BC not implemented yet");
1355  }
1356  }
1357  }
1358  }
1359 
1360  // Set up a list of element ids and edge ids that link to the
1361  // boundary conditions
1363  Array<OneD, int> &ElmtID, Array<OneD,int> &VertID)
1364  {
1365  map<int, int> VertGID;
1366  int i,n,id;
1367  int bid,cnt,Vid;
1368  int nbcs = m_bndConditions.num_elements();
1369 
1370  // make sure arrays are of sufficient length
1371  if (ElmtID.num_elements() != nbcs)
1372  {
1373  ElmtID = Array<OneD, int>(nbcs,-1);
1374  }
1375  else
1376  {
1377  fill(ElmtID.get(), ElmtID.get()+nbcs, -1);
1378  }
1379 
1380  if (VertID.num_elements() != nbcs)
1381  {
1382  VertID = Array<OneD, int>(nbcs);
1383  }
1384 
1385  // setup map of all global ids along boundary
1386  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1387  {
1388  Vid = m_bndCondExpansions[n]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
1389  VertGID[Vid] = cnt++;
1390  }
1391 
1392  // Loop over elements and find verts that match;
1394  for(cnt = n = 0; n < GetExpSize(); ++n)
1395  {
1396  exp = (*m_exp)[n];
1397  for(i = 0; i < exp->GetNverts(); ++i)
1398  {
1399  id = exp->GetGeom()->GetVid(i);
1400 
1401  if (VertGID.count(id) > 0)
1402  {
1403  bid = VertGID.find(id)->second;
1404  ASSERTL1(ElmtID[bid] == -1,"Edge already set");
1405  ElmtID[bid] = n;
1406  VertID[bid] = i;
1407  cnt ++;
1408  }
1409  }
1410  }
1411 
1412  ASSERTL1(cnt == nbcs,"Failed to visit all boundary condtiions");
1413  }
1414 
1416  boost::shared_ptr<ExpList> &result,
1417  const bool DeclareCoeffPhysArrays)
1418  {
1419  int n, cnt, nq;
1420  int offsetOld, offsetNew;
1421  std::vector<unsigned int> eIDs;
1422 
1423  Array<OneD, int> ElmtID,EdgeID;
1424  GetBoundaryToElmtMap(ElmtID,EdgeID);
1425 
1426  // Skip other boundary regions
1427  for (cnt = n = 0; n < i; ++n)
1428  {
1429  cnt += m_bndCondExpansions[n]->GetExpSize();
1430  }
1431 
1432  // Populate eIDs with information from BoundaryToElmtMap
1433  for (n = 0; n < m_bndCondExpansions[i]->GetExpSize(); ++n)
1434  {
1435  eIDs.push_back(ElmtID[cnt+n]);
1436  }
1437 
1438  // Create expansion list
1439  result =
1441  (*this, eIDs, DeclareCoeffPhysArrays);
1442 
1443  // Copy phys and coeffs to new explist
1444  if( DeclareCoeffPhysArrays)
1445  {
1446  Array<OneD, NekDouble> tmp1, tmp2;
1447  for (n = 0; n < result->GetExpSize(); ++n)
1448  {
1449  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
1450  offsetOld = GetPhys_Offset(ElmtID[cnt+n]);
1451  offsetNew = result->GetPhys_Offset(n);
1452  Vmath::Vcopy(nq, tmp1 = GetPhys()+ offsetOld, 1,
1453  tmp2 = result->UpdatePhys()+ offsetNew, 1);
1454 
1455  nq = GetExp(ElmtID[cnt+n])->GetNcoeffs();
1456  offsetOld = GetCoeff_Offset(ElmtID[cnt+n]);
1457  offsetNew = result->GetCoeff_Offset(n);
1458  Vmath::Vcopy(nq, tmp1 = GetCoeffs()+ offsetOld, 1,
1459  tmp2 = result->UpdateCoeffs()+ offsetNew, 1);
1460  }
1461  }
1462  }
1463 
1464  /**
1465  * @brief Reset this field, so that geometry information can be updated.
1466  */
1468  {
1469  ExpList::v_Reset();
1470 
1471  // Reset boundary condition expansions.
1472  for (int n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1473  {
1474  m_bndCondExpansions[n]->Reset();
1475  }
1476  }
1477 
1478  /**
1479  * Search through the edge expansions and identify which ones
1480  * have Robin/Mixed type boundary conditions. If find a Robin
1481  * boundary then store the edge id of the boundary condition
1482  * and the array of points of the physical space boundary
1483  * condition which are hold the boundary condition primitive
1484  * variable coefficient at the quatrature points
1485  *
1486  * \return std map containing the robin boundary condition
1487  * info using a key of the element id
1488  *
1489  * There is a next member to allow for more than one Robin
1490  * boundary condition per element
1491  */
1492  map<int, RobinBCInfoSharedPtr> DisContField1D::v_GetRobinBCInfo(void)
1493  {
1494  int i;
1495  map<int, RobinBCInfoSharedPtr> returnval;
1496  Array<OneD, int> ElmtID,VertID;
1497  GetBoundaryToElmtMap(ElmtID,VertID);
1498 
1499  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1500  {
1501  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1503  {
1504  int elmtid;
1505 
1506  Array<OneD, NekDouble> x0(1);
1507  Array<OneD, NekDouble> x1(1);
1508  Array<OneD, NekDouble> x2(1);
1509  Array<OneD, NekDouble> coeffphys(1);
1510 
1511  m_bndCondExpansions[i]->GetCoords(x0, x1, x2);
1512 
1513  coeffphys[0] = (boost::static_pointer_cast<SpatialDomains
1515  ->m_robinPrimitiveCoeff).Evaluate(x0[0],x1[0],x2[0],0.0);
1516 
1517  RobinBCInfoSharedPtr rInfo =
1519  AllocateSharedPtr(VertID[i],coeffphys);
1520 
1521  elmtid = ElmtID[i];
1522  // make link list if necessary (not likely in
1523  // 1D but needed in 2D & 3D)
1524  if(returnval.count(elmtid) != 0)
1525  {
1526  rInfo->next = returnval.find(elmtid)->second;
1527  }
1528  returnval[elmtid] = rInfo;
1529  }
1530  }
1531 
1532  return returnval;
1533  }
1534  } // end of namespace
1535 } //end of namespace
virtual void v_GetBndElmtExpansion(int i, boost::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:914
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
Definition: ExpList.h:1938
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1537
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:2076
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
Evaluate all boundary conditions at a given time..
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1477
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2233
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fn, Array< OneD, NekDouble > &outarray)
bool isLocal
Flag specifying if this entity is local to this partition.
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n...
Definition: ExpList.h:2084
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:248
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Array< OneD, SpatialDomains::BoundaryConditionShPtr > m_bndConditions
An array which contains the information about the boundary condition on the different boundary region...
boost::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo()
void FindPeriodicVertices(const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Generate a associative map of periodic vertices in a mesh.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2067
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
STL namespace.
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
void SetUpDG(const std::string &variable)
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2286
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1015
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:1524
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2046
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
boost::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:49
boost::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:1294
void SetAdjacentElementExp(int vertex, Expansion1DSharedPtr &v)
Definition: Expansion0D.h:115
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:49
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
BasisManagerT & BasisManager(void)
GlobalLinSysSharedPtr GetGlobalBndLinSys(const GlobalLinSysKey &mkey)
For a given key, returns the associated global linear system.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1047
void SetBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable, Array< OneD, MultiRegions::ExpListSharedPtr > &bndCondExpansions, Array< OneD, SpatialDomains::BoundaryConditionShPtr > &bndConditions)
Populates the list of boundary condition expansions.
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
The base class for all shapes.
Definition: StdExpansion.h:69
boost::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2143
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1036
DisContField1D()
Default constructor.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:227
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1024
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:976
std::map< int, Composite >::const_iterator CompositeMapConstIter
Definition: MeshGraph.h:117
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:217
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
This class is the abstraction of a global discontinuous two- dimensional spectral/hp element expansio...
NekMatrix< InnerMatrixType, BlockMatrixTag > Transpose(NekMatrix< InnerMatrixType, BlockMatrixTag > &rhs)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:972
void GenerateBoundaryConditionExpansion(const SpatialDomains::MeshGraphSharedPtr &graph1D, const SpatialDomains::BoundaryConditions &bcs, const std::string variable)
Discretises the boundary conditions.
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:969
AssemblyMapDGSharedPtr m_traceMap
Local to global DG mapping for trace space.
Defines a specification for a set of points.
Definition: Points.h:58
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Defines a list of flags.
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:226
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:114
static const NekDouble kNekUnsetDouble
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Solve the Helmholtz equation.
Describe a linear system.
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:115
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:2037
bool IsLeftAdjacentVertex(const int n, const int e)
GlobalLinSysMapShPtr m_globalBndMat
Global boundary matrix.
boost::shared_ptr< ExpList0D > ExpList0DSharedPtr
Shared pointer to an ExpList0D object.
Definition: ExpList0D.h:54
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
static AssemblyMapSharedPtr NullAssemblyMapSharedPtr
Definition: AssemblyMap.h:55
boost::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:225
std::vector< bool > & GetNegatedFluxNormal(void)
Array< OneD, MultiRegions::ExpListSharedPtr > m_bndCondExpansions
Discretised boundary conditions.
virtual void v_Reset()
Reset this field, so that geometry information can be updated.
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:2015
boost::shared_ptr< MeshGraph1D > MeshGraph1DSharedPtr
Definition: MeshGraph1D.h:86
std::set< int > m_boundaryVerts
A set storing the global IDs of any boundary edges.
boost::shared_ptr< BoundaryConditions > BoundaryConditionsSharedPtr
Definition: Conditions.h:284
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1641
This class is the abstraction of a one-dimensional multi-elemental expansions which is merely a colle...
Definition: ExpList1D.h:61
ExpListSharedPtr m_trace
Trace space storage for points between elements.
std::vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:3059
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:52
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:238
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
boost::shared_ptr< Basis > BasisSharedPtr
PeriodicMap m_periodicVerts
A map which identifies groups of periodic vertices.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:737
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:442
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
boost::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:215
Describes the specification for a Basis.
Definition: Basis.h:50
boost::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:60
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &VertID)