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, StdRegions::StdExpansionSharedPtr> >
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  // Number of solution points of each element
885  int nLocalSolutionPts;
886 
887  // Initial index of each element
888  int phys_offset;
889 
890  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
891  &elmtToTrace = m_traceMap->GetElmtToTrace();
892 
893  // Basis shared pointer
894  LibUtilities::BasisSharedPtr Basis;
895 
896  // Set forward and backard state to zero
897  Vmath::Zero(Fwd.num_elements(), Fwd, 1);
898  Vmath::Zero(Bwd.num_elements(), Bwd, 1);
899 
900  int cnt;
901 
902  // Loop on the elements
903  for (cnt = n = 0; n < nElements; ++n)
904  {
905  exp1D = (*m_exp)[n]->as<LocalRegions::Expansion1D>();
906 
907  // Set the offset of each element
908  phys_offset = GetPhys_Offset(n);
909 
910  // Set the number of solution points of each element
911  nLocalSolutionPts = (*m_exp)[n]->GetNumPoints(0);
912 
913  Basis = (*m_exp)[n]->GetBasis(0);
914 
915  for(v = 0; v < 2; ++v, ++cnt)
916  {
917  int offset = m_trace->GetPhys_Offset(elmtToTrace[n][v]->GetElmtId());
918 
919  if (m_leftAdjacentVerts[cnt])
920  {
921  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
922  Fwd[offset]);
923  }
924  else
925  {
926  (*m_exp)[n]->GetVertexPhysVals(v, field + phys_offset,
927  Bwd[offset]);
928  }
929  }
930  }
931 
932  // Fill boundary conditions into missing elements.
933  int id = 0;
934 
935  for(cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
936  {
937  if (m_bndConditions[n]->GetBoundaryConditionType() ==
938  SpatialDomains::eDirichlet)
939  {
940  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
941  Bwd[id] = m_bndCondExpansions[n]->GetPhys()[0]; //this is not getting the correct value?
942  cnt++;
943  }
944  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
945  SpatialDomains::eNeumann ||
946  m_bndConditions[n]->GetBoundaryConditionType() ==
947  SpatialDomains::eRobin)
948  {
949  ASSERTL0((m_bndCondExpansions[n]->GetPhys())[0]==0.0,
950  "Method not set up for non-zero Neumann "
951  "boundary condition");
952  id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
953  Bwd[id] = Fwd[id];
954 
955  cnt++;
956  }
957  else if (m_bndConditions[n]->GetBoundaryConditionType() ==
958  SpatialDomains::eNotDefined)
959  {
960  // Do nothing
961  }
962  else if (m_bndConditions[n]->GetBoundaryConditionType() !=
963  SpatialDomains::ePeriodic)
964  {
965  ASSERTL0(false,
966  "Method not set up for this boundary condition.");
967  }
968  }
969 
970  // Copy any periodic boundary conditions.
971  for (n = 0; n < m_periodicFwdCopy.size(); ++n)
972  {
973  Bwd[m_periodicBwdCopy[n]] = Fwd[m_periodicFwdCopy[n]];
974  }
975 
976  // Do parallel exchange for forwards/backwards spaces.
977  m_traceMap->UniversalTraceAssemble(Fwd);
978  m_traceMap->UniversalTraceAssemble(Bwd);
979 
980  }
981 
982 
983  void DisContField1D::v_ExtractTracePhys(
984  Array<OneD, NekDouble> &outarray)
985  {
986  ASSERTL1(m_physState == true,"local physical space is not true ");
987  v_ExtractTracePhys(m_phys, outarray);
988  }
989 
990  /**
991  * @brief This method extracts the trace (verts in 1D) from the field @a
992  * inarray and puts the values in @a outarray.
993  *
994  * It assumes the field is C0 continuous so that it can overwrite the
995  * edge data when visited by the two adjacent elements.
996  *
997  * @param inarray An array containing the 1D data from which we wish
998  * to extract the edge data.
999  * @param outarray The resulting edge information.
1000  *
1001  * This will not work for non-boundary expansions
1002  */
1003  void DisContField1D::v_ExtractTracePhys(
1004  const Array<OneD, const NekDouble> &inarray,
1005  Array<OneD, NekDouble> &outarray)
1006  {
1007  // Loop over elemente and collect forward expansion
1008  int nexp = GetExpSize();
1009  int n,p,offset,phys_offset;
1010 
1011  ASSERTL1(outarray.num_elements() >= m_trace->GetExpSize(),
1012  "input array is of insufficient length");
1013 
1014  for (n = 0; n < nexp; ++n)
1015  {
1016  phys_offset = GetPhys_Offset(n);
1017 
1018  for (p = 0; p < (*m_exp)[n]->GetNverts(); ++p)
1019  {
1020  offset = m_trace->GetPhys_Offset(
1021  (m_traceMap->GetElmtToTrace())[n][p]->GetElmtId());
1022  (*m_exp)[n]->GetVertexPhysVals(p, inarray + phys_offset,
1023  outarray[offset]);
1024  }
1025  }
1026  }
1027 
1028  void DisContField1D::v_AddTraceIntegral(
1029  const Array<OneD, const NekDouble> &Fn,
1030  Array<OneD, NekDouble> &outarray)
1031  {
1032  int n,offset, t_offset;
1033 
1034  Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
1035  &elmtToTrace = m_traceMap->GetElmtToTrace();
1036 
1037  // Basis shared pointer
1038  LibUtilities::BasisSharedPtr Basis;
1039 
1040  vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
1041 
1042  for (n = 0; n < GetExpSize(); ++n)
1043  {
1044  // Basis definition on each element
1045  Basis = (*m_exp)[n]->GetBasis(0);
1046 
1047  // Number of coefficients on each element
1048  int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1049 
1050  offset = GetCoeff_Offset(n);
1051 
1052  // Implementation for every points except Gauss points
1053  if (Basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
1054  {
1055  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1056  if(negatedFluxNormal[2*n])
1057  {
1058  outarray[offset] -= Fn[t_offset];
1059  }
1060  else
1061  {
1062  outarray[offset] += Fn[t_offset];
1063  }
1064 
1065  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1066 
1067  if(negatedFluxNormal[2*n+1])
1068  {
1069  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] -= Fn[t_offset];
1070  }
1071  else
1072  {
1073  outarray[offset+(*m_exp)[n]->GetVertexMap(1)] += Fn[t_offset];
1074  }
1075 
1076  }
1077  else
1078  {
1079 #if 0
1080  DNekMatSharedPtr m_Ixm;
1081  LibUtilities::BasisSharedPtr BASE;
1082  const LibUtilities::PointsKey
1083  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1084  const LibUtilities::BasisKey
1085  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1086 
1087  BASE = LibUtilities::BasisManager()[BS_k];
1088 
1089  Array<OneD, NekDouble> coords(3, 0.0);
1090 
1091  int j;
1092 
1093  for(p = 0; p < 2; ++p)
1094  {
1095  NekDouble vertnorm = 0.0;
1096  for (int i=0; i<((*m_exp)[n]->
1097  GetVertexNormal(p)).num_elements(); i++)
1098  {
1099  vertnorm += ((*m_exp)[n]->GetVertexNormal(p))[i][0];
1100  coords[0] = vertnorm ;
1101  }
1102 
1103  t_offset = GetTrace()->GetPhys_Offset(n+p);
1104 
1105  if (vertnorm >= 0.0)
1106  {
1107  m_Ixm = BASE->GetI(coords);
1108 
1109 
1110  for (j = 0; j < e_ncoeffs; j++)
1111  {
1112  outarray[offset + j] +=
1113  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1114  }
1115  }
1116 
1117  if (vertnorm < 0.0)
1118  {
1119  m_Ixm = BASE->GetI(coords);
1120 
1121  for (j = 0; j < e_ncoeffs; j++)
1122  {
1123  outarray[offset + j] -=
1124  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1125  }
1126  }
1127  }
1128 #else
1129  int j;
1130  static DNekMatSharedPtr m_Ixm, m_Ixp;
1131  static int sav_ncoeffs = 0;
1132  if(!m_Ixm.get() || e_ncoeffs != sav_ncoeffs)
1133  {
1134  LibUtilities::BasisSharedPtr BASE;
1135  const LibUtilities::PointsKey
1136  BS_p(e_ncoeffs,LibUtilities::eGaussGaussLegendre);
1137  const LibUtilities::BasisKey
1138  BS_k(LibUtilities::eGauss_Lagrange,e_ncoeffs,BS_p);
1139 
1140  BASE = LibUtilities::BasisManager()[BS_k];
1141 
1142  Array<OneD, NekDouble> coords(1, 0.0);
1143 
1144  coords[0] = -1.0;
1145  m_Ixm = BASE->GetI(coords);
1146 
1147  coords[0] = 1.0;
1148  m_Ixp = BASE->GetI(coords);
1149 
1150  sav_ncoeffs = e_ncoeffs;
1151  }
1152 
1153  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
1154  if(negatedFluxNormal[2*n])
1155  {
1156  for (j = 0; j < e_ncoeffs; j++)
1157  {
1158  outarray[offset + j] -=
1159  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1160  }
1161  }
1162  else
1163  {
1164  for (j = 0; j < e_ncoeffs; j++)
1165  {
1166  outarray[offset + j] +=
1167  (m_Ixm->GetPtr())[j] * Fn[t_offset];
1168  }
1169  }
1170 
1171  t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][1]->GetElmtId());
1172  if (negatedFluxNormal[2*n+1])
1173  {
1174  for (j = 0; j < e_ncoeffs; j++)
1175  {
1176  outarray[offset + j] -=
1177  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1178  }
1179  }
1180  else
1181  {
1182  for (j = 0; j < e_ncoeffs; j++)
1183  {
1184  outarray[offset + j] +=
1185  (m_Ixp->GetPtr())[j] * Fn[t_offset];
1186  }
1187  }
1188 #endif
1189  }
1190  }
1191  }
1192 
1193 
1194  void DisContField1D::v_HelmSolve(
1195  const Array<OneD, const NekDouble> &inarray,
1196  Array<OneD, NekDouble> &outarray,
1197  const FlagList &flags,
1198  const StdRegions::ConstFactorMap &factors,
1199  const StdRegions::VarCoeffMap &varcoeff,
1200  const Array<OneD, const NekDouble> &dirForcing)
1201  {
1202  int i,n,cnt,nbndry;
1203  int nexp = GetExpSize();
1204  Array<OneD,NekDouble> f(m_ncoeffs);
1205  DNekVec F(m_ncoeffs,f,eWrapper);
1206  Array<OneD,NekDouble> e_f, e_l;
1207 
1208  //----------------------------------
1209  // Setup RHS Inner product
1210  //----------------------------------
1211  IProductWRTBase(inarray,f);
1212  Vmath::Neg(m_ncoeffs,f,1);
1213 
1214  //----------------------------------
1215  // Solve continuous Boundary System
1216  //----------------------------------
1217  int GloBndDofs = m_traceMap->GetNumGlobalBndCoeffs();
1218  int NumDirBCs = m_traceMap->GetNumLocalDirBndCoeffs();
1219  int e_ncoeffs,id;
1220 
1221  GlobalMatrixKey HDGLamToUKey(
1222  StdRegions::eHybridDGLamToU,
1223  NullAssemblyMapSharedPtr,
1224  factors,
1225  varcoeff);
1226 
1227  const DNekScalBlkMatSharedPtr &HDGLamToU =
1228  GetBlockMatrix(HDGLamToUKey);
1229 
1230  // Retrieve global trace space storage, \Lambda, from trace expansion
1231  Array<OneD,NekDouble> BndSol = Array<OneD,NekDouble>
1232  (m_traceMap->GetNumLocalBndCoeffs());
1233 
1234 
1235  Array<OneD,NekDouble> BndRhs(GloBndDofs,0.0);
1236  // Zero trace space
1237  Vmath::Zero(GloBndDofs,BndSol,1);
1238 
1239  int LocBndCoeffs = m_traceMap->GetNumLocalBndCoeffs();
1240  Array<OneD, NekDouble> loc_lambda(LocBndCoeffs);
1241  DNekVec LocLambda(LocBndCoeffs,loc_lambda,eWrapper);
1242 
1243  //----------------------------------
1244  // Evaluate Trace Forcing
1245  //----------------------------------
1246  // Determing <u_lam,f> terms using HDGLamToU matrix
1247  for (cnt = n = 0; n < nexp; ++n)
1248  {
1249  nbndry = (*m_exp)[n]->NumDGBndryCoeffs();
1250 
1251  e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
1252  e_f = f+m_coeff_offset[n];
1253  e_l = loc_lambda + cnt;
1254 
1255  // use outarray as tmp space
1256  DNekVec Floc (nbndry, e_l, eWrapper);
1257  DNekVec ElmtFce (e_ncoeffs, e_f, eWrapper);
1258  Floc = Transpose(*(HDGLamToU->GetBlock(n,n)))*ElmtFce;
1259 
1260  cnt += nbndry;
1261  }
1262 
1263  // Assemble into global operator
1264  m_traceMap->AssembleBnd(loc_lambda,BndRhs);
1265 
1266  cnt = 0;
1267  // Copy Dirichlet boundary conditions into trace space
1268  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1269  {
1270  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1271  SpatialDomains::eDirichlet)
1272  {
1273  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1274  BndSol[id] = m_bndCondExpansions[i]->GetCoeff(0);
1275  }
1276  else
1277  {
1278  id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i);
1279  BndRhs[id] += m_bndCondExpansions[i]->GetCoeff(0);
1280  }
1281  }
1282 
1283  //----------------------------------
1284  // Solve trace problem
1285  //----------------------------------
1286  if (GloBndDofs - NumDirBCs > 0)
1287  {
1288  GlobalLinSysKey key(StdRegions::eHybridDGHelmBndLam,
1289  m_traceMap,factors);
1290  GlobalLinSysSharedPtr LinSys = GetGlobalBndLinSys(key);
1291  LinSys->Solve(BndRhs,BndSol,m_traceMap);
1292  }
1293 
1294  //----------------------------------
1295  // Internal element solves
1296  //----------------------------------
1297  GlobalMatrixKey invHDGhelmkey(StdRegions::eInvHybridDGHelmholtz,
1298  NullAssemblyMapSharedPtr,
1299  factors);
1300 
1301  const DNekScalBlkMatSharedPtr& InvHDGHelm =
1302  GetBlockMatrix(invHDGhelmkey);
1303  DNekVec out(m_ncoeffs,outarray,eWrapper);
1304  Vmath::Zero(m_ncoeffs,outarray,1);
1305 
1306  // get local trace solution from BndSol
1307  m_traceMap->GlobalToLocalBnd(BndSol,loc_lambda);
1308 
1309  // out = u_f + u_lam = (*InvHDGHelm)*f + (LamtoU)*Lam
1310  out = (*InvHDGHelm)*F + (*HDGLamToU)*LocLambda;
1311  }
1312 
1313  /**
1314  * Based on the expression \f$g(x,t)\f$ for the boundary conditions,
1315  * this function evaluates the boundary conditions for all boundaries
1316  * at time-level \a t.
1317  * @param time The time at which the boundary conditions
1318  * should be evaluated.
1319  * @param bndCondExpansions List of boundary expansions.
1320  * @param bndConditions Information about the boundary conditions.
1321  */
1322  void DisContField1D::v_EvaluateBoundaryConditions(
1323  const NekDouble time,
1324  const std::string varName,
1325  const NekDouble x2_in,
1326  const NekDouble x3_in)
1327  {
1328  int i;
1329 
1330  Array<OneD, NekDouble> x0(1);
1331  Array<OneD, NekDouble> x1(1);
1332  Array<OneD, NekDouble> x2(1);
1333 
1334  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1335  {
1336  if (time == 0.0 || m_bndConditions[i]->GetUserDefined() ==
1337  SpatialDomains::eTimeDependent ||
1338  m_bndConditions[i]->GetUserDefined() ==
1339  SpatialDomains::eQinflow ||
1340  m_bndConditions[i]->GetUserDefined() ==
1341  SpatialDomains::eRCRterminal )
1342  {
1343  m_bndCondExpansions[i]->GetCoords(x0, x1, x2);
1344 
1345  if (x2_in != NekConstants::kNekUnsetDouble && x3_in !=
1346  NekConstants::kNekUnsetDouble)
1347  {
1348  x1[0] = x2_in;
1349  x2[0] = x3_in;
1350  }
1351 
1352  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1353  SpatialDomains::eDirichlet)
1354  {
1355  m_bndCondExpansions[i]->SetCoeff(0,
1356  (boost::static_pointer_cast<SpatialDomains
1357  ::DirichletBoundaryCondition>(m_bndConditions[i])
1358  ->m_dirichletCondition).Evaluate(x0[0],x1[0],x2[0],time));
1359  m_bndCondExpansions[i]->SetPhys(0,m_bndCondExpansions[i]->GetCoeff(0));
1360  }
1361  else if (m_bndConditions[i]->GetBoundaryConditionType()
1362  == SpatialDomains::eNeumann)
1363  {
1364  m_bndCondExpansions[i]->SetCoeff(0,
1365  (boost::static_pointer_cast<SpatialDomains
1366  ::NeumannBoundaryCondition>(m_bndConditions[i])
1367  ->m_neumannCondition).Evaluate(x0[0],x1[0],x2[0],time));
1368  }
1369  else if (m_bndConditions[i]->GetBoundaryConditionType()
1370  == SpatialDomains::eRobin)
1371  {
1372  m_bndCondExpansions[i]->SetCoeff(0,
1373  (boost::static_pointer_cast<SpatialDomains
1374  ::RobinBoundaryCondition>(m_bndConditions[i])
1375  ->m_robinFunction).Evaluate(x0[0],x1[0],x2[0],time));
1376 
1377  m_bndCondExpansions[i]->SetPhys(0,
1378  (boost::static_pointer_cast<SpatialDomains
1379  ::RobinBoundaryCondition>(m_bndConditions[i])
1380  ->m_robinPrimitiveCoeff).Evaluate(x0[0],x1[0],x2[0],time));
1381 
1382  }
1383  else if (m_bndConditions[i]->GetBoundaryConditionType()
1384  == SpatialDomains::eNotDefined)
1385  {
1386  }
1387  else
1388  {
1389  ASSERTL0(false, "This type of BC not implemented yet");
1390  }
1391  }
1392  }
1393  }
1394 
1395  // Set up a list of element ids and edge ids that link to the
1396  // boundary conditions
1397  void DisContField1D::v_GetBoundaryToElmtMap(
1398  Array<OneD, int> &ElmtID, Array<OneD,int> &VertID)
1399  {
1400  map<int, int> VertGID;
1401  int i,n,id;
1402  int bid,cnt,Vid;
1403  int nbcs = m_bndConditions.num_elements();
1404 
1405  // make sure arrays are of sufficient length
1406  if (ElmtID.num_elements() != nbcs)
1407  {
1408  ElmtID = Array<OneD, int>(nbcs,-1);
1409  }
1410  else
1411  {
1412  fill(ElmtID.get(), ElmtID.get()+nbcs, -1);
1413  }
1414 
1415  if (VertID.num_elements() != nbcs)
1416  {
1417  VertID = Array<OneD, int>(nbcs);
1418  }
1419 
1420  // setup map of all global ids along boundary
1421  for (cnt = n = 0; n < m_bndCondExpansions.num_elements(); ++n)
1422  {
1423  Vid = m_bndCondExpansions[n]->GetExp(0)->GetGeom()->GetVertex(0)->GetVid();
1424  VertGID[Vid] = cnt++;
1425  }
1426 
1427  // Loop over elements and find verts that match;
1428  LocalRegions::ExpansionSharedPtr exp;
1429  for(cnt = n = 0; n < GetExpSize(); ++n)
1430  {
1431  exp = (*m_exp)[n];
1432  for(i = 0; i < exp->GetNverts(); ++i)
1433  {
1434  id = exp->GetGeom()->GetVid(i);
1435 
1436  if (VertGID.count(id) > 0)
1437  {
1438  bid = VertGID.find(id)->second;
1439  ASSERTL1(ElmtID[bid] == -1,"Edge already set");
1440  ElmtID[bid] = n;
1441  VertID[bid] = i;
1442  cnt ++;
1443  }
1444  }
1445  }
1446 
1447  ASSERTL1(cnt == nbcs,"Failed to visit all boundary condtiions");
1448  }
1449 
1450  /**
1451  * Search through the edge expansions and identify which ones
1452  * have Robin/Mixed type boundary conditions. If find a Robin
1453  * boundary then store the edge id of the boundary condition
1454  * and the array of points of the physical space boundary
1455  * condition which are hold the boundary condition primitive
1456  * variable coefficient at the quatrature points
1457  *
1458  * \return std map containing the robin boundary condition
1459  * info using a key of the element id
1460  *
1461  * There is a next member to allow for more than one Robin
1462  * boundary condition per element
1463  */
1464  map<int, RobinBCInfoSharedPtr> DisContField1D::v_GetRobinBCInfo(void)
1465  {
1466  int i;
1467  map<int, RobinBCInfoSharedPtr> returnval;
1468  Array<OneD, int> ElmtID,VertID;
1469  GetBoundaryToElmtMap(ElmtID,VertID);
1470 
1471  for (i = 0; i < m_bndCondExpansions.num_elements(); ++i)
1472  {
1473  MultiRegions::ExpListSharedPtr locExpList;
1474 
1475  if (m_bndConditions[i]->GetBoundaryConditionType() ==
1476  SpatialDomains::eRobin)
1477  {
1478  int elmtid;
1479  Array<OneD, NekDouble> Array_tmp;
1480 
1481  locExpList = m_bndCondExpansions[i];
1482 
1483  RobinBCInfoSharedPtr rInfo =
1484  MemoryManager<RobinBCInfo>::
1485  AllocateSharedPtr(
1486  VertID[i],Array_tmp = locExpList->GetPhys());
1487 
1488  elmtid = ElmtID[i];
1489  // make link list if necessary (not likely in
1490  // 1D but needed in 2D & 3D)
1491  if(returnval.count(elmtid) != 0)
1492  {
1493  rInfo->next = returnval.find(elmtid)->second;
1494  }
1495  returnval[elmtid] = rInfo;
1496  }
1497  }
1498 
1499  return returnval;
1500  }
1501  } // end of namespace
1502 } //end of namespace
Generated on Sun Mar 15 2015 00:11:30 for Nektar++ by   doxygen 1.8.1.2