Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpList1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList1D.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: Expansion list 1D definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iomanip>
37 #include <MultiRegions/ExpList1D.h>
38 #include <LibUtilities/Polylib/Polylib.h>
39 #include <LocalRegions/SegExp.h>
40 #include <LocalRegions/QuadExp.h>
41 #include <LocalRegions/Expansion2D.h>
42 #include <SpatialDomains/MeshGraph2D.h>
43 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
44 #include <LibUtilities/Foundations/Interp.h>
45 
46 namespace Nektar
47 {
48  namespace MultiRegions
49  {
50  /**
51  * @class ExpList1D
52  * This multi-elemental expansion, which does not exhibit any coupling
53  * between the expansion on the separate elements, can be formulated
54  * as, \f[u^{\delta}(x_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
55  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(x_i).\f]
56  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
57  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
58  *
59  * Instances of this class may be optionally constructed which use
60  * generalised segment expansions (LocalRegions#GenSegExp), rather than
61  * the standard segment expansions (LocalRegions#SegExp).
62  * LocalRegions#GenSegExp provides additional spatial
63  * data including segment normals and is enabled using the \a
64  * UseGenSegExp flag.
65  *
66  * This class inherits all its variables and member functions from the
67  * base class MultiRegions#ExpList.
68  */
69 
70  /**
71  * Assumes use of standard segment expansions only. All data storage
72  * areas are initialised to empty arrays by the default ExpList
73  * constructor.
74  */
75  ExpList1D::ExpList1D():
76  ExpList()
77  {
78  SetExpType(e1D);
79  }
80 
81 
82  /**
83  * Creates an identical copy of another ExpList1D object.
84  */
85  ExpList1D::ExpList1D(const ExpList1D &In, const bool DeclareCoeffPhysArrays):
86  ExpList(In,DeclareCoeffPhysArrays)
87  {
88  SetExpType(e1D);
89  }
90 
91 
92  /**
93  * After initialising the data inherited through MultiRegions#ExpList,
94  * populate the expansion list from the segments defined in the supplied
95  * SpatialDomains#MeshGraph1D. All expansions in the graph are defined
96  * using the same LibUtilities#BasisKey which overrides that specified
97  * in \a graph1D.
98  *
99  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
100  * for details.
101  * @deprecated This constructor is no longer required since
102  * the basis key is now available from the graph.
103  * @param Ba BasisKey describing quadrature points and
104  * number of modes.
105  * @param graph1D Domain and expansion definitions.
106  */
107  ExpList1D::ExpList1D(const LibUtilities::SessionReaderSharedPtr &pSession,
108  const LibUtilities::BasisKey &Ba,
109  const SpatialDomains::MeshGraphSharedPtr &graph1D):
110  ExpList(pSession,graph1D)
111  {
112  SetExpType(e1D);
113 
114  int id=0;
115  LocalRegions::SegExpSharedPtr seg;
116  SpatialDomains::SegGeomSharedPtr SegmentGeom;
117 
118  const SpatialDomains::ExpansionMap &expansions
119  = graph1D->GetExpansions();
120 
121  // For each element in the mesh, create a segment expansion using
122  // the supplied BasisKey and segment geometry.
123  SpatialDomains::ExpansionMap::const_iterator expIt;
124  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
125  {
126  if ((SegmentGeom = boost::dynamic_pointer_cast<
127  SpatialDomains::SegGeom>(
128  expIt->second->m_geomShPtr)))
129  {
130  seg = MemoryManager<LocalRegions::SegExp>
131  ::AllocateSharedPtr(Ba,SegmentGeom);
132  seg->SetElmtId(id++);
133  (*m_exp).push_back(seg);
134  }
135  else
136  {
137  ASSERTL0(false,"dynamic cast to a SegGeom failed");
138  }
139  }
140 
141  // Setup Default optimisation information.
142  int nel = GetExpSize();
143  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
144  ::AllocateSharedPtr(nel);
145 
146  // Allocate storage for data and populate element offset lists.
147  SetCoeffPhysOffsets();
148 
149  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
150  m_phys = Array<OneD, NekDouble>(m_npoints);
151 
152  ReadGlobalOptimizationParameters();
153  }
154 
155 
156  /**
157  * Given a mesh \a graph1D, containing information about the domain and
158  * the spectral/hp element expansion, this constructor fills the list
159  * of local expansions \texttt{m_exp} with the proper expansions,
160  * calculates the total number of quadrature points \f$x_i\f$ and local
161  * expansion coefficients \f$\hat{u}^e_n\f$ and allocates memory for
162  * the arrays #m_coeffs and #m_phys.
163  *
164  * For each element its corresponding LibUtilities#BasisKey is
165  * retrieved and this is used to construct either a standard segment
166  * (LocalRegions#SegExp) or a generalised segment
167  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
168  * Finally, ExpList#SetCoeffPhys is called to initialise the data
169  * storage areas and set up the offset arrays.
170  *
171  * @param graph1D A mesh, containing information about the
172  * domain and the spectral/hp element expansion.
173  * @param UseGenSegExp If true, create general segment expansions
174  * instead of just normal segment expansions.
175  */
176  ExpList1D::ExpList1D(const LibUtilities::SessionReaderSharedPtr &pSession,
177  const SpatialDomains::MeshGraphSharedPtr &graph1D,
178  const bool DeclareCoeffPhysArrays):
179  ExpList(pSession,graph1D)
180  {
181  SetExpType(e1D);
182 
183  int id=0;
184  LocalRegions::SegExpSharedPtr seg;
185  SpatialDomains::SegGeomSharedPtr SegmentGeom;
186 
187  // Retrieve the list of expansions
188  const SpatialDomains::ExpansionMap &expansions
189  = graph1D->GetExpansions();
190 
191  // Process each expansion in the graph
192  SpatialDomains::ExpansionMap::const_iterator expIt;
193  for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
194  {
195  // Retrieve basis key from expansion
196  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
197 
198  if ((SegmentGeom = boost::dynamic_pointer_cast<
199  SpatialDomains::SegGeom>(
200  expIt->second->m_geomShPtr)))
201  {
202  seg = MemoryManager<LocalRegions::SegExp>
203  ::AllocateSharedPtr(bkey, SegmentGeom);
204 
205  // Assign next ID
206  seg->SetElmtId(id++);
207 
208  // Add the expansion
209  (*m_exp).push_back(seg);
210  }
211  else
212  {
213  ASSERTL0(false,"dynamic cast to a SegGeom failed");
214  }
215  }
216 
217  // Setup Default optimisation information.
218  int nel = GetExpSize();
219  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
220  ::AllocateSharedPtr(nel);
221 
222  // set up offset arrays.
223  SetCoeffPhysOffsets();
224 
225  if(DeclareCoeffPhysArrays)
226  {
227  // Set up m_coeffs, m_phys.
228  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
229  m_phys = Array<OneD, NekDouble>(m_npoints);
230  }
231 
232  ReadGlobalOptimizationParameters();
233  }
234 
235 
236 
237  /**
238  * Given a mesh \a graph1D, and the spectral/hp element
239  * expansion as well as and separate information about a \a
240  * domain, this constructor fills the list of local
241  * expansions \texttt{m_exp} with the proper expansions,
242  * calculates the total number of quadrature points \f$x_i\f$
243  * and local expansion coefficients \f$\hat{u}^e_n\f$ and
244  * allocates memory for the arrays #m_coeffs and #m_phys.
245  *
246  * For each element its corresponding LibUtilities#BasisKey is
247  * retrieved and this is used to construct either a standard segment
248  * (LocalRegions#SegExp) or a generalised segment
249  * (LocalRegions#GenSegExp) which is stored in the list #m_exp.
250  * Finally, ExpList#SetCoeffPhys is called to initialise the data
251  * storage areas and set up the offset arrays.
252  *
253  * @param graph1D A mesh, containing information about the
254  * domain and the spectral/hp element expansion.
255  * @param UseGenSegExp If true, create general segment expansions
256  * instead of just normal segment expansions.
257  */
258  ExpList1D::ExpList1D(const LibUtilities::SessionReaderSharedPtr &pSession,
259  const SpatialDomains::MeshGraphSharedPtr &graph1D,
260  const SpatialDomains::CompositeMap &domain,
261  const bool DeclareCoeffPhysArrays,
262  const std::string var,
263  bool SetToOneSpaceDimension):
264  ExpList(pSession,graph1D)
265  {
266  int j,id=0;
267  LocalRegions::SegExpSharedPtr seg;
268  SpatialDomains::SegGeomSharedPtr SegmentGeom;
269  SpatialDomains::Composite comp;
270  SpatialDomains::CompositeMap::const_iterator compIt;
271 
272  // Retrieve the list of expansions
273  const SpatialDomains::ExpansionMap &expansions
274  = graph1D->GetExpansions(var);
275 
276  // Process each composite region in domain
277  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
278  {
279  comp = compIt->second;
280 
281  // Process each expansion in the graph
282  for(j = 0; j < compIt->second->size(); ++j)
283  {
284  SpatialDomains::ExpansionMap::const_iterator expIt;
285 
286  if((SegmentGeom = boost::dynamic_pointer_cast<
287  SpatialDomains::SegGeom>(
288  (*compIt->second)[j])))
289  {
290  // Retrieve basis key from expansion and define expansion
291  if((expIt = expansions.find(SegmentGeom->GetGlobalID())) != expansions.end())
292  {
293  LibUtilities::BasisKey bkey = expIt->second->m_basisKeyVector[0];
294 
295  if(SetToOneSpaceDimension)
296  {
297  SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
298  SegmentGeom->GenerateOneSpaceDimGeom();
299 
300  seg = MemoryManager<LocalRegions::SegExp>
301  ::AllocateSharedPtr(bkey, OneDSegmentGeom);
302  }
303  else
304  {
305  seg = MemoryManager<LocalRegions::SegExp>
306  ::AllocateSharedPtr(bkey, SegmentGeom);
307  }
308  }
309  else
310  {
311  ASSERTL0(false,"Failed to find basis key");
312  }
313 
314  }
315  else
316  {
317  ASSERTL0(false,"Failed to dynamic cast geometry to SegGeom");
318  }
319 
320 
321  // Assign next ID
322  seg->SetElmtId(id++);
323 
324  // Add the expansion
325  (*m_exp).push_back(seg);
326  }
327  }
328 
329  // Setup Default optimisation information.
330  int nel = GetExpSize();
331  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
332  ::AllocateSharedPtr(nel);
333 
334  // set up offset arrays.
335  SetCoeffPhysOffsets();
336 
337  if(DeclareCoeffPhysArrays)
338  {
339  // Set up m_coeffs, m_phys.
340  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
341  m_phys = Array<OneD, NekDouble>(m_npoints);
342  }
343 
344  ReadGlobalOptimizationParameters();
345  }
346 
347 
348  /**
349  * Fills the list of local expansions with the segments from the 2D
350  * mesh specified by \a domain. This CompositeMap contains a list of
351  * Composites which define the Neumann boundary.
352  * @see ExpList1D#ExpList1D(SpatialDomains::MeshGraph1D&, bool)
353  * for details.
354  * @param domain A domain, comprising of one or more composite
355  * regions,
356  * @param graph2D A mesh, containing information about the
357  * domain and the spectral/hp element expansion.
358  * @param UseGenSegExp If true, create general segment expansions
359  * instead of just normal segment expansions.
360  */
361  ExpList1D::ExpList1D(const SpatialDomains::CompositeMap &domain,
362  const SpatialDomains::MeshGraphSharedPtr &graph2D,
363  const bool DeclareCoeffPhysArrays,
364  const std::string variable):
365  ExpList()
366  {
367  SetExpType(e1D);
368 
369  int j, id=0;
370  SpatialDomains::Composite comp;
371  SpatialDomains::CompositeMap::const_iterator compIt;
372  SpatialDomains::SegGeomSharedPtr SegmentGeom;
373  LocalRegions::SegExpSharedPtr seg;
374 
375  // Process each composite region.
376  for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
377  {
378  comp = compIt->second;
379  // Process each expansion in the region.
380  for(j = 0; j < compIt->second->size(); ++j)
381  {
382  if((SegmentGeom = boost::dynamic_pointer_cast<
383  SpatialDomains::SegGeom>(
384  (*compIt->second)[j])))
385  {
386  // Retrieve the basis key from the expansion.
387  LibUtilities::BasisKey bkey
388  = boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(graph2D)->GetEdgeBasisKey(SegmentGeom, variable);
389 
390  seg = MemoryManager<LocalRegions::SegExp>
391  ::AllocateSharedPtr(bkey, SegmentGeom);
392 
393  // Add the segment to the expansion list.
394  seg->SetElmtId(id++);
395  (*m_exp).push_back(seg);
396  }
397  else
398  {
399  ASSERTL0(false,"dynamic cast to a SegGeom failed");
400  }
401  }
402 
403  }
404 
405  // Setup Default optimisation information.
406  int nel = GetExpSize();
407  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
408  ::AllocateSharedPtr(nel);
409 
410  // Allocate storage for data and populate element offset lists.
411  SetCoeffPhysOffsets();
412 
413  // Set up m_coeffs, m_phys.
414  if(DeclareCoeffPhysArrays)
415  {
416  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
417  m_phys = Array<OneD, NekDouble>(m_npoints);
418  }
419  }
420 
421  /**
422  * Store expansions for the trace space expansions used in
423  * DisContField2D.
424  *
425  * @param bndConstraint Array of ExpList1D objects each containing a
426  * 1D spectral/hp element expansion on a single
427  * boundary region.
428  * @param bndCond Array of BoundaryCondition objects which contain
429  * information about the boundary conditions on the
430  * different boundary regions.
431  * @param locexp Complete domain expansion list.
432  * @param graph2D 2D mesh corresponding to the expansion list.
433  * @param periodicEdges List of periodic edges.
434  * @param UseGenSegExp If true, create general segment expansions
435  * instead of just normal segment expansions.
436  */
437  ExpList1D::ExpList1D(
438  const LibUtilities::SessionReaderSharedPtr &pSession,
439  const Array<OneD,const ExpListSharedPtr> &bndConstraint,
440  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr> &bndCond,
441  const LocalRegions::ExpansionVector &locexp,
442  const SpatialDomains::MeshGraphSharedPtr &graph2D,
443  const PeriodicMap &periodicEdges,
444  const bool DeclareCoeffPhysArrays,
445  const std::string variable):
446  ExpList()
447  {
448  int i, j, id, elmtid = 0;
449  set<int> edgesDone;
450 
451  SpatialDomains::Geometry1DSharedPtr segGeom;
452  SpatialDomains::Geometry2DSharedPtr ElGeom;
453  LocalRegions::SegExpSharedPtr seg;
454  LocalRegions::SegExpSharedPtr seg_tmp;
455  LocalRegions::Expansion1DSharedPtr exp1D;
456  LocalRegions::Expansion2DSharedPtr exp2D;
457 
458  SetExpType(e1D);
459 
460  map<int,int> EdgeDone;
461  map<int,int> NormalSet;
462 
463  LocalRegions::SegExpSharedPtr Seg;
464 
465  // First loop over boundary conditions to renumber
466  // Dirichlet boundaries
467  for(i = 0; i < bndCond.num_elements(); ++i)
468  {
469  if(bndCond[i]->GetBoundaryConditionType()
470  == SpatialDomains::eDirichlet)
471  {
472  for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
473  {
474  LibUtilities::BasisKey bkey = bndConstraint[i]
475  ->GetExp(j)->GetBasis(0)->GetBasisKey();
476  exp1D = bndConstraint[i]->GetExp(j)->
477  as<LocalRegions::Expansion1D>();
478  segGeom = exp1D->GetGeom1D();
479 
480  seg = MemoryManager<LocalRegions::SegExp>
481  ::AllocateSharedPtr(bkey, segGeom);
482  edgesDone.insert(segGeom->GetEid());
483 
484  seg->SetElmtId(elmtid++);
485  (*m_exp).push_back(seg);
486  }
487  }
488  }
489 
490  map<int, pair<SpatialDomains::Geometry1DSharedPtr,
491  LibUtilities::BasisKey> > edgeOrders;
492  map<int, pair<SpatialDomains::Geometry1DSharedPtr,
493  LibUtilities::BasisKey> >::iterator it;
494 
495  for(i = 0; i < locexp.size(); ++i)
496  {
497  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
498 
499  for(j = 0; j < locexp[i]->GetNedges(); ++j)
500  {
501  segGeom = exp2D->GetGeom2D()->GetEdge(j);
502  id = segGeom->GetEid();
503  // Ignore Dirichlet edges
504  if (edgesDone.count(id) != 0)
505  {
506  continue;
507  }
508 
509  it = edgeOrders.find(id);
510 
511  if (it == edgeOrders.end())
512  {
513  edgeOrders.insert(std::make_pair(id, std::make_pair(
514  segGeom, locexp[i]->DetEdgeBasisKey(j))));
515  }
516  else // variable modes/points
517  {
518  LibUtilities::BasisKey edge
519  = locexp[i]->DetEdgeBasisKey(j);
520  LibUtilities::BasisKey existing
521  = it->second.second;
522 
523  int np1 = edge .GetNumPoints();
524  int np2 = existing.GetNumPoints();
525  int nm1 = edge .GetNumModes ();
526  int nm2 = existing.GetNumModes ();
527 
528  if (np2 >= np1 && nm2 >= nm1)
529  {
530  continue;
531  }
532  else if (np2 < np1 && nm2 < nm1)
533  {
534  it->second.second = edge;
535  }
536  else
537  {
538  ASSERTL0(false,
539  "inappropriate number of points/modes (max "
540  "num of points is not set with max order)");
541  }
542  }
543  }
544  }
545 
546  LibUtilities::CommSharedPtr vComm =
547  pSession->GetComm()->GetRowComm();
548  int nproc = vComm->GetSize(); // number of processors
549  int edgepr = vComm->GetRank(); // ID processor
550 
551  if (nproc > 1)
552  {
553  int eCnt = 0;
554 
555  // Count the number of edges on each partition
556  for(i = 0; i < locexp.size(); ++i)
557  {
558  eCnt += locexp[i]->GetNedges();
559  }
560 
561  // Set up the offset and the array that will contain the list of
562  // edge IDs, then reduce this across processors.
563  Array<OneD, int> edgesCnt(nproc, 0);
564  edgesCnt[edgepr] = eCnt;
565  vComm->AllReduce(edgesCnt, LibUtilities::ReduceSum);
566 
567  // Set up offset array.
568  int totEdgeCnt = Vmath::Vsum(nproc, edgesCnt, 1);
569  Array<OneD, int> eTotOffsets(nproc,0);
570  for (i = 1; i < nproc; ++i)
571  {
572  eTotOffsets[i] = eTotOffsets[i-1] + edgesCnt[i-1];
573  }
574 
575  // Local list of the edges per element
576  Array<OneD, int> EdgesTotID(totEdgeCnt, 0);
577  Array<OneD, int> EdgesTotNm(totEdgeCnt, 0);
578  Array<OneD, int> EdgesTotPnts(totEdgeCnt, 0);
579 
580  int cntr = eTotOffsets[edgepr];
581 
582  for(i = 0; i < locexp.size(); ++i)
583  {
584  exp2D = locexp[i]->as<LocalRegions::Expansion2D>();
585 
586  int nedges = locexp[i]->GetNedges();
587 
588  for(j = 0; j < nedges; ++j, ++cntr)
589  {
590  LibUtilities::BasisKey bkeyEdge =
591  locexp[i]->DetEdgeBasisKey(j);
592  EdgesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
593  EdgesTotNm [cntr] = bkeyEdge.GetNumModes();
594  EdgesTotPnts[cntr] = bkeyEdge.GetNumPoints();
595  }
596  }
597 
598  vComm->AllReduce(EdgesTotID, LibUtilities::ReduceSum);
599  vComm->AllReduce(EdgesTotNm, LibUtilities::ReduceSum);
600  vComm->AllReduce(EdgesTotPnts, LibUtilities::ReduceSum);
601 
602  for (i = 0; i < totEdgeCnt; ++i)
603  {
604  it = edgeOrders.find(EdgesTotID[i]);
605 
606  if (it == edgeOrders.end())
607  {
608  continue;
609  }
610 
611  LibUtilities::BasisKey existing
612  = it->second.second;
613  LibUtilities::BasisKey edge(
614  existing.GetBasisType(), EdgesTotNm[i],
615  LibUtilities::PointsKey(EdgesTotPnts[i],
616  existing.GetPointsType()));
617 
618 
619  int np1 = edge .GetNumPoints();
620  int np2 = existing.GetNumPoints();
621  int nm1 = edge .GetNumModes ();
622  int nm2 = existing.GetNumModes ();
623 
624  if (np2 >= np1 && nm2 >= nm1)
625  {
626  continue;
627  }
628  else if (np2 < np1 && nm2 < nm1)
629  {
630  it->second.second = edge;
631  }
632  else
633  {
634  ASSERTL0(false,
635  "inappropriate number of points/modes (max "
636  "num of points is not set with max order)");
637  }
638  }
639  }
640 
641  for (it = edgeOrders.begin(); it != edgeOrders.end(); ++it)
642  {
643  seg = MemoryManager<LocalRegions::SegExp>
644  ::AllocateSharedPtr(it->second.second, it->second.first);
645  seg->SetElmtId(elmtid++);
646  (*m_exp).push_back(seg);
647  }
648 
649  // Setup Default optimisation information.
650  int nel = GetExpSize();
651  m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
652  ::AllocateSharedPtr(nel);
653 
654  // Set up offset information and array sizes
655  SetCoeffPhysOffsets();
656 
657  // Set up m_coeffs, m_phys.
658  if(DeclareCoeffPhysArrays)
659  {
660  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
661  m_phys = Array<OneD, NekDouble>(m_npoints);
662  }
663  }
664 
665  /**
666  * Each expansion (local element) is processed in turn to
667  * determine the number of coefficients and physical data
668  * points it contributes to the domain. Three arrays,
669  * #m_coeff_offset, #m_phys_offset and #m_offset_elmt_id, are
670  * also initialised and updated to store the data offsets of
671  * each element in the #m_coeffs and #m_phys arrays, and the
672  * element id that each consecutive block is associated
673  * respectively.
674  */
675  void ExpList1D::SetCoeffPhysOffsets()
676  {
677  int i;
678 
679  // Set up offset information and array sizes
680  m_coeff_offset = Array<OneD,int>(m_exp->size());
681  m_phys_offset = Array<OneD,int>(m_exp->size());
682  m_offset_elmt_id = Array<OneD,int>(m_exp->size());
683 
684  m_ncoeffs = m_npoints = 0;
685 
686  for(i = 0; i < m_exp->size(); ++i)
687  {
688  m_coeff_offset[i] = m_ncoeffs;
689  m_phys_offset [i] = m_npoints;
690  m_offset_elmt_id[i] = i;
691  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
692  m_npoints += (*m_exp)[i]->GetTotPoints();
693  }
694  }
695 
696  /**
697  *
698  */
699  ExpList1D::~ExpList1D()
700  {
701  }
702 
703  /**
704  * To perform post-processing on the entire domain use \a elmtId = 0.
705  * @param kernel The post-processing kernel.
706  * @param inarray The set of evaluation points.
707  * @param outarray Contains the resulting post-processed
708  * solution for element \a elmId.
709  * @param h The mesh spacing.
710  * @param elmId Optionally specifies which element to perform
711  * the post-processing on (0=whole domain).
712  */
713  void ExpList1D::PostProcess(LibUtilities::KernelSharedPtr kernel,
714  Array<OneD,NekDouble> &inarray,
715  Array<OneD,NekDouble> &outarray,
716  NekDouble h,
717  int elmId)
718 
719  {
720  int i,j,r;
721 
722  // get the local element expansion of the elmId element
723  StdRegions::StdExpansionSharedPtr elmExp = GetExp(elmId);
724 
725  // Get the quadrature points and weights required for integration
726  int quad_npoints = elmExp->GetTotPoints();
727  LibUtilities::PointsKey quadPointsKey(quad_npoints,
728  elmExp->GetPointsType(0));
729  Array<OneD,NekDouble> quad_points
730  = LibUtilities::PointsManager()[quadPointsKey]->GetZ();
731  Array<OneD,NekDouble> quad_weights
732  = LibUtilities::PointsManager()[quadPointsKey]->GetW();
733 
734  // Declare variable for the local kernel breaks
735  int kernel_width = kernel->GetKernelWidth();
736  Array<OneD,NekDouble> local_kernel_breaks(kernel_width+1);
737 
738  // Declare variable for the transformed quadrature points
739  Array<OneD,NekDouble> mapped_quad_points(quad_npoints);
740 
741  // For each evaluation point
742  for(i = 0; i < inarray.num_elements(); i++)
743  {
744  // Move the center of the kernel to the current point
745  kernel->MoveKernelCenter(inarray[i],local_kernel_breaks);
746 
747  // Find the mesh breaks under the kernel support
748  Array<OneD,NekDouble> mesh_breaks;
749  kernel->FindMeshUnderKernel(local_kernel_breaks,h,mesh_breaks);
750 
751  // Sort the total breaks for integration purposes
752  int total_nbreaks = local_kernel_breaks.num_elements() +
753  mesh_breaks.num_elements();
754  // number of the total breaks
755  Array<OneD,NekDouble> total_breaks(total_nbreaks);
756  kernel->Sort(local_kernel_breaks,mesh_breaks,total_breaks);
757 
758  // Integrate the product of kernel and function over the total
759  // breaks
760  NekDouble integral_value = 0.0;
761  for(j = 0; j < total_breaks.num_elements()-1; j++)
762  {
763  NekDouble a = total_breaks[j];
764  NekDouble b = total_breaks[j+1];
765 
766  // Map the quadrature points to the appropriate interval
767  for(r = 0; r < quad_points.num_elements(); r++)
768  {
769  mapped_quad_points[r]
770  = (quad_points[r] + 1.0) * 0.5 * (b - a) + a;
771  }
772 
773  // Evaluate the function at the transformed quadrature
774  // points
775  Array<OneD,NekDouble> u_value(quad_npoints);
776  Array<OneD,NekDouble> coeffs = GetCoeffs();
777 
778  PeriodicEval(coeffs,mapped_quad_points,h,
779  elmExp->GetBasisNumModes(0),u_value);
780 
781  // Evaluate the kernel at the transformed quadrature points
782  Array<OneD,NekDouble> k_value(quad_npoints);
783  kernel->EvaluateKernel(mapped_quad_points,h,k_value);
784 
785  // Integrate
786  for(r = 0; r < quad_npoints; r++)
787  {
788  integral_value += (b - a) * 0.5 * k_value[r]
789  * u_value[r] * quad_weights[r];
790  }
791  }
792  outarray[i] = integral_value/h;
793  }
794  }
795 
796 
797  /**
798  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of an expansion,
799  * periodically evaluate the spectral/hp expansion
800  * \f$u^{\delta}(\boldsymbol{x})\f$ at arbitrary points.
801  * @param inarray1 An array of size \f$N_{\mathrm{eof}}\f$
802  * containing the local coefficients
803  * \f$\hat{u}_n^e\f$.
804  * @param inarray2 Contains the set of evaluation points.
805  * @param h The mesh spacing.
806  * @param nmodes The number of polynomial modes for each element
807  * (we consider that each element has the same
808  * number of polynomial modes).
809  * @param outarray Contains the resulting values at the
810  * evaluation points
811  */
812  void ExpList1D::PeriodicEval(Array<OneD,NekDouble> &inarray1,
813  Array<OneD,NekDouble> &inarray2,
814  NekDouble h, int nmodes,
815  Array<OneD,NekDouble> &outarray)
816  {
817  int i,j,r;
818 
819  // Get the number of elements in the domain
820  int num_elm = GetExpSize();
821 
822  // initializing the outarray
823  for(i = 0; i < outarray.num_elements(); i++)
824  {
825  outarray[i] = 0.0;
826  }
827 
828  // Make a copy for further modification
829  int x_size = inarray2.num_elements();
830  Array<OneD,NekDouble> x_values_cp(x_size);
831 
832  // Determining the element to which the x belongs
833  Array<OneD,int> x_elm(x_size);
834  for(i = 0; i < x_size; i++ )
835  {
836  x_elm[i] = (int)floor(inarray2[i]/h);
837  }
838 
839  // Clamp indices periodically
840  for(i = 0; i < x_size; i++)
841  {
842  while(x_elm[i] < 0)
843  {
844  x_elm[i] += num_elm;
845  }
846  while(x_elm[i] >= num_elm)
847  {
848  x_elm[i] -= num_elm ;
849  }
850  }
851 
852  // Map the values of x to [-1 1] on its interval
853  for(i = 0; i < x_size; i++)
854  {
855  x_values_cp[i] = (inarray2[i]/h - floor(inarray2[i]/h))*2 - 1.0;
856  }
857 
858  // Evaluate the jocobi polynomials
859  // (Evaluating the base at some points other than the quadrature
860  // points). Should it be added to the base class????
861  Array<TwoD,NekDouble> jacobi_poly(nmodes,x_size);
862  for(i = 0; i < nmodes; i++)
863  {
864  Polylib::jacobfd(x_size,x_values_cp.get(),
865  jacobi_poly.get()+i*x_size,NULL,i,0.0,0.0);
866  }
867 
868  // Evaluate the function values
869  for(r = 0; r < nmodes; r++)
870  {
871  for(j = 0; j < x_size; j++)
872  {
873  int index = ((x_elm[j])*nmodes)+r;
874  outarray[j] += inarray1[index]*jacobi_poly[r][j];
875  }
876  }
877 
878  }
879 
880 
881  /**
882  * Sets up the normals on all edges of expansions in the domain.
883  * @param locexp Complete list of domain expansions.
884  */
885 // void ExpList1D::SetUpPhysNormals(
886 // const StdRegions::StdExpansionVector &locexp)
887 // {
888 // map<int, int> EdgeGID;
889 // int i,cnt,n,id;
890 //
891 // // setup map of all global ids along boundary
892 // for(cnt = i = 0; i < (*m_exp).size(); ++i)
893 // {
894 // id = (*m_exp)[i]->GetGeom1D()->GetEid();
895 // EdgeGID[id] = cnt++;
896 // }
897 //
898 // // Loop over elements and find edges that match;
899 // for(cnt = n = 0; n < locexp.size(); ++n)
900 // {
901 // for(i = 0; i < locexp[n]->GetNedges(); ++i)
902 // {
903 // id = locexp[n]->GetGeom2D()->GetEid(i);
904 //
905 // if(EdgeGID.count(id) > 0)
906 // {
907 // (*m_exp)[EdgeGID.find(id)->second]
908 // ->SetUpPhysNormals(locexp[n],i);
909 // }
910 // }
911 // }
912 // }
913 
914  //croth
915  void ExpList1D::v_SetUpPhysNormals()
916  {
917  int i, j;
918  for (i = 0; i < m_exp->size(); ++i)
919  {
920  for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
921  {
922  (*m_exp)[i]->ComputeVertexNormal(j);
923  }
924  }
925  }
926 
927  /**
928  * Upwind the left and right states given by the Arrays Fwd and Bwd
929  * using the vector quantity Vec and ouput the upwinded value in the
930  * array upwind.
931  *
932  * @param Vec Velocity field.
933  * @param Fwd Left state.
934  * @param Bwd Right state.
935  * @param Upwind Output vector.
936  */
937  void ExpList1D::v_Upwind(
938  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
939  const Array<OneD, const NekDouble> &Fwd,
940  const Array<OneD, const NekDouble> &Bwd,
941  Array<OneD, NekDouble> &Upwind)
942  {
943  int i,j,k,e_npoints,offset;
944  Array<OneD,NekDouble> normals;
945  NekDouble Vn;
946 
947  // Assume whole array is of same coordimate dimension
948  int coordim = GetCoordim(0);
949 
950  ASSERTL1(Vec.num_elements() >= coordim,
951  "Input vector does not have sufficient dimensions to "
952  "match coordim");
953 
954  // Process each expansion
955  for(i = 0; i < m_exp->size(); ++i)
956  {
957  // Get the number of points in the expansion and the normals.
958  e_npoints = (*m_exp)[i]->GetNumPoints(0);
959  normals = (*m_exp)[i]->GetPhysNormals();
960 
961  // Get the physical data offset of the expansion in m_phys.
962  offset = m_phys_offset[i];
963 
964  // Compute each data point.
965  for(j = 0; j < e_npoints; ++j)
966  {
967  // Calculate normal velocity.
968  Vn = 0.0;
969  for(k = 0; k < coordim; ++k)
970  {
971  Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
972  }
973 
974  // Upwind based on direction of normal velocity.
975  if(Vn > 0.0)
976  {
977  Upwind[offset + j] = Fwd[offset + j];
978  }
979  else
980  {
981  Upwind[offset + j] = Bwd[offset + j];
982  }
983  }
984  }
985  }
986 
987  /**
988  * One-dimensional upwind.
989  * \see ExpList1D::Upwind(
990  * const Array<OneD, const Array<OneD, NekDouble> >,
991  * const Array<OneD, const NekDouble>,
992  * const Array<OneD, const NekDouble>,
993  * Array<OneD, NekDouble>, int)
994  *
995  * @param Vn Velocity field.
996  * @param Fwd Left state.
997  * @param Bwd Right state.
998  * @param Upwind Output vector.
999  */
1000  void ExpList1D::v_Upwind(
1001  const Array<OneD, const NekDouble> &Vn,
1002  const Array<OneD, const NekDouble> &Fwd,
1003  const Array<OneD, const NekDouble> &Bwd,
1004  Array<OneD, NekDouble> &Upwind)
1005  {
1006  int i,j,e_npoints,offset;
1007  Array<OneD,NekDouble> normals;
1008 
1009  // Process each expansion.
1010  for(i = 0; i < m_exp->size(); ++i)
1011  {
1012  // Get the number of points and the data offset.
1013  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1014  offset = m_phys_offset[i];
1015 
1016  // Process each point in the expansion.
1017  for(j = 0; j < e_npoints; ++j)
1018  {
1019  // Upwind based on one-dimensional velocity.
1020  if(Vn[offset + j] > 0.0)
1021  {
1022  Upwind[offset + j] = Fwd[offset + j];
1023  }
1024  else
1025  {
1026  Upwind[offset + j] = Bwd[offset + j];
1027  }
1028  }
1029  }
1030  }
1031 
1032 
1033  /**
1034  * For each local element, copy the normals stored in the element list
1035  * into the array \a normals.
1036  * @param normals Multidimensional array in which to copy normals
1037  * to. Must have dimension equal to or larger than
1038  * the spatial dimension of the elements.
1039  */
1040  void ExpList1D::v_GetNormals(
1041  Array<OneD, Array<OneD, NekDouble> > &normals)
1042  {
1043  int i,j,k,e_npoints,offset;
1044  SpatialDomains::Geometry1DSharedPtr segGeom;
1045  Array<OneD,Array<OneD,NekDouble> > locnormals;
1046  Array<OneD,Array<OneD,NekDouble> > locnormals2;
1047  Array<OneD,Array<OneD,NekDouble> > Norms;
1048  // Assume whole array is of same coordinate dimension
1049  int coordim = GetCoordim(0);
1050 
1051  ASSERTL1(normals.num_elements() >= coordim,
1052  "Output vector does not have sufficient dimensions to "
1053  "match coordim");
1054 
1055  for (i = 0; i < m_exp->size(); ++i)
1056  {
1057  LocalRegions::Expansion1DSharedPtr loc_exp = (*m_exp)[i]->as<LocalRegions::Expansion1D>();
1058 
1059  LocalRegions::Expansion2DSharedPtr loc_elmt =
1060  loc_exp->GetLeftAdjacentElementExp();
1061 
1062  int edgeNumber = loc_exp->GetLeftAdjacentElementEdge();
1063 
1064  // Get the number of points and normals for this expansion.
1065  e_npoints = (*m_exp)[i]->GetNumPoints(0);
1066 
1067  locnormals = loc_elmt->GetEdgeNormal(edgeNumber);
1068  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
1069  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
1070 
1071  if (e_nmodes != loc_nmodes)
1072  {
1073  if (loc_exp->GetRightAdjacentElementEdge() >= 0)
1074  {
1075  LocalRegions::Expansion2DSharedPtr loc_elmt =
1076  loc_exp->GetRightAdjacentElementExp();
1077 
1078  int EdgeNumber = loc_exp->GetRightAdjacentElementEdge();
1079  // Serial case: right element is connected so we can
1080  // just grab that normal.
1081  locnormals = loc_elmt->GetEdgeNormal(EdgeNumber);
1082 
1083  offset = m_phys_offset[i];
1084 
1085  // Process each point in the expansion.
1086  for (j = 0; j < e_npoints; ++j)
1087  {
1088  // Process each spatial dimension and copy the values
1089  // into the output array.
1090  for (k = 0; k < coordim; ++k)
1091  {
1092  normals[k][offset + j] = -locnormals[k][j];
1093  }
1094  }
1095  }
1096  else
1097  {
1098  // Parallel case: need to interpolate normal.
1099  Array<OneD, Array<OneD, NekDouble> > normal(coordim);
1100 
1101  for (int p = 0; p < coordim; ++p)
1102  {
1103  normal[p] = Array<OneD, NekDouble>(e_npoints,0.0);
1104  LibUtilities::PointsKey to_key =
1105  loc_exp->GetBasis(0)->GetPointsKey();
1106  LibUtilities::PointsKey from_key =
1107  loc_elmt->GetBasis(0)->GetPointsKey();
1108  LibUtilities::Interp1D(from_key,
1109  locnormals[p],
1110  to_key,
1111  normal[p]);
1112  }
1113 
1114  offset = m_phys_offset[i];
1115 
1116  // Process each point in the expansion.
1117  for (j = 0; j < e_npoints; ++j)
1118  {
1119  // Process each spatial dimension and copy the values
1120  // into the output array.
1121  for (k = 0; k < coordim; ++k)
1122  {
1123  normals[k][offset + j] = normal[k][j];
1124  }
1125  }
1126  }
1127  }
1128  else
1129  {
1130  // Get the physical data offset for this expansion.
1131  offset = m_phys_offset[i];
1132 
1133  // Process each point in the expansion.
1134  for (j = 0; j < e_npoints; ++j)
1135  {
1136  // Process each spatial dimension and copy the values
1137  // into the output array.
1138  for (k = 0; k < coordim; ++k)
1139  {
1140  normals[k][offset + j] = locnormals[k][j];
1141  }
1142  }
1143  }
1144  }
1145  }
1146 
1147  /**
1148  *
1149  */
1150  void ExpList1D::v_ReadGlobalOptimizationParameters()
1151  {
1152 // Array<OneD, int> NumShape(1,0);
1153 // NumShape[0] = GetExpSize();
1154 //
1155 // int one = 1;
1156 // m_globalOptParam = MemoryManager<NekOptimize::GlobalOptParam>
1157 // ::AllocateSharedPtr(m_session,one,NumShape);
1158  }
1159 
1160 
1161  /**
1162  * Writes out the header for a <PIECE> VTK XML segment describing the
1163  * geometric information which comprises this element. This includes
1164  * vertex coordinates for each quadrature point, vertex connectivity
1165  * information, cell types and cell offset data.
1166  *
1167  * @param outfile Output stream to write data to.
1168  */
1169  void ExpList1D::v_WriteVtkPieceHeader(std::ofstream &outfile, int expansion)
1170  {
1171  int i,j;
1172  int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1173  int ntot = nquad0;
1174  int ntotminus = (nquad0-1);
1175 
1176  Array<OneD,NekDouble> coords[3];
1177  coords[0] = Array<OneD,NekDouble>(ntot);
1178  coords[1] = Array<OneD,NekDouble>(ntot);
1179  coords[2] = Array<OneD,NekDouble>(ntot);
1180  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1181 
1182  outfile << " <Piece NumberOfPoints=\""
1183  << ntot << "\" NumberOfCells=\""
1184  << ntotminus << "\">" << endl;
1185  outfile << " <Points>" << endl;
1186  outfile << " <DataArray type=\"Float32\" "
1187  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1188  outfile << " ";
1189  for (i = 0; i < ntot; ++i)
1190  {
1191  for (j = 0; j < 3; ++j)
1192  {
1193  outfile << setprecision(8) << scientific
1194  << (float)coords[j][i] << " ";
1195  }
1196  outfile << endl;
1197  }
1198  outfile << endl;
1199  outfile << " </DataArray>" << endl;
1200  outfile << " </Points>" << endl;
1201  outfile << " <Cells>" << endl;
1202  outfile << " <DataArray type=\"Int32\" "
1203  << "Name=\"connectivity\" format=\"ascii\">" << endl;
1204  for (i = 0; i < nquad0-1; ++i)
1205  {
1206  outfile << i << " " << i+1 << endl;
1207  }
1208  outfile << endl;
1209  outfile << " </DataArray>" << endl;
1210  outfile << " <DataArray type=\"Int32\" "
1211  << "Name=\"offsets\" format=\"ascii\">" << endl;
1212  for (i = 0; i < ntotminus; ++i)
1213  {
1214  outfile << i*2+2 << " ";
1215  }
1216  outfile << endl;
1217  outfile << " </DataArray>" << endl;
1218  outfile << " <DataArray type=\"UInt8\" "
1219  << "Name=\"types\" format=\"ascii\">" << endl;
1220  for (i = 0; i < ntotminus; ++i)
1221  {
1222  outfile << "3 ";
1223  }
1224  outfile << endl;
1225  outfile << " </DataArray>" << endl;
1226  outfile << " </Cells>" << endl;
1227  outfile << " <PointData>" << endl;
1228  }
1229 
1230  } //end of namespace
1231 } //end of namespace
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2