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 
37 #include <MultiRegions/DisContField1D.h>
38 #include <StdRegions/StdSegExp.h>
39 #include <LibUtilities/Foundations/ManagerAccess.h>
40 #include <SpatialDomains/MeshGraph1D.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  */
56  DisContField1D::DisContField1D():
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  */
76  DisContField1D::DisContField1D(
77  const LibUtilities::SessionReaderSharedPtr &pSession,
78  const SpatialDomains::MeshGraphSharedPtr &graph1D,
79  const std::string &variable,
80  const bool SetUpJustDG)
81  : ExpList1D(pSession,graph1D),
82  m_bndCondExpansions(),
83  m_bndConditions()
84  {
85  SpatialDomains::BoundaryConditions bcs(m_session, graph1D);
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
102  if (m_trace != NullExpListSharedPtr)
103  {
104  return;
105  }
106 
107  SpatialDomains::MeshGraph1DSharedPtr graph1D =
108  boost::dynamic_pointer_cast<SpatialDomains::MeshGraph1D>(
109  m_graph);
110 
111  m_globalBndMat = MemoryManager<GlobalLinSysMap>::AllocateSharedPtr();
112 
113  ExpList0DSharedPtr trace = MemoryManager<ExpList0D>::
114  AllocateSharedPtr(
115  m_bndCondExpansions,
116  m_bndConditions,
117  *m_exp,graph1D,
118  m_periodicVerts);
119 
120  m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
121 
122  m_traceMap = MemoryManager<AssemblyMapDG>::
123  AllocateSharedPtr(m_session, graph1D, trace, *this,
124  m_bndCondExpansions, m_bndConditions,
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;
135  LocalRegions::Expansion0DSharedPtr exp0d;
136  LocalRegions::Expansion1DSharedPtr exp1d;
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 
158  SetUpPhysNormals();
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() !=
166  SpatialDomains::ePeriodic)
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  {
189  m_leftAdjacentVerts[cnt] = IsLeftAdjacentVertex(n, v);
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  {
200  PeriodicMap::iterator it = m_periodicVerts.find(
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  {
255  set<int>::iterator it;
256  LocalRegions::Expansion0DSharedPtr traceEl =
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();
273  PeriodicMap::iterator pIt = m_periodicVerts.find(
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.
308  SpatialDomains::BoundaryConditionsSharedPtr DisContField1D::GetDomainBCs(
309  const SpatialDomains::CompositeMap &domain,
310  const SpatialDomains::BoundaryConditions &Allbcs,
311  const std::string &variable)
312  {
313  SpatialDomains::BoundaryConditionsSharedPtr returnval;
314 
315  returnval = MemoryManager<SpatialDomains::BoundaryConditions>::AllocateSharedPtr();
316 
317  SpatialDomains::BoundaryRegion::iterator bregionIt;
318  map<int,int> GeometryToRegionsMap;
319 
320  SpatialDomains::BoundaryRegionCollection::const_iterator it;
321 
322  const SpatialDomains::BoundaryRegionCollection &bregions
323  = Allbcs.GetBoundaryRegions();
324  const SpatialDomains::BoundaryConditionCollection &bconditions
325  = Allbcs.GetBoundaryConditions();
326 
327  // Set up a map of all boundary regions
328  for(it = bregions.begin(); it != bregions.end(); ++it)
329  {
330  SpatialDomains::BoundaryRegion::iterator bregionIt;
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 
341  SpatialDomains::CompositeMapConstIter domIt;
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 
366  map<int,SpatialDomains::GeometrySharedPtr>::iterator regIt;
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  {
387  SpatialDomains::BoundaryRegionShPtr breg(MemoryManager<SpatialDomains::BoundaryRegion>::AllocateSharedPtr());
388 
389  // Set up Composite (GemetryVector) to contain vertex and put into bRegion
390  SpatialDomains::Composite gvec(MemoryManager<SpatialDomains::GeometryVector>::AllocateSharedPtr());
391  gvec->push_back(regIt->second);
392  (*breg)[regIt->first] = gvec;
393 
394  returnval->AddBoundaryRegions(bregions.size()+numNewBc,breg);
395 
396  SpatialDomains::BoundaryConditionMapShPtr bCondition = MemoryManager<SpatialDomains::BoundaryConditionMap>::AllocateSharedPtr();
397 
398  // Set up just boundary condition for this variable.
399  SpatialDomains::BoundaryConditionShPtr notDefinedCondition(MemoryManager<SpatialDomains::NotDefinedBoundaryCondition>::AllocateSharedPtr(m_session, "0"));
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  */
418  DisContField1D::DisContField1D(
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  */
443  DisContField1D::DisContField1D(const DisContField1D &In):
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  */
463  DisContField1D::DisContField1D(const ExpList1D &In):
464  ExpList1D(In)
465  {
466  }
467 
468  /**
469  *
470  */
471  DisContField1D::~DisContField1D()
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  */
485  void DisContField1D::GenerateBoundaryConditionExpansion(
486  const SpatialDomains::MeshGraphSharedPtr &graph1D,
487  const SpatialDomains::BoundaryConditions &bcs,
488  const std::string variable)
489  {
490  int cnt = 0;
491 
492  const SpatialDomains::BoundaryRegionCollection &bregions
493  = bcs.GetBoundaryRegions();
494  const SpatialDomains::BoundaryConditionCollection &bconditions
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() !=
504  SpatialDomains::ePeriodic )
505  {
506  SpatialDomains::BoundaryRegion::iterator bregionIt;
507  for (bregionIt = it->second->begin();
508  bregionIt != it->second->end(); bregionIt++)
509  {
510  cnt += bregionIt->second->size();
511  }
512  }
513  }
514 
515  m_bndCondExpansions
516  = Array<OneD,MultiRegions::ExpListSharedPtr>(cnt);
517 
518  m_bndConditions
519  = Array<OneD,SpatialDomains::BoundaryConditionShPtr>(cnt);
520 
521  SetBoundaryConditionExpansion(graph1D,bcs,variable,
522  m_bndCondExpansions,
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  */
535  void DisContField1D::FindPeriodicVertices(
536  const SpatialDomains::BoundaryConditions &bcs,
537  const std::string variable)
538  {
539  const SpatialDomains::BoundaryRegionCollection &bregions
540  = bcs.GetBoundaryRegions();
541  const SpatialDomains::BoundaryConditionCollection &bconditions
542  = bcs.GetBoundaryConditions();
543 
544  SpatialDomains::MeshGraph1DSharedPtr graph1D
545  = boost::dynamic_pointer_cast<
546  SpatialDomains::MeshGraph1D>(m_graph);
547  SpatialDomains::BoundaryRegionCollection::const_iterator it;
548 
549  LibUtilities::CommSharedPtr vComm =
550  m_session->GetComm()->GetRowComm();
551 
552  int i, region1ID, region2ID;
553 
554  SpatialDomains::BoundaryConditionShPtr locBCond;
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()
565  != SpatialDomains::ePeriodic)
566  {
567  continue;
568  }
569 
570  int id = (*(it->second->begin()->second))[0]->GetGlobalID();
571 
572  BregionToVertMap[it->first] = id;
573  }
574 
575  map<int,int>::iterator iit;
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()
618  != SpatialDomains::ePeriodic)
619  {
620  continue;
621  }
622 
623  // Identify periodic boundary region IDs.
624  region1ID = it->first;
625  region2ID = boost::static_pointer_cast<
626  SpatialDomains::PeriodicBoundaryCondition>(
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],
636  StdRegions::eNoOrientation,
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  */
656  void DisContField1D::SetBoundaryConditionExpansion(
657  const SpatialDomains::MeshGraphSharedPtr &graph1D,
658  const SpatialDomains::BoundaryConditions &bcs,
659  const std::string variable,
660  Array<OneD, MultiRegions::ExpListSharedPtr>
661  &bndCondExpansions,
662  Array<OneD, SpatialDomains
663  ::BoundaryConditionShPtr> &bndConditions)
664  {
665  int k;
666  int cnt = 0;
667 
668  const SpatialDomains::BoundaryRegionCollection &bregions
669  = bcs.GetBoundaryRegions();
670  const SpatialDomains::BoundaryConditionCollection &bconditions
671  = bcs.GetBoundaryConditions();
672  SpatialDomains::BoundaryRegionCollection::const_iterator it;
673 
674  MultiRegions::ExpList0DSharedPtr locPointExp;
675  SpatialDomains::BoundaryConditionShPtr locBCond;
676  SpatialDomains::PointGeomSharedPtr vert;
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() ==
686  SpatialDomains::eDirichlet)
687 
688  {
689  SpatialDomains::BoundaryRegion::iterator bregionIt;
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
696  <SpatialDomains::PointGeom>(
697  (*bregionIt->second)[k])))
698  {
699  locPointExp
700  = MemoryManager<MultiRegions::ExpList0D>
701  ::AllocateSharedPtr(vert);
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  {
722  case SpatialDomains::eNeumann:
723  case SpatialDomains::eRobin:
724  case SpatialDomains::eNotDefined: // presume this will be reused as Neuman, Robin or Dirichlet later
725  {
726  SpatialDomains::BoundaryRegion::iterator bregionIt;
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
733  <SpatialDomains::PointGeom>(
734  (*bregionIt->second)[k])))
735  {
736  locPointExp
737  = MemoryManager<MultiRegions::ExpList0D>
738  ::AllocateSharedPtr(vert);
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
751  case SpatialDomains::eDirichlet:
752  case SpatialDomains::ePeriodic:
753  break;
754  default:
755  ASSERTL0(false,"This type of BC not implemented yet");
756  break;
757  }
758  }
759  }
760 
761  /**
762  *
763  */
764  GlobalLinSysSharedPtr DisContField1D::GetGlobalBndLinSys(
765  const GlobalLinSysKey &mkey)
766  {
767  ASSERTL0(mkey.GetMatrixType() == StdRegions::eHybridDGHelmBndLam,
768  "Routine currently only tested for HybridDGHelmholtz");
769 
770  ASSERTL1(mkey.GetGlobalSysSolnType() != eDirectFullMatrix,
771  "Full matrix global systems are not supported for HDG "
772  "expansions");
773 
774  ASSERTL1(mkey.GetGlobalSysSolnType()
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 
796  vector<bool> &DisContField1D::GetNegatedFluxNormal(void)
797  {
798  if(m_negatedFluxNormal.size() == 0)
799  {
800  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
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 
811  LocalRegions::Expansion0DSharedPtr vertExp =
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  */
835  void DisContField1D::v_GetFwdBwdTracePhys(Array<OneD, NekDouble> &Fwd,
836  Array<OneD, NekDouble> &Bwd)
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 
868  void DisContField1D::v_GetFwdBwdTracePhys(
869  const Array<OneD, const NekDouble> &field,
870  Array<OneD, NekDouble> &Fwd,
871  Array<OneD, NekDouble> &Bwd)
872  {
873  // Expansion casts
874  LocalRegions::Expansion1DSharedPtr exp1D;
875  LocalRegions::Expansion1DSharedPtr exp1DFirst;
876  LocalRegions::Expansion1DSharedPtr exp1DLast;
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 
887  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
888  &elmtToTrace = m_traceMap->GetElmtToTrace();
889 
890  // Basis shared pointer
891  LibUtilities::BasisSharedPtr Basis;
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() ==
932  SpatialDomains::eDirichlet)
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() ==
939  SpatialDomains::eNeumann ||
940  m_bndConditions[n]->GetBoundaryConditionType() ==
941  SpatialDomains::eRobin)
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() ==
952  SpatialDomains::eNotDefined)
953  {
954  // Do nothing
955  }
956  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
957  SpatialDomains::ePeriodic)
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 
977  void DisContField1D::v_ExtractTracePhys(
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  */
997  void DisContField1D::v_ExtractTracePhys(
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 
1022  void DisContField1D::v_AddTraceIntegral(
1023  const Array<OneD, const NekDouble> &Fn,
1024  Array<OneD, NekDouble> &outarray)
1025  {
1026  int n,offset, t_offset;
1027 
1028  Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
1029  &elmtToTrace = m_traceMap->GetElmtToTrace();
1030 
1031  // Basis shared pointer
1032  LibUtilities::BasisSharedPtr Basis;
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;
1075  LibUtilities::BasisSharedPtr BASE;
1076  const LibUtilities::PointsKey
1077  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1078  const LibUtilities::BasisKey
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  {
1128  LibUtilities::BasisSharedPtr BASE;
1129  const LibUtilities::PointsKey
1130  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1131  const LibUtilities::BasisKey
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 
1188  void DisContField1D::v_HelmSolve(
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();
1198  Array<OneD,NekDouble> f(m_ncoeffs);
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(
1216  StdRegions::eHybridDGLamToU,
1217  NullAssemblyMapSharedPtr,
1218  factors,
1219  varcoeff);
1220 
1221  const DNekScalBlkMatSharedPtr &HDGLamToU =
1222  GetBlockMatrix(HDGLamToUKey);
1223 
1224  // Retrieve global trace space storage, \Lambda, from trace expansion
1225  Array<OneD,NekDouble> BndSol = Array<OneD,NekDouble>
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() ==
1265  SpatialDomains::eDirichlet)
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  {
1282  GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam,
1283  m_traceMap,factors);
1284  GlobalLinSysSharedPtr LinSys = GetGlobalBndLinSys(key);
1285  LinSys->Solve(BndRhs,BndSol,m_traceMap);
1286  }
1287 
1288  //----------------------------------
1289  // Internal element solves
1290  //----------------------------------
1291  GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
1292  NullAssemblyMapSharedPtr,
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  */
1316  void DisContField1D::v_EvaluateBoundaryConditions(
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]->GetUserDefined() ==
1331  SpatialDomains::eTimeDependent ||
1332  m_bndConditions[i]->GetUserDefined() ==
1333  SpatialDomains::eQinflow ||
1334  m_bndConditions[i]->GetUserDefined() ==
1335  SpatialDomains::eRCRterminal )
1336  {
1337  m_bndCondExpansions[i]->GetCoords(x0, x1, x2);
1338 
1339  if (x2_in != NekConstants::kNekUnsetDouble && x3_in !=
1340  NekConstants::kNekUnsetDouble)
1341  {
1342  x1[0] = x2_in;
1343  x2[0] = x3_in;
1344  }
1345 
1346  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1347  SpatialDomains::eDirichlet)
1348  {
1349  m_bndCondExpansions[i]->SetCoeff(0,
1350  (boost::static_pointer_cast<SpatialDomains
1351  ::DirichletBoundaryCondition>(m_bndConditions[i])
1352  ->m_dirichletCondition).Evaluate(x0[0],x1[0],x2[0],time));
1353  m_bndCondExpansions[i]->SetPhys(0,m_bndCondExpansions[i]->GetCoeff(0));
1354  }
1355  else if (m_bndConditions[i]->GetBoundaryConditionType()
1356  == SpatialDomains::eNeumann)
1357  {
1358  m_bndCondExpansions[i]->SetCoeff(0,
1359  (boost::static_pointer_cast<SpatialDomains
1360  ::NeumannBoundaryCondition>(m_bndConditions[i])
1361  ->m_neumannCondition).Evaluate(x0[0],x1[0],x2[0],time));
1362  }
1363  else if (m_bndConditions[i]->GetBoundaryConditionType()
1364  == SpatialDomains::eRobin)
1365  {
1366  m_bndCondExpansions[i]->SetCoeff(0,
1367  (boost::static_pointer_cast<SpatialDomains
1368  ::RobinBoundaryCondition>(m_bndConditions[i])
1369  ->m_robinFunction).Evaluate(x0[0],x1[0],x2[0],time));
1370 
1371  m_bndCondExpansions[i]->SetPhys(0,
1372  (boost::static_pointer_cast<SpatialDomains
1373  ::RobinBoundaryCondition>(m_bndConditions[i])
1374  ->m_robinPrimitiveCoeff).Evaluate(x0[0],x1[0],x2[0],time));
1375 
1376  }
1377  else if (m_bndConditions[i]->GetBoundaryConditionType()
1378  == SpatialDomains::eNotDefined)
1379  {
1380  }
1381  else
1382  {
1383  ASSERTL0(false, "This type of BC not implemented yet");
1384  }
1385  }
1386  }
1387  }
1388 
1389  // Set up a list of element ids and edge ids that link to the
1390  // boundary conditions
1391  void DisContField1D::v_GetBoundaryToElmtMap(
1392  Array<OneD, int> &ElmtID, Array<OneD,int> &VertID)
1393  {
1394  map<int, int> VertGID;
1395  int i,n,id;
1396  int bid,cnt,Vid;
1397  int nbcs = m_bndConditions.num_elements();
1398 
1399  // make sure arrays are of sufficient length
1400  if (ElmtID.num_elements() != nbcs)
1401  {
1402  ElmtID = Array<OneD, int>(nbcs,-1);
1403  }
1404  else
1405  {
1406  fill(ElmtID.get(), ElmtID.get()+nbcs, -1);
1407  }
1408 
1409  if (VertID.num_elements() != nbcs)
1410  {
1411  VertID = Array<OneD, int>(nbcs);
1412  }
1413 
1414  // setup map of all global ids along boundary
1415  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1416  {
1417  Vid = m_bndCondExpansions[n]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
1418  VertGID[Vid] = cnt++;
1419  }
1420 
1421  // Loop over elements and find verts that match;
1422  LocalRegions::ExpansionSharedPtr exp;
1423  for(cnt = n = 0; n < GetExpSize(); ++n)
1424  {
1425  exp = (*m_exp)[n];
1426  for(i = 0; i < exp->GetNverts(); ++i)
1427  {
1428  id = exp->GetGeom()->GetVid(i);
1429 
1430  if (VertGID.count(id) > 0)
1431  {
1432  bid = VertGID.find(id)->second;
1433  ASSERTL1(ElmtID[bid] == -1,"Edge already set");
1434  ElmtID[bid] = n;
1435  VertID[bid] = i;
1436  cnt ++;
1437  }
1438  }
1439  }
1440 
1441  ASSERTL1(cnt == nbcs,"Failed to visit all boundary condtiions");
1442  }
1443 
1444  /**
1445  * Search through the edge expansions and identify which ones
1446  * have Robin/Mixed type boundary conditions. If find a Robin
1447  * boundary then store the edge id of the boundary condition
1448  * and the array of points of the physical space boundary
1449  * condition which are hold the boundary condition primitive
1450  * variable coefficient at the quatrature points
1451  *
1452  * \return std map containing the robin boundary condition
1453  * info using a key of the element id
1454  *
1455  * There is a next member to allow for more than one Robin
1456  * boundary condition per element
1457  */
1458  map<int, RobinBCInfoSharedPtr> DisContField1D::v_GetRobinBCInfo(void)
1459  {
1460  int i;
1461  map<int, RobinBCInfoSharedPtr> returnval;
1462  Array<OneD, int> ElmtID,VertID;
1463  GetBoundaryToElmtMap(ElmtID,VertID);
1464 
1465  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1466  {
1467  MultiRegions::ExpListSharedPtr locExpList;
1468 
1469  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1470  SpatialDomains::eRobin)
1471  {
1472  int elmtid;
1473  Array<OneD, NekDouble> Array_tmp;
1474 
1475  locExpList = m_bndCondExpansions[i];
1476 
1477  RobinBCInfoSharedPtr rInfo =
1478  MemoryManager<RobinBCInfo>::
1479  AllocateSharedPtr(
1480  VertID[i],Array_tmp = locExpList->GetPhys());
1481 
1482  elmtid = ElmtID[i];
1483  // make link list if necessary (not likely in
1484  // 1D but needed in 2D & 3D)
1485  if(returnval.count(elmtid) != 0)
1486  {
1487  rInfo->next = returnval.find(elmtid)->second;
1488  }
1489  returnval[elmtid] = rInfo;
1490  }
1491  }
1492 
1493  return returnval;
1494  }
1495  } // end of namespace
1496 } //end of namespace
Generated on Wed Mar 4 2015 16:33:13 for Nektar++ by   doxygen 1.8.1.2