Nektar++
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 namespace Nektar
43 {
44  namespace MultiRegions
45  {
46  /**
47  * @class DisContField1D
48  * This class augments the list of local expansions inherited from
49  * ExpList1D with boundary conditions. Inter-element boundaries are
50  * handled using an discontinuous Galerkin scheme.
51  */
52 
53  /**
54  * Constructs an empty expansion list with no boundary conditions.
55  */
57  ExpList1D(),
58  m_bndCondExpansions(),
59  m_bndConditions()
60  {
61  }
62 
63 
64  /**
65  * An expansion list for the boundary expansions is generated first for
66  * the field. These are subsequently evaluated for time zero. The trace
67  * map is then constructed.
68  * @param graph1D A mesh, containing information about the domain
69  * and the spectral/hp element expansions.
70  * @param bcs Information about the enforced boundary
71  * conditions.
72  * @param variable The session variable associated with the
73  * boundary conditions to enforce.
74  * @param solnType Type of global system to use.
75  */
79  const std::string &variable,
80  const bool SetUpJustDG)
81  : ExpList1D(pSession,graph1D),
82  m_bndCondExpansions(),
83  m_bndConditions()
84  {
86 
87  GenerateBoundaryConditionExpansion(graph1D,bcs,variable);
88  EvaluateBoundaryConditions(0.0, variable);
89  ApplyGeomInfo();
90  FindPeriodicVertices(bcs,variable);
91 
92  if(SetUpJustDG)
93  {
94  SetUpDG(variable);
95  }
96 
97  }
98 
99  void DisContField1D::SetUpDG(const std::string &variable)
100  {
101  // Check for multiple calls
103  {
104  return;
105  }
106 
108  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph1D>(
109  m_graph);
110 
112 
117  *m_exp,graph1D,
119 
120  m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
121 
123  AllocateSharedPtr(m_session, graph1D, trace, *this,
125  m_periodicVerts, variable);
126 
127  // Scatter trace points to 1D elements. For each element, we find
128  // the trace point associated to each vertex. The element then
129  // retains a pointer to the trace space points, to ensure
130  // uniqueness of normals when retrieving from two adjoining
131  // elements which do not lie in a plane.
132 
133  int ElmtPointGeom = 0;
134  int TracePointGeom = 0;
137  for (int i = 0; i < m_exp->size(); ++i)
138  {
139  exp1d = (*m_exp)[i]->as<LocalRegions::Expansion1D>();
140  for (int j = 0; j < exp1d->GetNverts(); ++j)
141  {
142  ElmtPointGeom = (exp1d->GetGeom1D())->GetVid(j);
143 
144  for (int k = 0; k < m_trace->GetExpSize(); ++k)
145  {
146  TracePointGeom = m_trace->GetExp(k)->GetGeom()->GetVid(0);
147 
148  if (TracePointGeom == ElmtPointGeom)
149  {
150  exp0d = m_trace->GetExp(k)->as<LocalRegions::Expansion0D>();
151  exp0d->SetAdjacentElementExp(j,exp1d);
152  break;
153  }
154  }
155  }
156  }
157 
159 
160  int cnt, n, e;
161 
162  // Identify boundary verts
163  for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
164  {
165  if (m_bndConditions[n]->GetBoundaryConditionType() !=
167  {
168  for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
169  {
170  m_boundaryVerts.insert(
171  m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e));
172  }
173  }
174  else
175  {
176  ASSERTL0(false,"Periodic verts need setting up");
177  }
178  cnt += m_bndCondExpansions[n]->GetExpSize();
179  }
180 
181  // Set up left-adjacent edge list.
182  m_leftAdjacentVerts.resize(2*((*m_exp).size()));
183 
184  // count size of trace
185  for (cnt = n = 0; n < m_exp->size(); ++n)
186  {
187  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v, ++cnt)
188  {
190  }
191  }
192 
193 
194  boost::unordered_map<int,pair<int,int> > perVertToExpMap;
195  boost::unordered_map<int,pair<int,int> >::iterator it2;
196  for (n = 0; n < m_exp->size(); ++n)
197  {
198  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
199  {
201  (*m_exp)[n]->GetGeom()->GetVid(v));
202 
203  if (it != m_periodicVerts.end())
204  {
205  perVertToExpMap[it->first] = make_pair(n,v);
206  }
207  }
208  }
209 
210 
211  // Set up mapping to copy Fwd of periodic bcs to Bwd of other edge.
212  for (n = 0; n < m_exp->size(); ++n)
213  {
214  for (int v = 0; v < (*m_exp)[n]->GetNverts(); ++v)
215  {
216  int vertGeomId = (*m_exp)[n]->GetGeom()->GetVid(v);
217 
218  // Check to see if this face is periodic.
219  PeriodicMap::iterator it = m_periodicVerts.find(vertGeomId);
220 
221  if (it != m_periodicVerts.end())
222  {
223  const PeriodicEntity &ent = it->second[0];
224  it2 = perVertToExpMap.find(ent.id);
225 
226  if (it2 == perVertToExpMap.end())
227  {
228  if (m_session->GetComm()->GetRowComm()->GetSize() > 1 &&
229  !ent.isLocal)
230  {
231  continue;
232  }
233  else
234  {
235  ASSERTL1(false, "Periodic vert not found!");
236  }
237  }
238 
239  int offset = m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
240  [n][v]->GetElmtId());
241  m_periodicFwdCopy.push_back(offset);
242 
243  int offset2 = m_trace->GetPhys_Offset((m_traceMap->GetElmtToTrace())
244  [it2->second.first]
245  [it2->second.second]->GetElmtId());
246  m_periodicBwdCopy.push_back(offset2);
247  }
248  }
249  }
250  }
251 
252 
253  bool DisContField1D::IsLeftAdjacentVertex(const int n, const int e)
254  {
257  m_traceMap->GetElmtToTrace()[n][e]->as<LocalRegions::Expansion0D>();
258 
259 
260  bool fwd = true;
261  if (traceEl->GetLeftAdjacentElementVertex () == -1 ||
262  traceEl->GetRightAdjacentElementVertex() == -1)
263  {
264  // Boundary edge (1 connected element). Do nothing in
265  // serial.
266  it = m_boundaryVerts.find(traceEl->GetElmtId());
267 
268  // If the edge does not have a boundary condition set on
269  // it, then assume it is a partition edge or periodic.
270  if (it == m_boundaryVerts.end())
271  {
272  int traceGeomId = traceEl->GetGeom0D()->GetGlobalID();
274  traceGeomId);
275 
276  if (pIt != m_periodicVerts.end() && !pIt->second[0].isLocal)
277  {
278  fwd = traceGeomId == min(traceGeomId,pIt->second[0].id);
279  }
280  else
281  {
282  int offset = m_trace->GetPhys_Offset(traceEl->GetElmtId());
283  fwd = m_traceMap->
284  GetTraceToUniversalMapUnique(offset) >= 0;
285  }
286  }
287  }
288  else if (traceEl->GetLeftAdjacentElementVertex () != -1 &&
289  traceEl->GetRightAdjacentElementVertex() != -1)
290  {
291  // Non-boundary edge (2 connected elements).
292  fwd = dynamic_cast<Nektar::StdRegions::StdExpansion*>
293  (traceEl->GetLeftAdjacentElementExp().get()) ==
294  (*m_exp)[n].get();
295  }
296  else
297  {
298  ASSERTL2(false, "Unconnected trace element!");
299  }
300 
301  return fwd;
302  }
303 
304 
305  // Given all boundary regions for the whole solution determine
306  // which ones (if any) are part of domain and ensure all other
307  // conditions are given as UserDefined Dirichlet.
309  const SpatialDomains::CompositeMap &domain,
311  const std::string &variable)
312  {
314 
316 
318  map<int,int> GeometryToRegionsMap;
319 
320  SpatialDomains::BoundaryRegionCollection::const_iterator it;
321 
323  = Allbcs.GetBoundaryRegions();
325  = Allbcs.GetBoundaryConditions();
326 
327  // Set up a map of all boundary regions
328  for(it = bregions.begin(); it != bregions.end(); ++it)
329  {
331  for (bregionIt = it->second->begin();
332  bregionIt != it->second->end(); bregionIt++)
333  {
334  // can assume that all regions only contain one point in 1D
335  // Really do not need loop above
336  int id = (*(bregionIt->second))[0]->GetGlobalID();
337  GeometryToRegionsMap[id] = it->first;
338  }
339  }
340 
342  map<int,SpatialDomains::GeometrySharedPtr> EndOfDomain;
343 
344  // Now find out which points in domain have only one vertex
345  for(domIt = domain.begin(); domIt != domain.end(); ++domIt)
346  {
347  SpatialDomains::Composite geomvector = domIt->second;
348  for(int i = 0; i < geomvector->size(); ++i)
349  {
350  for(int j = 0; j < 2; ++j)
351  {
352  int vid = (*geomvector)[i]->GetVid(j);
353  if(EndOfDomain.count(vid) == 0)
354  {
355  EndOfDomain[vid] = (*geomvector)[i]->GetVertex(j);
356  }
357  else
358  {
359  EndOfDomain.erase(vid);
360  }
361  }
362  }
363  }
364  ASSERTL1(EndOfDomain.size() == 2,"Did not find two ends of domain");
365 
367  int numNewBc = 1;
368  for(regIt = EndOfDomain.begin(); regIt != EndOfDomain.end(); ++regIt)
369  {
370  if(GeometryToRegionsMap.count(regIt->first) != 0) // Set up boundary condition up
371  {
372  map<int,int>::iterator iter = GeometryToRegionsMap.find(regIt->first);
373  ASSERTL1(iter != GeometryToRegionsMap.end(),"Failied to find GeometryToRegionMap");
374  int regionId = iter->second;
375  SpatialDomains::BoundaryRegionCollection::const_iterator bregionsIter = bregions.find(regionId);
376  ASSERTL1(bregionsIter != bregions.end(),"Failed to find boundary region");
377  SpatialDomains::BoundaryRegionShPtr breg = bregionsIter->second;
378  returnval->AddBoundaryRegions (regionId,breg);
379 
380  SpatialDomains::BoundaryConditionCollection::const_iterator bconditionsIter = bconditions.find(regionId);
381  ASSERTL1(bconditionsIter != bconditions.end(),"Failed to find boundary collection");
382  SpatialDomains::BoundaryConditionMapShPtr bcond = bconditionsIter->second;
383  returnval->AddBoundaryConditions(regionId,bcond);
384  }
385  else // Set up an undefined region.
386  {
388 
389  // Set up Composite (GemetryVector) to contain vertex and put into bRegion
391  gvec->push_back(regIt->second);
392  (*breg)[regIt->first] = gvec;
393 
394  returnval->AddBoundaryRegions(bregions.size()+numNewBc,breg);
395 
397 
398  // Set up just boundary condition for this variable.
400  (*bCondition)[variable] = notDefinedCondition;
401 
402  returnval->AddBoundaryConditions(bregions.size()+numNewBc,bCondition);
403  ++numNewBc;
404 
405  }
406  }
407 
408  return returnval;
409  }
410 
411  /**
412  * Constructor for use in multidomain computations where a
413  * domain list can be passed instead of graph1D
414  *
415  * @param domain Subdomain specified in the inputfile from
416  * which the DisContField1D is set up
417  */
419  const LibUtilities::SessionReaderSharedPtr &pSession,
420  const SpatialDomains::MeshGraphSharedPtr &graph1D,
421  const SpatialDomains::CompositeMap &domain,
422  const SpatialDomains::BoundaryConditions &Allbcs,
423  const std::string &variable,
424  bool SetToOneSpaceDimension):
425  ExpList1D(pSession,graph1D,domain, true,variable,SetToOneSpaceDimension),
426  m_bndCondExpansions(),
427  m_bndConditions()
428  {
429  SpatialDomains::BoundaryConditionsSharedPtr DomBCs = GetDomainBCs(domain,Allbcs,variable);
430 
431  GenerateBoundaryConditionExpansion(graph1D,*DomBCs,variable);
432  EvaluateBoundaryConditions(0.0, variable);
433  ApplyGeomInfo();
434  FindPeriodicVertices(*DomBCs,variable);
435 
436  SetUpDG(variable);
437  }
438 
439  /**
440  * Constructs a field as a copy of an existing field.
441  * @param In Existing DisContField1D object to copy.
442  */
444  ExpList1D(In),
445  m_bndCondExpansions(In.m_bndCondExpansions),
446  m_bndConditions(In.m_bndConditions),
447  m_globalBndMat(In.m_globalBndMat),
448  m_trace(In.m_trace),
449  m_traceMap(In.m_traceMap),
450  m_boundaryVerts(In.m_boundaryVerts),
451  m_periodicVerts(In.m_periodicVerts),
452  m_periodicFwdCopy(In.m_periodicFwdCopy),
453  m_periodicBwdCopy(In.m_periodicBwdCopy),
454  m_leftAdjacentVerts(In.m_leftAdjacentVerts)
455  {
456  }
457 
458 
459  /**
460  * Constructs a field as a copy of an existing explist1D field.
461  * @param In Existing ExpList1D object to copy.
462  */
464  ExpList1D(In)
465  {
466  }
467 
468  /**
469  *
470  */
472  {
473  }
474 
475 
476  /**
477  * Generate the boundary condition expansion list
478  * @param graph1D A mesh, containing information about the domain
479  * and the spectral/hp element expansions.
480  * @param bcs Information about the enforced boundary
481  * conditions.
482  * @param variable The session variable associated with the
483  * boundary conditions to enforce.
484  */
486  const SpatialDomains::MeshGraphSharedPtr &graph1D,
488  const std::string variable)
489  {
490  int cnt = 0;
491 
493  = bcs.GetBoundaryRegions();
495  = bcs.GetBoundaryConditions();
496  SpatialDomains::BoundaryRegionCollection::const_iterator it;
497 
498  // count the number of non-periodic boundary points
499  for (it = bregions.begin(); it != bregions.end(); ++it)
500  {
501  const SpatialDomains::BoundaryConditionShPtr boundaryCondition =
502  GetBoundaryCondition(bconditions, it->first, variable);
503  if (boundaryCondition->GetBoundaryConditionType() !=
505  {
507  for (bregionIt = it->second->begin();
508  bregionIt != it->second->end(); bregionIt++)
509  {
510  cnt += bregionIt->second->size();
511  }
512  }
513  }
514 
517 
520 
521  SetBoundaryConditionExpansion(graph1D,bcs,variable,
523  m_bndConditions);
524  }
525 
526 
527  /**
528  * @param graph1D A mesh containing information about the domain
529  * and the spectral/hp element expansion.
530  * @param bcs Information about the boundary conditions.
531  * @param variable Specifies the field.
532  * @param periodicVertices Map into which the list of periodic
533  * vertices is placed.
534  */
537  const std::string variable)
538  {
540  = bcs.GetBoundaryRegions();
542  = bcs.GetBoundaryConditions();
543 
545  = boost::dynamic_pointer_cast<
547  SpatialDomains::BoundaryRegionCollection::const_iterator it;
548 
550  m_session->GetComm()->GetRowComm();
551 
552  int i, region1ID, region2ID;
553 
555 
556  map<int,int> BregionToVertMap;
557 
558  // Construct list of all periodic Region and their global vertex on
559  // this process.
560  for (it = bregions.begin(); it != bregions.end(); ++it)
561  {
562  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
563 
564  if (locBCond->GetBoundaryConditionType()
566  {
567  continue;
568  }
569 
570  int id = (*(it->second->begin()->second))[0]->GetGlobalID();
571 
572  BregionToVertMap[it->first] = id;
573  }
574 
576  set<int> islocal;
577 
578  int n = vComm->GetSize();
579  int p = vComm->GetRank();
580 
581  Array<OneD, int> nregions(n, 0);
582  nregions[p] = BregionToVertMap.size();
583  vComm->AllReduce(nregions, LibUtilities::ReduceSum);
584 
585  int totRegions = Vmath::Vsum(n, nregions, 1);
586 
587  Array<OneD, int> regOffset(n, 0);
588 
589  for (i = 1; i < n; ++i)
590  {
591  regOffset[i] = regOffset[i-1] + nregions[i-1];
592  }
593 
594  Array<OneD, int> bregmap(totRegions, 0);
595  Array<OneD, int> bregid (totRegions, 0);
596  for(i = regOffset[p], iit = BregionToVertMap.begin();
597  iit != BregionToVertMap.end(); ++iit, ++i)
598  {
599  bregid [i] = iit->first;
600  bregmap[i] = iit->second;
601  islocal.insert(iit->first);
602  }
603 
604  vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
605  vComm->AllReduce(bregid, LibUtilities::ReduceSum);
606 
607  for (int i = 0; i < totRegions; ++i)
608  {
609  BregionToVertMap[bregid[i]] = bregmap[i];
610  }
611 
612  // Construct list of all periodic pairs local to this process.
613  for (it = bregions.begin(); it != bregions.end(); ++it)
614  {
615  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
616 
617  if (locBCond->GetBoundaryConditionType()
619  {
620  continue;
621  }
622 
623  // Identify periodic boundary region IDs.
624  region1ID = it->first;
625  region2ID = boost::static_pointer_cast<
627  locBCond)->m_connectedBoundaryRegion;
628 
629  ASSERTL0(BregionToVertMap.count(region1ID) != 0,
630  "Cannot determine vertex of region1ID");
631 
632  ASSERTL0(BregionToVertMap.count(region2ID) != 0,
633  "Cannot determine vertex of region2ID");
634 
635  PeriodicEntity ent(BregionToVertMap[region2ID],
637  islocal.count(region2ID) != 0);
638 
639  m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
640  }
641  }
642 
643 
644  /**
645  * @param graph1D A mesh containing information about the domain
646  * and the Spectral/hp element expansion.
647  * @param bcs Information about the boundary conditions.
648  * @param variable Specifies the field.
649  * @param bndCondExpansions Array of ExpList1D objects each
650  * containing a 1D spectral/hp element expansion
651  * on a single boundary region.
652  * @param bncConditions Array of BoundaryCondition objects which
653  * contain information about the boundary
654  * conditions on the different boundary regions.
655  */
657  const SpatialDomains::MeshGraphSharedPtr &graph1D,
659  const std::string variable,
661  &bndCondExpansions,
662  Array<OneD, SpatialDomains
663  ::BoundaryConditionShPtr> &bndConditions)
664  {
665  int k;
666  int cnt = 0;
667 
669  = bcs.GetBoundaryRegions();
671  = bcs.GetBoundaryConditions();
672  SpatialDomains::BoundaryRegionCollection::const_iterator it;
673 
677 
678  cnt = 0;
679  // list Dirichlet boundaries first
680  for (it = bregions.begin(); it != bregions.end(); ++it)
681  {
682  locBCond = GetBoundaryCondition(
683  bconditions, it->first, variable);
684 
685  if (locBCond->GetBoundaryConditionType() ==
687 
688  {
690  for (bregionIt = it->second->begin();
691  bregionIt != it->second->end(); bregionIt++)
692  {
693  for (k = 0; k < bregionIt->second->size(); k++)
694  {
695  if ((vert = boost::dynamic_pointer_cast
697  (*bregionIt->second)[k])))
698  {
699  locPointExp
702  bndCondExpansions[cnt] = locPointExp;
703  bndConditions[cnt++] = locBCond;
704  }
705  else
706  {
707  ASSERTL0(false,
708  "dynamic cast to a vertex failed");
709  }
710  }
711  }
712  }
713  } // end if Dirichlet
714 
715  // then, list the other (non-periodic) boundaries
716  for (it = bregions.begin(); it != bregions.end(); ++it)
717  {
718  locBCond = GetBoundaryCondition(bconditions, it->first, variable);
719 
720  switch(locBCond->GetBoundaryConditionType())
721  {
724  case SpatialDomains::eNotDefined: // presume this will be reused as Neuman, Robin or Dirichlet later
725  {
727  for (bregionIt = it->second->begin();
728  bregionIt != it->second->end(); bregionIt++)
729  {
730  for (k = 0; k < bregionIt->second->size(); k++)
731  {
732  if((vert = boost::dynamic_pointer_cast
734  (*bregionIt->second)[k])))
735  {
736  locPointExp
739  bndCondExpansions[cnt] = locPointExp;
740  bndConditions[cnt++] = locBCond;
741  }
742  else
743  {
744  ASSERTL0(false,
745  "dynamic cast to a vertex failed");
746  }
747  }
748  }
749  }
750  // do nothing for these types
753  break;
754  default:
755  ASSERTL0(false,"This type of BC not implemented yet");
756  break;
757  }
758  }
759  }
760 
761  /**
762  *
763  */
765  const GlobalLinSysKey &mkey)
766  {
768  "Routine currently only tested for HybridDGHelmholtz");
769 
771  "Full matrix global systems are not supported for HDG "
772  "expansions");
773 
775  ==m_traceMap->GetGlobalSysSolnType(),
776  "The local to global map is not set up for the requested "
777  "solution type");
778 
779  GlobalLinSysSharedPtr glo_matrix;
780  GlobalLinSysMap::iterator matrixIter = m_globalBndMat->find(mkey);
781 
782  if (matrixIter == m_globalBndMat->end())
783  {
784  glo_matrix = GenGlobalBndLinSys(mkey,m_traceMap);
785  (*m_globalBndMat)[mkey] = glo_matrix;
786  }
787  else
788  {
789  glo_matrix = matrixIter->second;
790  }
791 
792  return glo_matrix;
793  }
794 
795 
797  {
798  if(m_negatedFluxNormal.size() == 0)
799  {
801  &elmtToTrace = m_traceMap->GetElmtToTrace();
802 
803  m_negatedFluxNormal.resize(2*GetExpSize());
804 
805  for(int i = 0; i < GetExpSize(); ++i)
806  {
807 
808  for(int v = 0; v < 2; ++v)
809  {
810 
812  elmtToTrace[i][v]->as<LocalRegions::Expansion0D>();
813 
814  if(vertExp->GetLeftAdjacentElementExp()->GetGeom()->GetGlobalID() != (*m_exp)[i]->GetGeom()->GetGlobalID())
815  {
816  m_negatedFluxNormal[2*i+v] = true;
817  }
818  else
819  {
820  m_negatedFluxNormal[2*i+v] = false;
821  }
822  }
823  }
824 
825  }
826 
827  return m_negatedFluxNormal;
828  }
829 
830  /**
831  * Generate the forward or backward state for each trace point.
832  * @param Fwd Forward state.
833  * @param Bwd Backward state.
834  */
837  {
838  v_GetFwdBwdTracePhys(m_phys,Fwd,Bwd);
839  }
840 
841 
842  /**
843  * @brief This method extracts the "forward" and "backward" trace data
844  * from the array @a field and puts the data into output vectors @a Fwd
845  * and @a Bwd.
846  *
847  * We first define the convention which defines "forwards" and
848  * "backwards". First an association is made between the vertex of each
849  * element and its corresponding vertex in the trace space using the
850  * mapping #m_traceMap. The element can either be left-adjacent or
851  * right-adjacent to this trace edge (see
852  * Expansion0D::GetLeftAdjacentElementExp). Boundary edges are never
853  * left-adjacent since elemental left-adjacency is populated first.
854  *
855  * If the element is left-adjacent we extract the vertex trace data from
856  * @a field into the forward trace space @a Fwd; otherwise, we place it
857  * in the backwards trace space @a Bwd. In this way, we form a unique
858  * set of trace normals since these are always extracted from
859  * left-adjacent elements.
860  *
861  * @param field is a NekDouble array which contains the 1D data
862  * from which we wish to extract the backward and forward
863  * orientated trace/edge arrays.
864  * @param Fwd The resulting forwards space.
865  * @param Bwd The resulting backwards space.
866  */
867 
869  const Array<OneD, const NekDouble> &field,
872  {
873  // Expansion casts
877 
878  // Counter variables
879  int n, v;
880 
881  // Number of elements
882  int nElements = GetExpSize();
883 
884  // Initial index of each element
885  int phys_offset;
886 
888  &elmtToTrace = m_traceMap->GetElmtToTrace();
889 
890  // Basis shared pointer
892 
893  // Set forward and backard state to zero
894  Vmath::Zero(Fwd.num_elements(), Fwd, 1);
895  Vmath::Zero(Bwd.num_elements(), Bwd, 1);
896 
897  int cnt;
898 
899  // Loop on the elements
900  for (cnt = n = 0; n < nElements; ++n)
901  {
902  exp1D = (*m_exp)[n]->as<LocalRegions::Expansion1D>();
903 
904  // Set the offset of each element
905  phys_offset = GetPhys_Offset(n);
906 
907  Basis = (*m_exp)[n]->GetBasis(0);
908 
909  for(v = 0; v < 2; ++v, ++cnt)
910  {
911  int offset = m_trace->GetPhys_Offset(elmtToTrace[n][v]->GetElmtId());
912 
913  if (m_leftAdjacentVerts[cnt])
914  {
915  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
916  Fwd[offset]);
917  }
918  else
919  {
920  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
921  Bwd[offset]);
922  }
923  }
924  }
925 
926  // Fill boundary conditions into missing elements.
927  int id = 0;
928 
929  for(cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
930  {
931  if (m_bndConditions[n]->GetBoundaryConditionType() ==
933  {
934  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
935  Bwd[id] = m_bndCondExpansions[n]->GetPhys()[0]; //this is not getting the correct value?
936  cnt++;
937  }
938  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
940  m_bndConditions[n]->GetBoundaryConditionType() ==
942  {
943  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[0]==0.0,
944  "Method not set up for non-zero Neumann "
945  "boundary condition");
946  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
947  Bwd[id] = Fwd[id];
948 
949  cnt++;
950  }
951  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
953  {
954  // Do nothing
955  }
956  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
958  {
959  ASSERTL0(false,
960  "Method not set up for this boundary condition.");
961  }
962  }
963 
964  // Copy any periodic boundary conditions.
965  for (n = 0; n < m_periodicFwdCopy.size(); ++n)
966  {
967  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
968  }
969 
970  // Do parallel exchange for forwards/backwards spaces.
971  m_traceMap->UniversalTraceAssemble(Fwd);
972  m_traceMap->UniversalTraceAssemble(Bwd);
973 
974  }
975 
976 
978  Array<OneD, NekDouble> &outarray)
979  {
980  ASSERTL1(m_physState == true,"local physical space is not true ");
981  v_ExtractTracePhys(m_phys, outarray);
982  }
983 
984  /**
985  * @brief This method extracts the trace (verts in 1D) from the field @a
986  * inarray and puts the values in @a outarray.
987  *
988  * It assumes the field is C0 continuous so that it can overwrite the
989  * edge data when visited by the two adjacent elements.
990  *
991  * @param inarray An array containing the 1D data from which we wish
992  * to extract the edge data.
993  * @param outarray The resulting edge information.
994  *
995  * This will not work for non-boundary expansions
996  */
998  const Array<OneD, const NekDouble> &inarray,
999  Array<OneD, NekDouble> &outarray)
1000  {
1001  // Loop over elemente and collect forward expansion
1002  int nexp = GetExpSize();
1003  int n,p,offset,phys_offset;
1004 
1005  ASSERTL1(outarray.num_elements() >= m_trace->GetExpSize(),
1006  "input array is of insufficient length");
1007 
1008  for (n = 0; n < nexp; ++n)
1009  {
1010  phys_offset = GetPhys_Offset(n);
1011 
1012  for (p = 0; p < (*m_exp)[n]->GetNverts(); ++p)
1013  {
1014  offset = m_trace->GetPhys_Offset(
1015  (m_traceMap->GetElmtToTrace())[n][p]->GetElmtId());
1016  (*m_exp)[n]->GetVertexPhysVals(p, inarray + phys_offset,
1017  outarray[offset]);
1018  }
1019  }
1020  }
1021 
1023  const Array<OneD, const NekDouble> &Fn,
1024  Array<OneD, NekDouble> &outarray)
1025  {
1026  int n,offset, t_offset;
1027 
1029  &elmtToTrace = m_traceMap->GetElmtToTrace();
1030 
1031  // Basis shared pointer
1033 
1034  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
1035 
1036  for (n = 0; n < GetExpSize(); ++n)
1037  {
1038  // Basis definition on each element
1039  Basis = (*m_exp)[n]->GetBasis(0);
1040 
1041  // Number of coefficients on each element
1042  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1043 
1044  offset = GetCoeff_Offset(n);
1045 
1046  // Implementation for every points except Gauss points
1047  if (Basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
1048  {
1049  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1050  if(negatedFluxNormal[2*n])
1051  {
1052  outarray[offset] -= Fn[t_offset];
1053  }
1054  else
1055  {
1056  outarray[offset] += Fn[t_offset];
1057  }
1058 
1059  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1060 
1061  if(negatedFluxNormal[2*n+1])
1062  {
1063  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] -= Fn[t_offset];
1064  }
1065  else
1066  {
1067  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] += Fn[t_offset];
1068  }
1069 
1070  }
1071  else
1072  {
1073 #if 0
1074  DNekMatSharedPtr m_Ixm;
1077  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1079  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1080 
1081  BASE = LibUtilities::BasisManager()[BS_k];
1082 
1083  Array<OneD, NekDouble> coords(3, 0.0);
1084 
1085  int j;
1086 
1087  for(p = 0; p < 2; ++p)
1088  {
1089  NekDouble vertnorm = 0.0;
1090  for (int i=0; i<((*m_exp)[n]->
1091  GetVertexNormal(p)).num_elements(); i++)
1092  {
1093  vertnorm += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
1094  coords[0] = vertnorm ;
1095  }
1096 
1097  t_offset = GetTrace()->GetPhys_Offset(n+p);
1098 
1099  if (vertnorm >= 0.0)
1100  {
1101  m_Ixm = BASE->GetI(coords);
1102 
1103 
1104  for (j = 0; j < e_ncoeffs; j++)
1105  {
1106  outarray[offset + j] +=
1107  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1108  }
1109  }
1110 
1111  if (vertnorm < 0.0)
1112  {
1113  m_Ixm = BASE->GetI(coords);
1114 
1115  for (j = 0; j < e_ncoeffs; j++)
1116  {
1117  outarray[offset + j] -=
1118  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1119  }
1120  }
1121  }
1122 #else
1123  int j;
1124  static DNekMatSharedPtr m_Ixm, m_Ixp;
1125  static int sav_ncoeffs = 0;
1126  if(!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
1127  {
1130  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1132  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1133 
1134  BASE = LibUtilities::BasisManager()[BS_k];
1135 
1136  Array<OneD, NekDouble> coords(1, 0.0);
1137 
1138  coords[0] = -1.0;
1139  m_Ixm = BASE->GetI(coords);
1140 
1141  coords[0] = 1.0;
1142  m_Ixp = BASE->GetI(coords);
1143 
1144  sav_ncoeffs = e_ncoeffs;
1145  }
1146 
1147  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1148  if(negatedFluxNormal[2*n])
1149  {
1150  for (j = 0; j < e_ncoeffs; j++)
1151  {
1152  outarray[offset + j] -=
1153  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1154  }
1155  }
1156  else
1157  {
1158  for (j = 0; j < e_ncoeffs; j++)
1159  {
1160  outarray[offset + j] +=
1161  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1162  }
1163  }
1164 
1165  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1166  if (negatedFluxNormal[2*n+1])
1167  {
1168  for (j = 0; j < e_ncoeffs; j++)
1169  {
1170  outarray[offset + j] -=
1171  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1172  }
1173  }
1174  else
1175  {
1176  for (j = 0; j < e_ncoeffs; j++)
1177  {
1178  outarray[offset + j] +=
1179  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1180  }
1181  }
1182 #endif
1183  }
1184  }
1185  }
1186 
1187 
1189  const Array<OneD, const NekDouble> &inarray,
1190  Array<OneD, NekDouble> &outarray,
1191  const FlagList &flags,
1192  const StdRegions::ConstFactorMap &factors,
1193  const StdRegions::VarCoeffMap &varcoeff,
1194  const Array<OneD, const NekDouble> &dirForcing)
1195  {
1196  int i,n,cnt,nbndry;
1197  int nexp = GetExpSize();
1199  DNekVec F(m_ncoeffs,f,eWrapper);
1200  Array<OneD,NekDouble> e_f, e_l;
1201 
1202  //----------------------------------
1203  // Setup RHS Inner product
1204  //----------------------------------
1205  IProductWRTBase(inarray,f);
1206  Vmath::Neg(m_ncoeffs,f,1);
1207 
1208  //----------------------------------
1209  // Solve continuous Boundary System
1210  //----------------------------------
1211  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
1212  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
1213  int e_ncoeffs,id;
1214 
1215  GlobalMatrixKey HDGLamToUKey(
1218  factors,
1219  varcoeff);
1220 
1221  const DNekScalBlkMatSharedPtr &HDGLamToU =
1222  GetBlockMatrix(HDGLamToUKey);
1223 
1224  // Retrieve global trace space storage, \Lambda, from trace expansion
1226  (m_traceMap->GetNumLocalBndCoeffs());
1227 
1228 
1229  Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
1230  // Zero trace space
1231  Vmath::Zero(GloBndDofs,BndSol,1);
1232 
1233  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
1234  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
1235  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
1236 
1237  //----------------------------------
1238  // Evaluate Trace Forcing
1239  //----------------------------------
1240  // Determing <u_lam,f> terms using HDGLamToU matrix
1241  for (cnt = n = 0; n < nexp; ++n)
1242  {
1243  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1244 
1245  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1246  e_f = f+m_coeff_offset[n];
1247  e_l = loc_lambda + cnt;
1248 
1249  // use outarray as tmp space
1250  DNekVec Floc (nbndry, e_l, eWrapper);
1251  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
1252  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1253 
1254  cnt += nbndry;
1255  }
1256 
1257  // Assemble into global operator
1258  m_traceMap->AssembleBnd(loc_lambda,BndRhs);
1259 
1260  cnt = 0;
1261  // Copy Dirichlet boundary conditions into trace space
1262  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1263  {
1264  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1266  {
1267  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1268  BndSol[id] = m_bndCondExpansions[i]->GetCoeff(0);
1269  }
1270  else
1271  {
1272  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1273  BndRhs[id] += m_bndCondExpansions[i]->GetCoeff(0);
1274  }
1275  }
1276 
1277  //----------------------------------
1278  // Solve trace problem
1279  //----------------------------------
1280  if (GloBndDofs - NumDirBCs > 0)
1281  {
1283  m_traceMap,factors);
1285  LinSys->Solve(BndRhs,BndSol,m_traceMap);
1286  }
1287 
1288  //----------------------------------
1289  // Internal element solves
1290  //----------------------------------
1293  factors);
1294 
1295  const DNekScalBlkMatSharedPtr& InvHDGHelm =
1296  GetBlockMatrix(invHDGhelmkey);
1297  DNekVec out(m_ncoeffs,outarray,eWrapper);
1298  Vmath::Zero(m_ncoeffs,outarray,1);
1299 
1300  // get local trace solution from BndSol
1301  m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1302 
1303  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
1304  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
1305  }
1306 
1307  /**
1308  * Based on the expression \f$g(x,t)\f$ for the boundary conditions,
1309  * this function evaluates the boundary conditions for all boundaries
1310  * at time-level \a t.
1311  * @param time The time at which the boundary conditions
1312  * should be evaluated.
1313  * @param bndCondExpansions List of boundary expansions.
1314  * @param bndConditions Information about the boundary conditions.
1315  */
1317  const NekDouble time,
1318  const std::string varName,
1319  const NekDouble x2_in,
1320  const NekDouble x3_in)
1321  {
1322  int i;
1323 
1324  Array<OneD, NekDouble> x0(1);
1325  Array<OneD, NekDouble> x1(1);
1326  Array<OneD, NekDouble> x2(1);
1327 
1328  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1329  {
1330  if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
1331  {
1332  m_bndCondExpansions[i]->GetCoords(x0, x1, x2);
1333 
1334  if (x2_in != NekConstants::kNekUnsetDouble && x3_in !=
1336  {
1337  x1[0] = x2_in;
1338  x2[0] = x3_in;
1339  }
1340 
1341  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1343  {
1344  m_bndCondExpansions[i]->SetCoeff(0,
1345  (boost::static_pointer_cast<SpatialDomains
1346  ::DirichletBoundaryCondition>(m_bndConditions[i])
1347  ->m_dirichletCondition).Evaluate(x0[0],x1[0],x2[0],time));
1348  m_bndCondExpansions[i]->SetPhys(0,m_bndCondExpansions[i]->GetCoeff(0));
1349  }
1350  else if (m_bndConditions[i]->GetBoundaryConditionType()
1352  {
1353  m_bndCondExpansions[i]->SetCoeff(0,
1354  (boost::static_pointer_cast<SpatialDomains
1355  ::NeumannBoundaryCondition>(m_bndConditions[i])
1356  ->m_neumannCondition).Evaluate(x0[0],x1[0],x2[0],time));
1357  }
1358  else if (m_bndConditions[i]->GetBoundaryConditionType()
1360  {
1361  m_bndCondExpansions[i]->SetCoeff(0,
1362  (boost::static_pointer_cast<SpatialDomains
1363  ::RobinBoundaryCondition>(m_bndConditions[i])
1364  ->m_robinFunction).Evaluate(x0[0],x1[0],x2[0],time));
1365 
1366  m_bndCondExpansions[i]->SetPhys(0,
1367  (boost::static_pointer_cast<SpatialDomains
1368  ::RobinBoundaryCondition>(m_bndConditions[i])
1369  ->m_robinPrimitiveCoeff).Evaluate(x0[0],x1[0],x2[0],time));
1370 
1371  }
1372  else if (m_bndConditions[i]->GetBoundaryConditionType()
1374  {
1375  }
1376  else
1377  {
1378  ASSERTL0(false, "This type of BC not implemented yet");
1379  }
1380  }
1381  }
1382  }
1383 
1384  // Set up a list of element ids and edge ids that link to the
1385  // boundary conditions
1387  Array<OneD, int> &ElmtID, Array<OneD,int> &VertID)
1388  {
1389  map<int, int> VertGID;
1390  int i,n,id;
1391  int bid,cnt,Vid;
1392  int nbcs = m_bndConditions.num_elements();
1393 
1394  // make sure arrays are of sufficient length
1395  if (ElmtID.num_elements() != nbcs)
1396  {
1397  ElmtID = Array<OneD, int>(nbcs,-1);
1398  }
1399  else
1400  {
1401  fill(ElmtID.get(), ElmtID.get()+nbcs, -1);
1402  }
1403 
1404  if (VertID.num_elements() != nbcs)
1405  {
1406  VertID = Array<OneD, int>(nbcs);
1407  }
1408 
1409  // setup map of all global ids along boundary
1410  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1411  {
1412  Vid = m_bndCondExpansions[n]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
1413  VertGID[Vid] = cnt++;
1414  }
1415 
1416  // Loop over elements and find verts that match;
1418  for(cnt = n = 0; n < GetExpSize(); ++n)
1419  {
1420  exp = (*m_exp)[n];
1421  for(i = 0; i < exp->GetNverts(); ++i)
1422  {
1423  id = exp->GetGeom()->GetVid(i);
1424 
1425  if (VertGID.count(id) > 0)
1426  {
1427  bid = VertGID.find(id)->second;
1428  ASSERTL1(ElmtID[bid] == -1,"Edge already set");
1429  ElmtID[bid] = n;
1430  VertID[bid] = i;
1431  cnt ++;
1432  }
1433  }
1434  }
1435 
1436  ASSERTL1(cnt == nbcs,"Failed to visit all boundary condtiions");
1437  }
1438 
1439  /**
1440  * @brief Reset this field, so that geometry information can be updated.
1441  */
1443  {
1444  ExpList::v_Reset();
1445 
1446  // Reset boundary condition expansions.
1447  for (int n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1448  {
1449  m_bndCondExpansions[n]->Reset();
1450  }
1451  }
1452 
1453  /**
1454  * Search through the edge expansions and identify which ones
1455  * have Robin/Mixed type boundary conditions. If find a Robin
1456  * boundary then store the edge id of the boundary condition
1457  * and the array of points of the physical space boundary
1458  * condition which are hold the boundary condition primitive
1459  * variable coefficient at the quatrature points
1460  *
1461  * \return std map containing the robin boundary condition
1462  * info using a key of the element id
1463  *
1464  * There is a next member to allow for more than one Robin
1465  * boundary condition per element
1466  */
1467  map<int, RobinBCInfoSharedPtr> DisContField1D::v_GetRobinBCInfo(void)
1468  {
1469  int i;
1470  map<int, RobinBCInfoSharedPtr> returnval;
1471  Array<OneD, int> ElmtID,VertID;
1472  GetBoundaryToElmtMap(ElmtID,VertID);
1473 
1474  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1475  {
1476  MultiRegions::ExpListSharedPtr locExpList;
1477 
1478  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1480  {
1481  int elmtid;
1482  Array<OneD, NekDouble> Array_tmp;
1483 
1484  locExpList = m_bndCondExpansions[i];
1485 
1486  RobinBCInfoSharedPtr rInfo =
1489  VertID[i],Array_tmp = locExpList->GetPhys());
1490 
1491  elmtid = ElmtID[i];
1492  // make link list if necessary (not likely in
1493  // 1D but needed in 2D & 3D)
1494  if(returnval.count(elmtid) != 0)
1495  {
1496  rInfo->next = returnval.find(elmtid)->second;
1497  }
1498  returnval[elmtid] = rInfo;
1499  }
1500  }
1501 
1502  return returnval;
1503  }
1504  } // end of namespace
1505 } //end of namespace
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:772
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1394
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n...
Definition: ExpList.h:1867
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:1343
void EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble=NekConstants::kNekUnsetDouble, const NekDouble=NekConstants::kNekUnsetDouble)
Definition: ExpList.h:2027
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
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:1875
const BoundaryConditionCollection & GetBoundaryConditions(void) const
Definition: Conditions.h:235
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 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.
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:54
SpatialDomains::BoundaryConditionsSharedPtr GetDomainBCs(const SpatialDomains::CompositeMap &domain, const SpatialDomains::BoundaryConditions &Allbcs, const std::string &variable)
vector< int > m_periodicFwdCopy
A vector indicating degress of freedom which need to be copied from forwards to backwards space in ca...
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
void SetUpDG(const std::string &variable)
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2080
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:926
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:1381
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:1837
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
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:1152
void SetAdjacentElementExp(int vertex, Expansion1DSharedPtr &v)
Definition: Expansion0D.h:106
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:47
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:958
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:53
The base class for all shapes.
Definition: StdExpansion.h:69
boost::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:1942
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:947
DisContField1D()
Default constructor.
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:225
boost::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: ExpList.h:1340
int id
Geometry ID of entity.
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:935
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:887
std::map< int, Composite >::const_iterator CompositeMapConstIter
Definition: MeshGraph.h:114
std::map< int, BoundaryRegionShPtr > BoundaryRegionCollection
Definition: Conditions.h:204
boost::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:74
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:883
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:880
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:382
double NekDouble
Defines a list of flags.
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:213
boost::shared_ptr< GeometryVector > Composite
Definition: MeshGraph.h:111
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)
Solve the Helmholtz equation.
static const NekDouble kNekUnsetDouble
Describe a linear system.
std::map< int, Composite > CompositeMap
Definition: MeshGraph.h:112
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Describes a matrix with ordering defined by a local to global map.
const Array< OneD, const NekDouble > & GetPhys() const
This function returns (a reference to) the array (implemented as m_phys) containing the function ev...
Definition: ExpList.h:1828
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:212
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:213
NekDouble GetCoeff(int i)
Get the i th value (coefficient) of m_coeffs.
Definition: ExpList.h:1806
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:269
boost::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1507
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.
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:2631
boost::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:51
const BoundaryRegionCollection & GetBoundaryRegions(void) const
Definition: Conditions.h:225
boost::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:206
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:714
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:359
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:191
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:432
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:53
boost::shared_ptr< BoundaryRegion > BoundaryRegionShPtr
Definition: Conditions.h:202
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)