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