Nektar++
ExpList.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpList.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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Expansion list definition
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iomanip>
36 
37 #include <boost/core/ignore_unused.hpp>
38 
39 #include <MultiRegions/ExpList.h>
40 #include <LibUtilities/Communication/Comm.h>
41 #include <MultiRegions/GlobalLinSys.h>
42 
43 #include <LocalRegions/PointExp.h>
44 #include <LocalRegions/SegExp.h>
45 #include <LocalRegions/TriExp.h>
46 #include <LocalRegions/QuadExp.h>
47 #include <LocalRegions/NodalTriExp.h>
48 #include <LocalRegions/HexExp.h>
49 #include <LocalRegions/PrismExp.h>
50 #include <LocalRegions/PyrExp.h>
51 #include <LocalRegions/TetExp.h>
52 #include <LocalRegions/MatrixKey.h> // for MatrixKey
53 #include <LocalRegions/Expansion3D.h>
54 #include <LibUtilities/Foundations/Interp.h>
55 
56 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h> // for AssemblyMapCG, etc
57 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h> // for AssemblyMapCG, etc
58 #include <MultiRegions/GlobalLinSysKey.h> // for GlobalLinSysKey
59 #include <MultiRegions/GlobalMatrix.h> // for GlobalMatrix, etc
60 #include <MultiRegions/GlobalMatrixKey.h> // for GlobalMatrixKey
61 
62 #include <LibUtilities/LinearAlgebra/SparseMatrixFwd.hpp>
63 #include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
64 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
65 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
66 #include <LibUtilities/Polylib/Polylib.h>
67 #include <LibUtilities/Foundations/Interp.h>
68 #include <LibUtilities/Foundations/PhysGalerkinProject.h>
69 
70 #include <Collections/CollectionOptimisation.h>
71 #include <Collections/Operator.h>
72 
73 #include <LibUtilities/BasicUtils/Timer.h>
74 #include <LibUtilities/BasicUtils/Likwid.hpp>
75 
76 using namespace std;
77 
78 namespace Nektar
79 {
80  namespace MultiRegions
81  {
82  /**
83  * @class ExpList
84  * All multi-elemental expansions \f$u^{\delta}(\boldsymbol{x})\f$ can
85  * be considered as the assembly of the various elemental contributions.
86  * On a discrete level, this yields,
87  * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
88  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
89  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
90  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
91  * As it is the lowest level class, it contains the definition of the
92  * common data and common routines to all multi-elemental expansions.
93  *
94  * The class stores a vector of expansions, \a m_exp, (each derived from
95  * StdRegions#StdExpansion) which define the constituent components of
96  * the domain. The coefficients from these expansions are concatenated
97  * in \a m_coeffs, while the expansion evaluated at the quadrature
98  * points is stored in \a m_phys.
99  */
100 
101  /**
102  * Creates an empty expansion list.
103  */
104  ExpList::ExpList(const ExpansionType type):
105  m_expType(type),
106  m_ncoeffs(0),
107  m_npoints(0),
108  m_physState(false),
109  m_exp(MemoryManager<LocalRegions::ExpansionVector>
110  ::AllocateSharedPtr()),
111  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
112  m_WaveSpace(false)
113  {
114  }
115 
116 
117  /*----------------------------------------------------------------*/
118  /* Copy Construtor */
119  /*-----------------------------------------------------------------*/
120 
121  /**
122  * Copies an existing expansion list.
123  * @param in Source expansion list.
124  */
125  ExpList::ExpList(const ExpList &in, const bool DeclareCoeffPhysArrays):
126  std::enable_shared_from_this<ExpList>(in),
127  m_expType(in.m_expType),
128  m_comm(in.m_comm),
129  m_session(in.m_session),
130  m_graph(in.m_graph),
131  m_ncoeffs(in.m_ncoeffs),
132  m_npoints(in.m_npoints),
133  m_physState(false),
134  m_exp(in.m_exp),
135  m_collections(in.m_collections),
136  m_collectionsDoInit(in.m_collectionsDoInit),
137  m_coll_coeff_offset(in.m_coll_coeff_offset),
138  m_coll_phys_offset(in.m_coll_phys_offset),
139  m_coeff_offset(in.m_coeff_offset),
140  m_phys_offset(in.m_phys_offset),
141  m_blockMat(in.m_blockMat),
142  m_WaveSpace(false)
143  {
144 
145  // Set up m_coeffs, m_phys and offset arrays.
146  // use this to keep memory declaration in one place
147  SetupCoeffPhys(DeclareCoeffPhysArrays, false);
148  }
149 
150  /**
151  * Copies the eIds elements from an existing expansion list.
152  * @param in Source expansion list.
153  * @param in elements that will be in the new exp list.
154  */
155  ExpList::ExpList(const ExpList &in,
156  const std::vector<unsigned int> &eIDs,
157  const bool DeclareCoeffPhysArrays,
158  const Collections::ImplementationType ImpType):
159  m_expType(in.m_expType),
160  m_comm(in.m_comm),
161  m_session(in.m_session),
162  m_graph(in.m_graph),
163  m_physState(false),
164  m_exp(MemoryManager<LocalRegions::ExpansionVector>
165  ::AllocateSharedPtr()),
166  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
167  m_WaveSpace(false)
168  {
169  for (int i=0; i < eIDs.size(); ++i)
170  {
171  (*m_exp).push_back( (*(in.m_exp))[eIDs[i]]);
172  }
173 
174  // Set up m_coeffs, m_phys and offset arrays.
175  SetupCoeffPhys(DeclareCoeffPhysArrays);
176 
177  // set up collections
178  CreateCollections(ImpType);
179  }
180 
181 
182  /**
183  * Given a meshgraph \a graph, containing information about
184  * the domain and the spectral/hp element expansion, this
185  * constructor fills the list of local expansions
186  * \texttt{m_exp} with the proper expansions, calculates the
187  * total number of quadrature points \f$x_i\f$ and local
188  * expansion coefficients \f$\hat{u}^e_n\f$ and allocates
189  * memory for the arrays #m_coeffs and #m_phys.
190  *
191  * @param pSession A session within information about expansion
192  *
193  * @param graph A meshgraph, containing information about the
194  * domain and the spectral/hp element expansion.
195  *
196  * @param DeclareCoeffPhysArrays Declare the coefficient and
197  * phys space arrays
198  *
199  * @param ImpType Detail about the implementation type to use
200  * in operators
201  */
202  ExpList::ExpList(const LibUtilities::SessionReaderSharedPtr &pSession,
203  const SpatialDomains::MeshGraphSharedPtr &graph,
204  const bool DeclareCoeffPhysArrays,
205  const std::string &var,
206  const Collections::ImplementationType ImpType):
207  m_comm(pSession->GetComm()),
208  m_session(pSession),
209  m_graph(graph),
210  m_physState(false),
211  m_exp(MemoryManager<LocalRegions::ExpansionVector>
212  ::AllocateSharedPtr()),
213  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
214  m_WaveSpace(false)
215  {
216  // Retrieve the list of expansions
217  const SpatialDomains::ExpansionInfoMap &expansions
218  = graph->GetExpansionInfo(var);
219 
220  // Initialise Expansionn Vector
221  InitialiseExpVector(expansions);
222 
223  // Setup phys coeff space
224  SetupCoeffPhys(DeclareCoeffPhysArrays);
225 
226  // Initialise collection
227  CreateCollections(ImpType);
228  }
229 
230  /**
231  * Given an expansion vector \a expansions, containing
232  * information about the domain and the spectral/hp element
233  * expansion, this constructor fills the list of local
234  * expansions \texttt{m_exp} with the proper expansions,
235  * calculates the total number of quadrature points
236  * \f$\boldsymbol{x}_i\f$ and local expansion coefficients
237  * \f$\hat{u}^e_n\f$ and allocates memory for the arrays
238  * #m_coeffs and #m_phys.
239  *
240  * @param pSession A session within information about expansion
241  * @param expansions A vector containing information about the
242  * domain and the spectral/hp element
243  * expansion.
244  * @param DeclareCoeffPhysArrays Declare the coefficient and
245  * phys space arrays
246  * @param ImpType Detail about the implementation type to use
247  * in operators
248  */
249  ExpList::ExpList(const LibUtilities::SessionReaderSharedPtr &pSession,
250  const SpatialDomains::ExpansionInfoMap &expansions,
251  const bool DeclareCoeffPhysArrays,
252  const Collections::ImplementationType ImpType):
253  m_comm(pSession->GetComm()),
254  m_session(pSession),
255  m_physState(false),
256  m_exp(MemoryManager<LocalRegions::ExpansionVector>
257  ::AllocateSharedPtr()),
258  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
259  m_WaveSpace(false)
260  {
261 
262  // Initialise expansion vector
263  InitialiseExpVector(expansions);
264 
265  // Set up m_coeffs, m_phys and offset arrays.
266  SetupCoeffPhys(DeclareCoeffPhysArrays);
267 
268  // Setup Collection
269  CreateCollections(ImpType);
270  }
271 
272  //----------------------------------------------------------------------
273  // 0D Expansion Constructors
274  //----------------------------------------------------------------------
275  ExpList::ExpList(const SpatialDomains::PointGeomSharedPtr &geom):
276  m_expType(e0D),
277  m_ncoeffs(1),
278  m_npoints(1),
279  m_physState(false),
280  m_exp(MemoryManager<LocalRegions::ExpansionVector>
281  ::AllocateSharedPtr()),
282  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
283  m_WaveSpace(false)
284  {
285  LocalRegions::PointExpSharedPtr Point =
286  MemoryManager<LocalRegions::PointExp>::AllocateSharedPtr(geom);
287  (*m_exp).push_back(Point);
288 
289  SetupCoeffPhys();
290  }
291 
292 
293  /**
294  * Store expansions for the trace space expansions used in
295  * DisContField
296  *
297  * @param pSession A session within information about expansion
298  * @param bndConstraint Array of ExpList1D objects each containing a
299  * 1D spectral/hp element expansion on a single
300  * boundary region.
301  * @param bndCond Array of BoundaryCondition objects which contain
302  * information about the boundary conditions on the
303  * different boundary regions.
304  * @param locexp Complete domain expansion list.
305  * @param graph mesh corresponding to the expansion list.
306  * @param DeclareCoeffPhysArrays Declare the coefficient and
307  * phys space arrays
308  * @param variable The variable name associated with the expansion
309  * @param ImpType Detail about the implementation type to use
310  * in operators
311  */
312  ExpList::ExpList(
313  const LibUtilities::SessionReaderSharedPtr &pSession,
314  const Array<OneD,const ExpListSharedPtr> &bndConstraint,
315  const Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
316  &bndCond,
317  const LocalRegions::ExpansionVector &locexp,
318  const SpatialDomains::MeshGraphSharedPtr &graph,
319  const bool DeclareCoeffPhysArrays,
320  const std::string variable,
321  const Collections::ImplementationType ImpType):
322  m_comm(pSession->GetComm()),
323  m_session(pSession),
324  m_graph(graph),
325  m_physState(false),
326  m_exp(MemoryManager<LocalRegions::ExpansionVector>
327  ::AllocateSharedPtr()),
328  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
329  m_WaveSpace(false)
330  {
331  boost::ignore_unused(variable,ImpType);
332  int i, j, id, elmtid = 0;
333  set<int> tracesDone;
334 
335  SpatialDomains::PointGeomSharedPtr PointGeom;
336  SpatialDomains::Geometry1DSharedPtr segGeom;
337  SpatialDomains::Geometry2DSharedPtr ElGeom;
338  SpatialDomains::Geometry2DSharedPtr FaceGeom;
339  SpatialDomains::QuadGeomSharedPtr QuadGeom;
340  SpatialDomains::TriGeomSharedPtr TriGeom;
341 
342  LocalRegions::ExpansionSharedPtr exp;
343  LocalRegions::Expansion0DSharedPtr exp0D;
344  LocalRegions::Expansion1DSharedPtr exp1D;
345  LocalRegions::Expansion2DSharedPtr exp2D;
346  LocalRegions::Expansion3DSharedPtr exp3D;
347 
348  // First loop over boundary conditions to reorder
349  // Dirichlet boundaries
350  for(i = 0; i < bndCond.size(); ++i)
351  {
352  if(bndCond[i]->GetBoundaryConditionType() ==
353  SpatialDomains::eDirichlet)
354  {
355  for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
356  {
357  if((exp0D = std::dynamic_pointer_cast<
358  LocalRegions::Expansion0D>(
359  bndConstraint[i]->GetExp(j))))
360  {
361  m_expType = e0D;
362 
363  PointGeom = exp0D->GetGeom()->GetVertex(0);
364  exp= MemoryManager<LocalRegions::PointExp>::
365  AllocateSharedPtr(PointGeom);
366  tracesDone.insert(PointGeom->GetVid());
367  }
368  else if((exp1D = std::dynamic_pointer_cast<
369  LocalRegions::Expansion1D>(
370  bndConstraint[i]->GetExp(j))))
371  {
372  m_expType = e1D;
373 
374  LibUtilities::BasisKey bkey = exp1D->
375  GetBasis(0)->GetBasisKey();
376  segGeom = exp1D->GetGeom1D();
377  exp = MemoryManager<LocalRegions::SegExp>
378  ::AllocateSharedPtr(bkey, segGeom);
379  tracesDone.insert(segGeom->GetGlobalID());
380 
381  }
382  else if ((exp2D = std::dynamic_pointer_cast
383  <LocalRegions::Expansion2D>(bndConstraint[i]->
384  GetExp(j))))
385  {
386  m_expType = e2D;
387 
388  LibUtilities::BasisKey bkey0 = exp2D
389  ->GetBasis(0)->GetBasisKey();
390  LibUtilities::BasisKey bkey1 = exp2D
391  ->GetBasis(1)->GetBasisKey();
392  FaceGeom = exp2D->GetGeom2D();
393 
394  //if face is a quad
395  if((QuadGeom = std::dynamic_pointer_cast<
396  SpatialDomains::QuadGeom>(FaceGeom)))
397  {
398  exp = MemoryManager<LocalRegions::QuadExp>
399  ::AllocateSharedPtr(bkey0, bkey1, QuadGeom);
400  tracesDone.insert(QuadGeom->GetGlobalID());
401  }
402  //if face is a triangle
403  else if((TriGeom = std::dynamic_pointer_cast<
404  SpatialDomains::TriGeom>(FaceGeom)))
405  {
406  exp = MemoryManager<LocalRegions::TriExp>
407  ::AllocateSharedPtr(bkey0, bkey1, TriGeom);
408  tracesDone.insert(TriGeom->GetGlobalID());
409  }
410  else
411  {
412  NEKERROR(ErrorUtil::efatal,"dynamic cast to a "
413  "proper face geometry failed");
414  }
415  }
416  // Assign next id
417  exp->SetElmtId(elmtid++);
418 
419  // Add the expansion
420  (*m_exp).push_back(exp);
421  }
422  }
423  }
424 
425  map<int, pair<SpatialDomains::Geometry1DSharedPtr,
426  LibUtilities::BasisKey> > edgeOrders;
427 
428  map<int, pair<SpatialDomains::Geometry2DSharedPtr,
429  pair<LibUtilities::BasisKey,
430  LibUtilities::BasisKey> > > faceOrders;
431 
432  for(i = 0; i < locexp.size(); ++i)
433  {
434  if((exp1D =
435  std::dynamic_pointer_cast<
436  LocalRegions::Expansion1D>(locexp[i])))
437  {
438  m_expType = e0D;
439 
440  for(j = 0; j < 2; ++j)
441  {
442  PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
443  id = PointGeom->GetVid();
444 
445  // Ignore Dirichlet edges
446  if (tracesDone.count(id) != 0)
447  {
448  continue;
449  }
450 
451  exp = MemoryManager<LocalRegions::PointExp>::
452  AllocateSharedPtr(PointGeom);
453  tracesDone.insert(id);
454  exp->SetElmtId(elmtid++);
455  (*m_exp).push_back(exp);
456  }
457  }
458  else if((exp2D =
459  std::dynamic_pointer_cast<
460  LocalRegions::Expansion2D>(locexp[i])))
461  {
462  m_expType = e1D;
463  for(j = 0; j < locexp[i]->GetNtraces(); ++j)
464  {
465  segGeom = exp2D->GetGeom2D()->GetEdge(j);
466  id = segGeom->GetGlobalID();
467  // Ignore Dirichlet edges
468  if (tracesDone.count(id) != 0)
469  {
470  continue;
471  }
472 
473  auto it = edgeOrders.find(id);
474 
475  if (it == edgeOrders.end())
476  {
477  edgeOrders.insert(std::make_pair(id, std::make_pair(
478  segGeom, locexp[i]->GetTraceBasisKey(j))));
479  }
480  else // variable modes/points
481  {
482  LibUtilities::BasisKey edge
483  = locexp[i]->GetTraceBasisKey(j);
484  LibUtilities::BasisKey existing
485  = it->second.second;
486 
487  int np1 = edge .GetNumPoints();
488  int np2 = existing.GetNumPoints();
489  int nm1 = edge .GetNumModes ();
490  int nm2 = existing.GetNumModes ();
491 
492  if (np2 >= np1 && nm2 >= nm1)
493  {
494  continue;
495  }
496  else if (np2 < np1 && nm2 < nm1)
497  {
498  it->second.second = edge;
499  }
500  else
501  {
502  NEKERROR(ErrorUtil::efatal,
503  "inappropriate number of points/modes (max"
504  "num of points is not set with max order)");
505  }
506  }
507  }
508  }
509  else if((exp3D =
510  dynamic_pointer_cast<
511  LocalRegions::Expansion3D>(locexp[i])))
512  {
513  m_expType = e2D;
514  for (j = 0; j < exp3D->GetNtraces(); ++j)
515  {
516  FaceGeom = exp3D->GetGeom3D()->GetFace(j);
517  id = FaceGeom->GetGlobalID();
518 
519  if(tracesDone.count(id) != 0)
520  {
521  continue;
522  }
523  auto it = faceOrders.find(id);
524 
525  if (it == faceOrders.end())
526  {
527  LibUtilities::BasisKey face_dir0
528  = locexp[i]->GetTraceBasisKey(j,0);
529  LibUtilities::BasisKey face_dir1
530  = locexp[i]->GetTraceBasisKey(j,1);
531 
532  faceOrders.insert(
533  std::make_pair(
534  id, std::make_pair(FaceGeom,
535  std::make_pair(face_dir0, face_dir1))));
536  }
537  else // variable modes/points
538  {
539  LibUtilities::BasisKey face0 =
540  locexp[i]->GetTraceBasisKey(j,0);
541  LibUtilities::BasisKey face1 =
542  locexp[i]->GetTraceBasisKey(j,1);
543  LibUtilities::BasisKey existing0 =
544  it->second.second.first;
545  LibUtilities::BasisKey existing1 =
546  it->second.second.second;
547 
548  int np11 = face0 .GetNumPoints();
549  int np12 = face1 .GetNumPoints();
550  int np21 = existing0.GetNumPoints();
551  int np22 = existing1.GetNumPoints();
552  int nm11 = face0 .GetNumModes ();
553  int nm12 = face1 .GetNumModes ();
554  int nm21 = existing0.GetNumModes ();
555  int nm22 = existing1.GetNumModes ();
556 
557  if ((np22 >= np12 || np21 >= np11) &&
558  (nm22 >= nm12 || nm21 >= nm11))
559  {
560  continue;
561  }
562  else if((np22 < np12 || np21 < np11) &&
563  (nm22 < nm12 || nm21 < nm11))
564  {
565  it->second.second.first = face0;
566  it->second.second.second = face1;
567  }
568  else
569  {
570  NEKERROR(ErrorUtil::efatal,
571  "inappropriate number of points/modes (max"
572  "num of points is not set with max order)");
573  }
574  }
575  }
576  }
577  }
578 
579  int nproc = m_comm->GetSize(); // number of processors
580  int tracepr = m_comm->GetRank(); // ID processor
581 
582  if (nproc > 1)
583  {
584  int tCnt = 0;
585 
586  // Count the number of traces on each partition
587  for(i = 0; i < locexp.size(); ++i)
588  {
589  tCnt += locexp[i]->GetNtraces();
590  }
591 
592  // Set up the offset and the array that will contain the list of
593  // edge IDs, then reduce this across processors.
594  Array<OneD, int> tracesCnt(nproc, 0);
595  tracesCnt[tracepr] = tCnt;
596  m_comm->AllReduce(tracesCnt, LibUtilities::ReduceSum);
597 
598  // Set up offset array.
599  int totTraceCnt = Vmath::Vsum(nproc, tracesCnt, 1);
600  Array<OneD, int> tTotOffsets(nproc,0);
601 
602  for (i = 1; i < nproc; ++i)
603  {
604  tTotOffsets[i] = tTotOffsets[i-1] + tracesCnt[i-1];
605  }
606 
607  // Local list of the edges per element
608  Array<OneD, int> TracesTotID(totTraceCnt, 0);
609  Array<OneD, int> TracesTotNm0(totTraceCnt, 0);
610  Array<OneD, int> TracesTotNm1(totTraceCnt, 0);
611  Array<OneD, int> TracesTotPnts0(totTraceCnt, 0);
612  Array<OneD, int> TracesTotPnts1(totTraceCnt, 0);
613 
614  int cntr = tTotOffsets[tracepr];
615 
616  for(i = 0; i < locexp.size(); ++i)
617  {
618  if((exp2D = locexp[i]->as<LocalRegions::Expansion2D>()))
619  {
620 
621  int nedges = locexp[i]->GetNtraces();
622 
623  for(j = 0; j < nedges; ++j, ++cntr)
624  {
625  LibUtilities::BasisKey bkeyEdge =
626  locexp[i]->GetTraceBasisKey(j);
627  TracesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
628  TracesTotNm0 [cntr] = bkeyEdge.GetNumModes();
629  TracesTotPnts0[cntr] = bkeyEdge.GetNumPoints();
630  }
631  }
632  else if((exp3D = locexp[i]->as<LocalRegions::Expansion3D>()))
633  {
634  int nfaces = locexp[i]->GetNtraces();
635 
636  for(j = 0; j < nfaces; ++j, ++cntr)
637  {
638  LibUtilities::BasisKey face_dir0
639  = locexp[i]->GetTraceBasisKey(j,0);
640  LibUtilities::BasisKey face_dir1
641  = locexp[i]->GetTraceBasisKey(j,1);
642 
643  TracesTotID[cntr] = exp3D->GetGeom3D()->GetFid(j);
644  TracesTotNm0[cntr] = face_dir0.GetNumModes ();
645  TracesTotNm1[cntr] = face_dir1.GetNumModes ();
646  TracesTotPnts0[cntr] = face_dir0.GetNumPoints();
647  TracesTotPnts1[cntr] = face_dir1.GetNumPoints();
648  }
649  }
650  }
651 
652  m_comm->AllReduce(TracesTotID, LibUtilities::ReduceSum);
653  m_comm->AllReduce(TracesTotNm0, LibUtilities::ReduceSum);
654  m_comm->AllReduce(TracesTotPnts0, LibUtilities::ReduceSum);
655  if(m_expType == e2D)
656  {
657  m_comm->AllReduce(TracesTotNm1, LibUtilities::ReduceSum);
658  m_comm->AllReduce(TracesTotPnts1, LibUtilities::ReduceSum);
659  }
660 
661  if(edgeOrders.size())
662  {
663  for (i = 0; i < totTraceCnt; ++i)
664  {
665  auto it = edgeOrders.find(TracesTotID[i]);
666 
667  if (it == edgeOrders.end())
668  {
669  continue;
670  }
671 
672  LibUtilities::BasisKey existing
673  = it->second.second;
674  LibUtilities::BasisKey edge(existing.GetBasisType(),
675  TracesTotNm0[i],
676  LibUtilities::PointsKey(
677  TracesTotPnts0[i],
678  existing.GetPointsType()));
679 
680  int np1 = edge .GetNumPoints();
681  int np2 = existing.GetNumPoints();
682  int nm1 = edge .GetNumModes ();
683  int nm2 = existing.GetNumModes ();
684 
685  if (np2 >= np1 && nm2 >= nm1)
686  {
687  continue;
688  }
689  else if (np2 < np1 && nm2 < nm1)
690  {
691  it->second.second = edge;
692  }
693  else
694  {
695  NEKERROR(ErrorUtil::efatal,
696  "inappropriate number of points/modes (max "
697  "num of points is not set with max order)");
698  }
699  }
700  }
701  else if(faceOrders.size())
702  {
703  for (i = 0; i < totTraceCnt; ++i)
704  {
705  auto it = faceOrders.find(TracesTotID[i]);
706 
707  if (it == faceOrders.end())
708  {
709  continue;
710  }
711 
712  LibUtilities::BasisKey existing0 =
713  it->second.second.first;
714  LibUtilities::BasisKey existing1 =
715  it->second.second.second;
716  LibUtilities::BasisKey face0(
717  existing0.GetBasisType(), TracesTotNm0[i],
718  LibUtilities::PointsKey(TracesTotPnts0[i],
719  existing0.GetPointsType()));
720  LibUtilities::BasisKey face1(
721  existing1.GetBasisType(), TracesTotNm1[i],
722  LibUtilities::PointsKey(TracesTotPnts1[i],
723  existing1.GetPointsType()));
724 
725  int np11 = face0 .GetNumPoints();
726  int np12 = face1 .GetNumPoints();
727  int np21 = existing0.GetNumPoints();
728  int np22 = existing1.GetNumPoints();
729  int nm11 = face0 .GetNumModes ();
730  int nm12 = face1 .GetNumModes ();
731  int nm21 = existing0.GetNumModes ();
732  int nm22 = existing1.GetNumModes ();
733 
734  if ((np22 >= np12 || np21 >= np11) &&
735  (nm22 >= nm12 || nm21 >= nm11))
736  {
737  continue;
738  }
739  else if((np22 < np12 || np21 < np11) &&
740  (nm22 < nm12 || nm21 < nm11))
741  {
742  it->second.second.first = face0;
743  it->second.second.second = face1;
744  }
745  else
746  {
747  NEKERROR(ErrorUtil::efatal,
748  "inappropriate number of points/modes (max "
749  "num of points is not set with max order)");
750  }
751  }
752  }
753  }
754 
755  if(edgeOrders.size())
756  {
757  for (auto &it : edgeOrders)
758  {
759  exp = MemoryManager<LocalRegions::SegExp>
760  ::AllocateSharedPtr(it.second.second, it.second.first);
761  exp->SetElmtId(elmtid++);
762  (*m_exp).push_back(exp);
763  }
764  }
765  else
766  {
767  for (auto &it : faceOrders)
768  {
769  FaceGeom = it.second.first;
770 
771  if ((QuadGeom = std::dynamic_pointer_cast<
772  SpatialDomains::QuadGeom>(FaceGeom)))
773  {
774  exp = MemoryManager<LocalRegions::QuadExp>
775  ::AllocateSharedPtr(it.second.second.first,
776  it.second.second.second,
777  QuadGeom);
778  }
779  else if ((TriGeom = std::dynamic_pointer_cast<
780  SpatialDomains::TriGeom>(FaceGeom)))
781  {
782  exp = MemoryManager<LocalRegions::TriExp>
783  ::AllocateSharedPtr(it.second.second.first,
784  it.second.second.second,
785  TriGeom);
786  }
787  exp->SetElmtId(elmtid++);
788  (*m_exp).push_back(exp);
789  }
790  }
791 
792  // Set up m_coeffs, m_phys and offset arrays.
793  SetupCoeffPhys(DeclareCoeffPhysArrays);
794 
795 
796  // Set up collections
797  if(m_expType != e0D)
798  {
799  CreateCollections(ImpType);
800  }
801  }
802 
803  /**
804  * Fills the list of local expansions with the trace from the mesh
805  * specified by \a domain. This CompositeMap contains a list of
806  * Composites which define the boundary. It is also used to set up
807  * expansion domains in the 1D Pulse Wave solver.
808  *
809  * @param pSession A session within information about expansion
810  * @param domain A domain, comprising of one or more composite
811  * regions,
812  * @param graph A mesh, containing information about the
813  * domain and the spectral/hp element expansion.
814  * @param DeclareCoeffPhysArrays Declare the coefficient and
815  * phys space arrays. Default is true.
816  * @param variable The variable name associated with the expansion
817  * @param SetToOneSpaceDimension Reduce to one space dimension expansion
818  * @param comm An optional communicator that can be used with the
819  * boundary expansion in case of more global
820  * parallel operations. Default to a Null Communicator
821  * @param ImpType Detail about the implementation type to use
822  * in operators. Default is eNoImpType.
823  *
824  */
825  ExpList::ExpList(const LibUtilities::SessionReaderSharedPtr &pSession,
826  const SpatialDomains::CompositeMap &domain,
827  const SpatialDomains::MeshGraphSharedPtr &graph,
828  const bool DeclareCoeffPhysArrays,
829  const std::string variable,
830  bool SetToOneSpaceDimension,
831  const LibUtilities::CommSharedPtr comm,
832  const Collections::ImplementationType ImpType):
833  m_comm(comm),
834  m_session(pSession),
835  m_graph(graph),
836  m_physState(false),
837  m_exp(MemoryManager<LocalRegions::ExpansionVector>
838  ::AllocateSharedPtr()),
839  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
840  m_WaveSpace(false)
841  {
842  int j, elmtid=0;
843  SpatialDomains::PointGeomSharedPtr PtGeom;
844  SpatialDomains::SegGeomSharedPtr SegGeom;
845  SpatialDomains::TriGeomSharedPtr TriGeom;
846  SpatialDomains::QuadGeomSharedPtr QuadGeom;
847 
848  LocalRegions::ExpansionSharedPtr exp;
849 
850  LibUtilities::PointsType TriNb;
851 
852  int meshdim = graph->GetMeshDimension();
853 
854  // Retrieve the list of expansions (needed of meshdim == 1
855  const SpatialDomains::ExpansionInfoMap &expansions
856  = graph->GetExpansionInfo(variable);
857 
858  // Retrieve the list of expansions
859  // Process each composite region.
860  for(auto &compIt : domain)
861  {
862  // Process each expansion in the region.
863  for(j = 0; j < compIt.second->m_geomVec.size(); ++j)
864  {
865  if((PtGeom = std::dynamic_pointer_cast <
866  SpatialDomains::PointGeom>(compIt.second->m_geomVec[j])))
867  {
868  m_expType = e0D;
869 
870  exp = MemoryManager<LocalRegions::PointExp>
871  ::AllocateSharedPtr(PtGeom);
872  }
873  else if((SegGeom = std::dynamic_pointer_cast<
874  SpatialDomains::SegGeom>(compIt.second->m_geomVec[j])))
875  {
876  m_expType = e1D;
877 
878  // Retrieve the basis key from the expansion.
879  LibUtilities::BasisKey bkey = LibUtilities::NullBasisKey;
880 
881  if(meshdim == 1)
882  {
883  auto expIt = expansions.find(SegGeom->GetGlobalID());
884  ASSERTL0(expIt != expansions.end(),
885  "Failed to find basis key");
886  bkey = expIt->second->m_basisKeyVector[0];
887  }
888  else
889  {
890  bkey = graph->GetEdgeBasisKey(SegGeom, variable);
891  }
892 
893  if(SetToOneSpaceDimension)
894  {
895  SpatialDomains::SegGeomSharedPtr OneDSegmentGeom =
896  SegGeom->GenerateOneSpaceDimGeom();
897 
898  exp = MemoryManager<LocalRegions::SegExp>
899  ::AllocateSharedPtr(bkey, OneDSegmentGeom);
900  }
901  else
902  {
903 
904  exp = MemoryManager<LocalRegions::SegExp>
905  ::AllocateSharedPtr(bkey, SegGeom);
906  }
907  }
908  else if ((TriGeom = std::dynamic_pointer_cast<
909  SpatialDomains::TriGeom>(compIt.second->m_geomVec[j])))
910  {
911  m_expType = e2D;
912 
913  LibUtilities::BasisKey TriBa
914  = graph->GetFaceBasisKey(TriGeom,0,variable);
915  LibUtilities::BasisKey TriBb
916  = graph->GetFaceBasisKey(TriGeom,1,variable);
917 
918  if (graph->GetExpansionInfo().begin()->second->
919  m_basisKeyVector[0].GetBasisType() ==
920  LibUtilities::eGLL_Lagrange)
921  {
922  NEKERROR(ErrorUtil::efatal,"This method needs sorting");
923  TriNb = LibUtilities::eNodalTriElec;
924 
925  exp = MemoryManager<LocalRegions::NodalTriExp>
926  ::AllocateSharedPtr(TriBa,TriBb,TriNb,
927  TriGeom);
928  }
929  else
930  {
931  exp = MemoryManager<LocalRegions::TriExp>
932  ::AllocateSharedPtr(TriBa, TriBb, TriGeom);
933  }
934  }
935  else if ((QuadGeom = std::dynamic_pointer_cast<
936  SpatialDomains::QuadGeom>(compIt.second->m_geomVec[j])))
937  {
938  m_expType = e2D;
939 
940  LibUtilities::BasisKey QuadBa
941  = graph->GetFaceBasisKey(QuadGeom, 0, variable);
942  LibUtilities::BasisKey QuadBb
943  = graph->GetFaceBasisKey(QuadGeom, 1, variable);
944 
945  exp = MemoryManager<LocalRegions::QuadExp>
946  ::AllocateSharedPtr(QuadBa, QuadBb, QuadGeom);
947  }
948  else
949  {
950  NEKERROR(ErrorUtil::efatal,
951  "dynamic cast to a Geom (possibly 3D) failed");
952  }
953 
954  exp->SetElmtId(elmtid++);
955  (*m_exp).push_back(exp);
956  }
957  }
958 
959  // Set up m_coeffs, m_phys and offset arrays.
960  SetupCoeffPhys(DeclareCoeffPhysArrays);
961 
962  if(m_expType != e0D)
963  {
964  CreateCollections(ImpType);
965  }
966  }
967 
968  /**
969  * Each expansion (local element) is processed in turn to
970  * determine the number of coefficients and physical data
971  * points it contributes to the domain. Twoe arrays,
972  * #m_coeff_offset are #m_phys_offset are also initialised and
973  * updated to store the data offsets of each element in the
974  * #m_coeffs and #m_phys arrays, and the element id that each
975  * consecutive block is associated respectively.
976  * Finally we initialise #m_coeffs and #m_phys
977  */
978  void ExpList::SetupCoeffPhys(bool DeclareCoeffPhysArrays,
979  bool SetupOffsets)
980  {
981  if(SetupOffsets)
982  {
983  int i;
984 
985  // Set up offset information and array sizes
986  m_coeff_offset = Array<OneD,int>(m_exp->size());
987  m_phys_offset = Array<OneD,int>(m_exp->size());
988 
989  m_ncoeffs = m_npoints = 0;
990 
991  for(i = 0; i < m_exp->size(); ++i)
992  {
993  m_coeff_offset[i] = m_ncoeffs;
994  m_phys_offset [i] = m_npoints;
995  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
996  m_npoints += (*m_exp)[i]->GetTotPoints();
997  }
998  }
999 
1000  if(DeclareCoeffPhysArrays)
1001  {
1002  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs, 0.0);
1003  m_phys = Array<OneD, NekDouble>(m_npoints, 0.0);
1004  }
1005 
1006  m_coeffsToElmt = Array<OneD, pair<int, int> > {size_t(m_ncoeffs)};
1007 
1008  for (int i = 0; i < m_exp->size(); ++i)
1009  {
1010  int coeffs_offset = m_coeff_offset[i];
1011 
1012  int loccoeffs = (*m_exp)[i]->GetNcoeffs();
1013 
1014  for(int j = 0; j < loccoeffs; ++j)
1015  {
1016  m_coeffsToElmt[coeffs_offset+j].first = i;
1017  m_coeffsToElmt[coeffs_offset+j].second = j;
1018  }
1019  }
1020  }
1021 
1022  void ExpList::InitialiseExpVector( const
1023  SpatialDomains::ExpansionInfoMap &expmap)
1024  {
1025 
1026  SpatialDomains::SegGeomSharedPtr SegmentGeom;
1027  SpatialDomains::TriGeomSharedPtr TriangleGeom;
1028  SpatialDomains::QuadGeomSharedPtr QuadrilateralGeom;
1029  SpatialDomains::TetGeomSharedPtr TetGeom;
1030  SpatialDomains::HexGeomSharedPtr HexGeom;
1031  SpatialDomains::PrismGeomSharedPtr PrismGeom;
1032  SpatialDomains::PyrGeomSharedPtr PyrGeom;
1033 
1034  int id=0;
1035  LocalRegions::ExpansionSharedPtr exp;
1036 
1037  m_expType = eNoType;
1038  // Process each expansion in the graph
1039  for (auto &expIt : expmap)
1040  {
1041  const SpatialDomains::ExpansionInfoShPtr expInfo = expIt.second;
1042 
1043  switch(expInfo->m_basisKeyVector.size())
1044  {
1045  case 1: // Segment Expansions
1046  {
1047  ASSERTL1(m_expType == e1D || m_expType == eNoType,
1048  "Cannot mix expansion dimensions in one vector");
1049  m_expType = e1D;
1050 
1051  if ((SegmentGeom = std::dynamic_pointer_cast<
1052  SpatialDomains::SegGeom>(expInfo->m_geomShPtr)))
1053  {
1054  // Retrieve basis key from expansion
1055  LibUtilities::BasisKey bkey=
1056  expInfo->m_basisKeyVector[0];
1057 
1058  exp = MemoryManager<LocalRegions::SegExp>
1059  ::AllocateSharedPtr(bkey, SegmentGeom);
1060  }
1061  else
1062  {
1063  NEKERROR(ErrorUtil::efatal,
1064  "dynamic cast to a 1D Geom failed");
1065  }
1066  }
1067  break;
1068  case 2:
1069  {
1070  ASSERTL1(m_expType == e2D || m_expType == eNoType,
1071  "Cannot mix expansion dimensions in one vector");
1072  m_expType = e2D;
1073 
1074  LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1075  LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1076 
1077  if ((TriangleGeom = std::dynamic_pointer_cast<SpatialDomains
1078  ::TriGeom>(expInfo->m_geomShPtr)))
1079  {
1080 
1081  // This is not elegantly implemented needs re-thinking.
1082  if (Ba.GetBasisType() == LibUtilities::eGLL_Lagrange)
1083  {
1084  LibUtilities::BasisKey newBa(LibUtilities::eOrtho_A,
1085  Ba.GetNumModes(),
1086  Ba.GetPointsKey());
1087 
1088  LibUtilities::PointsType TriNb
1089  = LibUtilities::eNodalTriElec;
1090  exp = MemoryManager<LocalRegions::NodalTriExp>
1091  ::AllocateSharedPtr(newBa,Bb,TriNb,TriangleGeom);
1092  }
1093  else
1094  {
1095  exp = MemoryManager<LocalRegions::TriExp>
1096  ::AllocateSharedPtr(Ba,Bb,TriangleGeom);
1097  }
1098  }
1099  else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
1100  SpatialDomains::QuadGeom>(expInfo->m_geomShPtr)))
1101  {
1102  exp = MemoryManager<LocalRegions::QuadExp>
1103  ::AllocateSharedPtr(Ba,Bb,QuadrilateralGeom);
1104  }
1105  else
1106  {
1107  NEKERROR(ErrorUtil::efatal,
1108  "dynamic cast to a 2D Geom failed");
1109  }
1110  }
1111  break;
1112  case 3:
1113  {
1114  ASSERTL1(m_expType == e3D || m_expType == eNoType,
1115  "Cannot mix expansion dimensions in one vector");
1116  m_expType = e3D;
1117 
1118  LibUtilities::BasisKey Ba = expInfo->m_basisKeyVector[0];
1119  LibUtilities::BasisKey Bb = expInfo->m_basisKeyVector[1];
1120  LibUtilities::BasisKey Bc = expInfo->m_basisKeyVector[2];
1121 
1122  if((TetGeom = std::dynamic_pointer_cast<
1123  SpatialDomains::TetGeom>(expInfo->m_geomShPtr)))
1124  {
1125 
1126  if(Ba.GetBasisType() == LibUtilities::eGLL_Lagrange ||
1127  Ba.GetBasisType() == LibUtilities::eGauss_Lagrange)
1128  {
1129  NEKERROR(ErrorUtil::efatal,
1130  "LocalRegions::NodalTetExp is not implemented "
1131  "yet");
1132  }
1133  else
1134  {
1135  exp = MemoryManager<LocalRegions::TetExp>
1136  ::AllocateSharedPtr(Ba,Bb,Bc,TetGeom);
1137  }
1138  }
1139  else if((PrismGeom =
1140  std::dynamic_pointer_cast<SpatialDomains
1141  ::PrismGeom>(expInfo->m_geomShPtr)))
1142  {
1143  exp = MemoryManager<LocalRegions::PrismExp>
1144  ::AllocateSharedPtr(Ba,Bb,Bc,PrismGeom);
1145  }
1146  else if((PyrGeom = std::dynamic_pointer_cast<
1147  SpatialDomains::PyrGeom>(expInfo->m_geomShPtr)))
1148  {
1149 
1150  exp = MemoryManager<LocalRegions::PyrExp>
1151  ::AllocateSharedPtr(Ba,Bb,Bc,PyrGeom);
1152  }
1153  else if((HexGeom = std::dynamic_pointer_cast<
1154  SpatialDomains::HexGeom>(expInfo->m_geomShPtr)))
1155  {
1156  exp = MemoryManager<LocalRegions::HexExp>
1157  ::AllocateSharedPtr(Ba,Bb,Bc,HexGeom);
1158  }
1159  else
1160  {
1161  NEKERROR(ErrorUtil::efatal,
1162  "dynamic cast to a Geom failed");
1163  }
1164  }
1165  break;
1166  default:
1167  NEKERROR(ErrorUtil::efatal,
1168  "Dimension of basis key is greater than 3");
1169  }
1170 
1171  // Assign next id
1172  exp->SetElmtId(id++);
1173 
1174  // Add the expansion
1175  (*m_exp).push_back(exp);
1176  }
1177  }
1178 
1179  /**
1180  *
1181  */
1182  ExpansionType ExpList::GetExpType(void)
1183  {
1184  return m_expType;
1185  }
1186 
1187  ExpList::~ExpList()
1188  {
1189  }
1190 
1191  /**
1192  * Retrieves the block matrix specified by \a bkey, and computes
1193  * \f$ y=Mx \f$.
1194  * @param gkey GlobalMatrixKey specifying the block matrix to
1195  * use in the matrix-vector multiply.
1196  * @param inarray Input vector \f$ x \f$.
1197  * @param outarray Output vector \f$ y \f$.
1198  */
1199  void ExpList::MultiplyByBlockMatrix(
1200  const GlobalMatrixKey &gkey,
1201  const Array<OneD,const NekDouble> &inarray,
1202  Array<OneD, NekDouble> &outarray)
1203  {
1204  // Retrieve the block matrix using the given key.
1205  const DNekScalBlkMatSharedPtr& blockmat = GetBlockMatrix(gkey);
1206  int nrows = blockmat->GetRows();
1207  int ncols = blockmat->GetColumns();
1208 
1209  // Create NekVectors from the given data arrays
1210  NekVector<NekDouble> in (ncols,inarray, eWrapper);
1211  NekVector< NekDouble> out(nrows,outarray,eWrapper);
1212 
1213  // Perform matrix-vector multiply.
1214  out = (*blockmat)*in;
1215  }
1216 
1217  /**
1218  * multiply the metric jacobi and quadrature weights
1219  */
1220  void ExpList::MultiplyByQuadratureMetric(
1221  const Array<OneD, const NekDouble> &inarray,
1222  Array<OneD, NekDouble> &outarray)
1223  {
1224  Array<OneD,NekDouble> e_outarray;
1225 
1226  for (int i = 0; i < (*m_exp).size(); ++i)
1227  {
1228  (*m_exp)[i]->MultiplyByQuadratureMetric(
1229  inarray+m_phys_offset[i],
1230  e_outarray = outarray + m_phys_offset[i]);
1231  }
1232  }
1233 
1234  /**
1235  * Divided by the metric jacobi and quadrature weights
1236  */
1237  void ExpList::DivideByQuadratureMetric(
1238  const Array<OneD, const NekDouble> &inarray,
1239  Array<OneD, NekDouble> &outarray)
1240  {
1241  Array<OneD,NekDouble> e_outarray;
1242 
1243  for (int i = 0; i < (*m_exp).size(); ++i)
1244  {
1245  (*m_exp)[i]->DivideByQuadratureMetric(
1246  inarray+m_phys_offset[i],
1247  e_outarray = outarray + m_phys_offset[i]);
1248  }
1249  }
1250 
1251  /**
1252  * The operation is evaluated locally for every element by the function
1253  * StdRegions#StdExpansion#IProductWRTBase.
1254  *
1255  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1256  * containing the values of the function
1257  * \f$f(\boldsymbol{x})\f$ at the quadrature
1258  * points \f$\boldsymbol{x}_i\f$.
1259  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1260  * used to store the result.
1261  */
1262  void ExpList::v_IProductWRTBase_IterPerExp(
1263  const Array<OneD, const NekDouble> &inarray,
1264  Array<OneD, NekDouble> &outarray)
1265  {
1266  // initialise if required
1267  if(m_collectionsDoInit[Collections::eIProductWRTBase])
1268  {
1269  for (int i = 0; i < m_collections.size(); ++i)
1270  {
1271  m_collections[i].Initialise(Collections::eIProductWRTBase);
1272  }
1273  m_collectionsDoInit[Collections::eIProductWRTBase] = false;
1274  }
1275 
1276  Array<OneD,NekDouble> tmp;
1277  for (int i = 0; i < m_collections.size(); ++i)
1278  {
1279 
1280  m_collections[i].ApplyOperator(Collections::eIProductWRTBase,
1281  inarray + m_coll_phys_offset[i],
1282  tmp = outarray + m_coll_coeff_offset[i]);
1283  }
1284  }
1285 
1286  /**
1287  * The operation is evaluated locally for every element by the function
1288  * StdRegions#StdExpansion#IProductWRTDerivBase.
1289  *
1290  * @param dir {0,1} is the direction in which the
1291  * derivative of the basis should be taken
1292  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1293  * containing the values of the function
1294  * \f$f(\boldsymbol{x})\f$ at the quadrature
1295  * points \f$\boldsymbol{x}_i\f$.
1296  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1297  * used to store the result.
1298  */
1299  void ExpList::IProductWRTDerivBase(const int dir,
1300  const Array<OneD, const NekDouble> &inarray,
1301  Array<OneD, NekDouble> &outarray)
1302  {
1303  int i;
1304 
1305  Array<OneD,NekDouble> e_outarray;
1306 
1307  for(i = 0; i < (*m_exp).size(); ++i)
1308  {
1309  (*m_exp)[i]->IProductWRTDerivBase(dir,inarray+m_phys_offset[i],
1310  e_outarray = outarray+m_coeff_offset[i]);
1311  }
1312  }
1313 
1314 
1315  /**
1316  * @brief Directional derivative along a given direction
1317  *
1318  */
1319  void ExpList::IProductWRTDirectionalDerivBase(
1320  const Array<OneD, const NekDouble> &direction,
1321  const Array<OneD, const NekDouble> &inarray,
1322  Array<OneD, NekDouble> &outarray)
1323  {
1324  int npts_e;
1325  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1326  int nq = direction.size()/coordim;
1327 
1328  Array<OneD, NekDouble> e_outarray;
1329  Array<OneD, NekDouble> e_MFdiv;
1330 
1331  Array<OneD, NekDouble> locdir;
1332 
1333  for(int i = 0; i < (*m_exp).size(); ++i)
1334  {
1335  npts_e = (*m_exp)[i]->GetTotPoints();
1336  locdir = Array<OneD, NekDouble>(npts_e*coordim);
1337 
1338  for (int k = 0; k<coordim; ++k)
1339  {
1340  Vmath::Vcopy(npts_e, &direction[k*nq+m_phys_offset[i]], 1,
1341  &locdir[k*npts_e], 1);
1342  }
1343 
1344  (*m_exp)[i]->IProductWRTDirectionalDerivBase(
1345  locdir,
1346  inarray+m_phys_offset[i],
1347  e_outarray = outarray+m_coeff_offset[i] );
1348  }
1349  }
1350 
1351 
1352  /**
1353  * The operation is evaluated locally for every element by the function
1354  * StdRegions#StdExpansion#IProductWRTDerivBase.
1355  *
1356  * @param inarray An array of arrays of size \f$Q_{\mathrm{tot}}\f$
1357  * containing the values of the function
1358  * \f$f(\boldsymbol{x})\f$ at the quadrature
1359  * points \f$\boldsymbol{x}_i\f$ in dir directions.
1360  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
1361  * used to store the result.
1362  */
1363  void ExpList::IProductWRTDerivBase(const Array<OneD, const Array<OneD, NekDouble> > &inarray,
1364  Array<OneD, NekDouble> &outarray)
1365  {
1366  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1367  // assume coord dimension defines the size of Deriv Base
1368  int dim = GetCoordim(0);
1369 
1370  ASSERTL1(inarray.size() >= dim,"inarray is not of sufficient dimension");
1371 
1372  // initialise if required
1373  if(m_collectionsDoInit[Collections::eIProductWRTDerivBase])
1374  {
1375  for (int i = 0; i < m_collections.size(); ++i)
1376  {
1377  m_collections[i].Initialise(Collections::eIProductWRTDerivBase);
1378  }
1379  m_collectionsDoInit[Collections::eIProductWRTDerivBase] = false;
1380  }
1381 
1382  LibUtilities::Timer timer;
1383 
1384  LIKWID_MARKER_START("IProductWRTDerivBase_coll");
1385  timer.Start();
1386 
1387  switch(dim)
1388  {
1389  case 1:
1390  for (int i = 0; i < m_collections.size(); ++i)
1391  {
1392  m_collections[i].ApplyOperator(
1393  Collections::eIProductWRTDerivBase,
1394  inarray[0] + m_coll_phys_offset[i],
1395  tmp0 = outarray + m_coll_coeff_offset[i]);
1396  }
1397  break;
1398  case 2:
1399  for (int i = 0; i < m_collections.size(); ++i)
1400  {
1401  m_collections[i].ApplyOperator(
1402  Collections::eIProductWRTDerivBase,
1403  inarray[0] + m_coll_phys_offset[i],
1404  tmp0 = inarray[1] + m_coll_phys_offset[i],
1405  tmp1 = outarray + m_coll_coeff_offset[i]);
1406  }
1407  break;
1408  case 3:
1409  for (int i = 0; i < m_collections.size(); ++i)
1410  {
1411  m_collections[i].ApplyOperator(
1412  Collections::eIProductWRTDerivBase,
1413  inarray[0] + m_coll_phys_offset[i],
1414  tmp0 = inarray[1] + m_coll_phys_offset[i],
1415  tmp1 = inarray[2] + m_coll_phys_offset[i],
1416  tmp2 = outarray + m_coll_coeff_offset[i]);
1417  }
1418  break;
1419  default:
1420  NEKERROR(ErrorUtil::efatal,"Dimension of inarray not correct");
1421  break;
1422  }
1423 
1424  timer.Stop();
1425  LIKWID_MARKER_STOP("IProductWRTDerivBase_coll");
1426 
1427  // Elapsed time
1428  timer.AccumulateRegion("Collections:IProductWRTDerivBase");
1429 
1430  }
1431  /**
1432  * Given a function \f$f(\boldsymbol{x})\f$ evaluated at
1433  * the quadrature points, this function calculates the
1434  * derivatives \f$\frac{d}{dx_1}\f$, \f$\frac{d}{dx_2}\f$
1435  * and \f$\frac{d}{dx_3}\f$ of the function
1436  * \f$f(\boldsymbol{x})\f$ at the same quadrature
1437  * points. The local distribution of the quadrature points
1438  * allows an elemental evaluation of the derivative. This
1439  * is done by a call to the function
1440  * StdRegions#StdExpansion#PhysDeriv.
1441  *
1442  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1443  * containing the values of the function
1444  * \f$f(\boldsymbol{x})\f$ at the quadrature
1445  * points \f$\boldsymbol{x}_i\f$.
1446  * @param out_d0 The discrete evaluation of the
1447  * derivative\f$\frac{d}{dx_1}\f$ will
1448  * be stored in this array of size
1449  * \f$Q_{\mathrm{tot}}\f$.
1450  * @param out_d1 The discrete evaluation of the
1451  * derivative\f$\frac{d}{dx_2}\f$ will be
1452  * stored in this array of size
1453  * \f$Q_{\mathrm{tot}}\f$. Note that if no
1454  * memory is allocated for \a out_d1, the
1455  * derivative \f$\frac{d}{dx_2}\f$ will not be
1456  * calculated.
1457  * @param out_d2 The discrete evaluation of the
1458  * derivative\f$\frac{d}{dx_3}\f$ will be
1459  * stored in this array of size
1460  * \f$Q_{\mathrm{tot}}\f$. Note that if no
1461  * memory is allocated for \a out_d2, the
1462  * derivative \f$\frac{d}{dx_3}\f$ will not be
1463  * calculated.
1464  */
1465  void ExpList::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
1466  Array<OneD, NekDouble> &out_d0,
1467  Array<OneD, NekDouble> &out_d1,
1468  Array<OneD, NekDouble> &out_d2)
1469  {
1470  Array<OneD, NekDouble> e_out_d0;
1471  Array<OneD, NekDouble> e_out_d1;
1472  Array<OneD, NekDouble> e_out_d2;
1473  int offset;
1474 
1475  // initialise if required
1476  if(m_collectionsDoInit[Collections::ePhysDeriv])
1477  {
1478  for (int i = 0; i < m_collections.size(); ++i)
1479  {
1480  m_collections[i].Initialise(Collections::ePhysDeriv);
1481  }
1482  m_collectionsDoInit[Collections::ePhysDeriv] = false;
1483  }
1484 
1485  LibUtilities::Timer timer;
1486  timer.Start();
1487  for (int i = 0; i < m_collections.size(); ++i)
1488  {
1489  offset = m_coll_phys_offset[i];
1490  e_out_d0 = out_d0 + offset;
1491  e_out_d1 = out_d1 + offset;
1492  e_out_d2 = out_d2 + offset;
1493 
1494 
1495  m_collections[i].ApplyOperator(Collections::ePhysDeriv,
1496  inarray + offset,
1497  e_out_d0,e_out_d1, e_out_d2);
1498 
1499 
1500  }
1501  timer.Stop();
1502  // Elapsed time
1503  timer.AccumulateRegion("Collections:PhysDeriv");
1504  }
1505 
1506  void ExpList::v_PhysDeriv(const int dir,
1507  const Array<OneD, const NekDouble> &inarray,
1508  Array<OneD, NekDouble> &out_d)
1509  {
1510  Direction edir = DirCartesianMap[dir];
1511  v_PhysDeriv(edir, inarray,out_d);
1512  }
1513 
1514  void ExpList::v_PhysDeriv(Direction edir, const Array<OneD, const NekDouble> &inarray,
1515  Array<OneD, NekDouble> &out_d)
1516  {
1517  int i;
1518  if(edir==MultiRegions::eS)
1519  {
1520  Array<OneD, NekDouble> e_out_ds;
1521  for(i=0; i<(*m_exp).size(); ++i)
1522  {
1523  e_out_ds = out_d + m_phys_offset[i];
1524  (*m_exp)[i]->PhysDeriv_s(inarray+m_phys_offset[i],e_out_ds);
1525  }
1526  }
1527  else if(edir==MultiRegions::eN)
1528  {
1529  Array<OneD, NekDouble > e_out_dn;
1530  for(i=0; i<(*m_exp).size(); i++)
1531  {
1532  e_out_dn = out_d +m_phys_offset[i];
1533  (*m_exp)[i]->PhysDeriv_n(inarray+m_phys_offset[i],e_out_dn);
1534  }
1535  }
1536  else
1537  {
1538 
1539  // initialise if required
1540  if(m_collectionsDoInit[Collections::ePhysDeriv])
1541  {
1542  for (int i = 0; i < m_collections.size(); ++i)
1543  {
1544  m_collections[i].Initialise(Collections::ePhysDeriv);
1545  }
1546  m_collectionsDoInit[Collections::ePhysDeriv] = false;
1547  }
1548 
1549  // convert enum into int
1550  int intdir= (int)edir;
1551  Array<OneD, NekDouble> e_out_d;
1552  int offset;
1553  for (int i = 0; i < m_collections.size(); ++i)
1554  {
1555  offset = m_coll_phys_offset[i];
1556  e_out_d = out_d + offset;
1557 
1558  m_collections[i].ApplyOperator(Collections::ePhysDeriv,
1559  intdir,
1560  inarray + offset,
1561  e_out_d);
1562  }
1563  }
1564  }
1565 
1566  void ExpList::v_CurlCurl(
1567  Array<OneD, Array<OneD, NekDouble> > &Vel,
1568  Array<OneD, Array<OneD, NekDouble> > &Q)
1569  {
1570  int nq = GetTotPoints();
1571  Array<OneD,NekDouble> Vx(nq);
1572  Array<OneD,NekDouble> Uy(nq);
1573  Array<OneD,NekDouble> Dummy(nq);
1574 
1575  bool halfMode = false;
1576  if ( GetExpType() == e3DH1D)
1577  {
1578  m_session->MatchSolverInfo("ModeType", "HalfMode",
1579  halfMode, false);
1580  }
1581 
1582  switch(m_expType)
1583  {
1584  case e2D:
1585  {
1586  PhysDeriv(xDir, Vel[yDir], Vx);
1587  PhysDeriv(yDir, Vel[xDir], Uy);
1588 
1589 
1590  Vmath::Vsub(nq, Vx, 1, Uy, 1, Dummy, 1);
1591 
1592  PhysDeriv(Dummy,Q[1],Q[0]);
1593 
1594  Vmath::Smul(nq, -1.0, Q[1], 1, Q[1], 1);
1595  }
1596  break;
1597 
1598  case e3D:
1599  case e3DH1D:
1600  case e3DH2D:
1601  {
1602  Array<OneD,NekDouble> Vz(nq);
1603  Array<OneD,NekDouble> Uz(nq);
1604  Array<OneD,NekDouble> Wx(nq);
1605  Array<OneD,NekDouble> Wy(nq);
1606 
1607  PhysDeriv(Vel[xDir], Dummy, Uy, Uz);
1608  PhysDeriv(Vel[yDir], Vx, Dummy, Vz);
1609  PhysDeriv(Vel[zDir], Wx, Wy, Dummy);
1610 
1611  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
1612  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
1613  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
1614 
1615  PhysDeriv(Q[0], Dummy, Uy, Uz);
1616  PhysDeriv(Q[1], Vx, Dummy, Vz);
1617  PhysDeriv(Q[2], Wx, Wy, Dummy);
1618 
1619  // For halfmode, need to change the sign of z derivatives
1620  if (halfMode)
1621  {
1622  Vmath::Neg(nq, Uz, 1);
1623  Vmath::Neg(nq, Vz, 1);
1624  }
1625 
1626  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
1627  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
1628  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
1629  }
1630  break;
1631  default:
1632  ASSERTL0(0,"Dimension not supported");
1633  break;
1634  }
1635  }
1636 
1637 
1638  void ExpList::v_PhysDirectionalDeriv(
1639  const Array<OneD, const NekDouble> &direction,
1640  const Array<OneD, const NekDouble> &inarray,
1641  Array<OneD, NekDouble> &outarray)
1642  {
1643  int npts_e;
1644  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1645  int nq = direction.size() / coordim;
1646 
1647  Array<OneD, NekDouble> e_outarray;
1648  Array<OneD, NekDouble> e_MFdiv;
1649  Array<OneD, NekDouble> locdir;
1650 
1651  for(int i = 0; i < (*m_exp).size(); ++i)
1652  {
1653  npts_e = (*m_exp)[i]->GetTotPoints();
1654  locdir = Array<OneD, NekDouble>(npts_e*coordim);
1655 
1656  for (int k = 0; k<coordim; ++k)
1657  {
1658  Vmath::Vcopy(npts_e, &direction[k*nq+m_phys_offset[i]], 1,
1659  &locdir[k*npts_e], 1);
1660  }
1661 
1662  (*m_exp)[i]->PhysDirectionalDeriv(
1663  inarray + m_phys_offset[i],
1664  locdir,
1665  e_outarray = outarray+m_phys_offset[i] );
1666  }
1667  }
1668 
1669 
1670  void ExpList::ExponentialFilter(
1671  Array<OneD, NekDouble> &array,
1672  const NekDouble alpha,
1673  const NekDouble exponent,
1674  const NekDouble cutoff)
1675  {
1676  Array<OneD,NekDouble> e_array;
1677 
1678  for(int i = 0; i < (*m_exp).size(); ++i)
1679  {
1680  (*m_exp)[i]->ExponentialFilter(
1681  e_array = array+m_phys_offset[i],
1682  alpha,
1683  exponent,
1684  cutoff);
1685  }
1686  }
1687 
1688  /**
1689  * The coefficients of the function to be acted upon
1690  * should be contained in the \param inarray. The
1691  * resulting coefficients are stored in \param outarray
1692  *
1693  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
1694  * containing the inner product.
1695  */
1696  void ExpList::MultiplyByElmtInvMass(
1697  const Array<OneD, const NekDouble> &inarray,
1698  Array<OneD, NekDouble> &outarray)
1699  {
1700  GlobalMatrixKey mkey(StdRegions::eInvMass);
1701  const DNekScalBlkMatSharedPtr& InvMass = GetBlockMatrix(mkey);
1702 
1703  // Inverse mass matrix
1704  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
1705  if(inarray.get() == outarray.get())
1706  {
1707  NekVector<NekDouble> in(m_ncoeffs,inarray); // copy data
1708  out = (*InvMass)*in;
1709  }
1710  else
1711  {
1712  NekVector<NekDouble> in(m_ncoeffs,inarray,eWrapper);
1713  out = (*InvMass)*in;
1714  }
1715  }
1716 
1717  /**
1718  * Given a function \f$u(\boldsymbol{x})\f$ defined at the
1719  * quadrature points, this function determines the
1720  * transformed elemental coefficients \f$\hat{u}_n^e\f$
1721  * employing a discrete elemental Galerkin projection from
1722  * physical space to coefficient space. For each element,
1723  * the operation is evaluated locally by the function
1724  * StdRegions#StdExpansion#IproductWRTBase followed by a
1725  * call to #MultiRegions#MultiplyByElmtInvMass.
1726  *
1727  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
1728  * containing the values of the function
1729  * \f$f(\boldsymbol{x})\f$ at the quadrature
1730  * points \f$\boldsymbol{x}_i\f$.
1731  * @param outarray The resulting coefficients
1732  * \f$\hat{u}_n^e\f$ will be stored in this
1733  * array of size \f$N_{\mathrm{eof}}\f$.
1734  */
1735  void ExpList::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
1736  Array<OneD, NekDouble> &outarray)
1737  {
1738  Array<OneD,NekDouble> f(m_ncoeffs);
1739 
1740  IProductWRTBase_IterPerExp(inarray,f);
1741  MultiplyByElmtInvMass(f,outarray);
1742 
1743  }
1744 
1745  void ExpList::v_FwdTrans_BndConstrained(
1746  const Array<OneD, const NekDouble>& inarray,
1747  Array<OneD, NekDouble> &outarray)
1748  {
1749  int i;
1750 
1751  Array<OneD,NekDouble> e_outarray;
1752 
1753  for(i= 0; i < (*m_exp).size(); ++i)
1754  {
1755  (*m_exp)[i]->FwdTrans_BndConstrained(inarray+m_phys_offset[i],
1756  e_outarray = outarray+m_coeff_offset[i]);
1757  }
1758  }
1759 
1760  /**
1761  * This function smooth a field after some calculaitons which have
1762  * been done elementally.
1763  *
1764  * @param field An array containing the field in physical space
1765  *
1766  */
1767  void ExpList::v_SmoothField(Array<OneD, NekDouble> &field)
1768  {
1769  boost::ignore_unused(field);
1770  // Do nothing unless the method is implemented in the appropriate
1771  // class, i.e. ContField1D,ContField2D, etc.
1772 
1773  // So far it has been implemented just for ContField2D and
1774  // ContField3DHomogeneous1D
1775 
1776  // Block in case users try the smoothing with a modal expansion.
1777  // Maybe a different techique for the smoothing require
1778  // implementation for modal basis.
1779 
1780  ASSERTL0((*m_exp)[0]->GetBasisType(0)
1781  == LibUtilities::eGLL_Lagrange ||
1782  (*m_exp)[0]->GetBasisType(0)
1783  == LibUtilities::eGauss_Lagrange,
1784  "Smoothing is currently not allowed unless you are using "
1785  "a nodal base for efficiency reasons. The implemented "
1786  "smoothing technique requires the mass matrix inversion "
1787  "which is trivial just for GLL_LAGRANGE_SEM and "
1788  "GAUSS_LAGRANGE_SEMexpansions.");
1789  }
1790 
1791 
1792  /**
1793  * This function assembles the block diagonal matrix
1794  * \f$\underline{\boldsymbol{M}}^e\f$, which is the
1795  * concatenation of the local matrices
1796  * \f$\boldsymbol{M}^e\f$ of the type \a mtype, that is
1797  *
1798  * \f[
1799  * \underline{\boldsymbol{M}}^e = \left[
1800  * \begin{array}{cccc}
1801  * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
1802  * 0 & \boldsymbol{M}^2 & 0 & 0 \\
1803  * 0 & 0 & \ddots & 0 \\
1804  * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array}\right].\f]
1805  *
1806  * @param mtype the type of matrix to be assembled
1807  * @param scalar an optional parameter
1808  * @param constant an optional parameter
1809  */
1810  const DNekScalBlkMatSharedPtr ExpList::GenBlockMatrix(
1811  const GlobalMatrixKey &gkey)
1812  {
1813  int i,cnt1;
1814  int n_exp = 0;
1815  DNekScalMatSharedPtr loc_mat;
1816  DNekScalBlkMatSharedPtr BlkMatrix;
1817  map<int,int> elmt_id;
1818  LibUtilities::ShapeType ShapeType = gkey.GetShapeType();
1819 
1820  if(ShapeType != LibUtilities::eNoShapeType)
1821  {
1822  for(i = 0 ; i < (*m_exp).size(); ++i)
1823  {
1824  if((*m_exp)[i]->DetShapeType()
1825  == ShapeType)
1826  {
1827  elmt_id[n_exp++] = i;
1828  }
1829  }
1830  }
1831  else
1832  {
1833  n_exp = (*m_exp).size();
1834  for(i = 0; i < n_exp; ++i)
1835  {
1836  elmt_id[i] = i;
1837  }
1838  }
1839 
1840  Array<OneD,unsigned int> nrows(n_exp);
1841  Array<OneD,unsigned int> ncols(n_exp);
1842 
1843  switch(gkey.GetMatrixType())
1844  {
1845  case StdRegions::eBwdTrans:
1846  {
1847  // set up an array of integers for block matrix construction
1848  for(i = 0; i < n_exp; ++i)
1849  {
1850  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1851  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1852  }
1853  }
1854  break;
1855  case StdRegions::eIProductWRTBase:
1856  {
1857  // set up an array of integers for block matrix construction
1858  for(i = 0; i < n_exp; ++i)
1859  {
1860  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1861  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1862  }
1863  }
1864  break;
1865  case StdRegions::eMass:
1866  case StdRegions::eInvMass:
1867  case StdRegions::eHelmholtz:
1868  case StdRegions::eLaplacian:
1869  case StdRegions::eInvHybridDGHelmholtz:
1870  {
1871  // set up an array of integers for block matrix construction
1872  for(i = 0; i < n_exp; ++i)
1873  {
1874  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1875  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1876  }
1877  }
1878  break;
1879 
1880  case StdRegions::eHybridDGLamToU:
1881  {
1882  // set up an array of integers for block matrix construction
1883  for(i = 0; i < n_exp; ++i)
1884  {
1885  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1886  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1887  }
1888  }
1889  break;
1890 
1891  case StdRegions::eHybridDGHelmBndLam:
1892  {
1893  // set up an array of integers for block matrix construction
1894  for(i = 0; i < n_exp; ++i)
1895  {
1896  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1897  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1898  }
1899  }
1900  break;
1901 
1902  default:
1903  {
1904  NEKERROR(ErrorUtil::efatal,
1905  "Global Matrix creation not defined for this "
1906  "type of matrix");
1907  }
1908  }
1909 
1910  MatrixStorage blkmatStorage = eDIAGONAL;
1911  BlkMatrix = MemoryManager<DNekScalBlkMat>
1912  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
1913 
1914  int nvarcoeffs = gkey.GetNVarCoeffs();
1915  int eid;
1916  Array<OneD, NekDouble> varcoeffs_wk;
1917 
1918  for(i = cnt1 = 0; i < n_exp; ++i)
1919  {
1920  // need to be initialised with zero size for non
1921  // variable coefficient case
1922  StdRegions::VarCoeffMap varcoeffs;
1923 
1924  eid = elmt_id[i];
1925  if(nvarcoeffs>0)
1926  {
1927  for (auto &x : gkey.GetVarCoeffs())
1928  {
1929  varcoeffs[x.first] = x.second + m_phys_offset[eid];
1930  }
1931  }
1932 
1933  LocalRegions::MatrixKey matkey(gkey.GetMatrixType(),
1934  (*m_exp)[eid]->DetShapeType(),
1935  *(*m_exp)[eid],
1936  gkey.GetConstFactors(),
1937  varcoeffs );
1938 
1939  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[elmt_id.find(i)->second])->GetLocMatrix(matkey);
1940  BlkMatrix->SetBlock(i,i,loc_mat);
1941  }
1942 
1943  return BlkMatrix;
1944  }
1945 
1946  const DNekScalBlkMatSharedPtr& ExpList::GetBlockMatrix(
1947  const GlobalMatrixKey &gkey)
1948  {
1949  auto matrixIter = m_blockMat->find(gkey);
1950 
1951  if(matrixIter == m_blockMat->end())
1952  {
1953  return ((*m_blockMat)[gkey] = GenBlockMatrix(gkey));
1954  }
1955  else
1956  {
1957  return matrixIter->second;
1958  }
1959  }
1960 
1961  void ExpList::GeneralMatrixOp_IterPerExp(
1962  const GlobalMatrixKey &gkey,
1963  const Array<OneD,const NekDouble> &inarray,
1964  Array<OneD, NekDouble> &outarray)
1965  {
1966  int nvarcoeffs = gkey.GetNVarCoeffs();
1967 
1968  if((nvarcoeffs == 0)&&(gkey.GetMatrixType() == StdRegions::eHelmholtz))
1969  {
1970  // initialise if required
1971  if(m_collections.size()&&m_collectionsDoInit[Collections::eHelmholtz])
1972  {
1973  for (int i = 0; i < m_collections.size(); ++i)
1974  {
1975  m_collections[i].Initialise(Collections::eHelmholtz, gkey.GetConstFactors());
1976  }
1977  m_collectionsDoInit[Collections::eHelmholtz] = false;
1978  }
1979  else
1980  {
1981  for (int i = 0; i < m_collections.size(); ++i)
1982  {
1983  m_collections[i].CheckFactors(Collections::eHelmholtz, gkey.GetConstFactors(),
1984  m_coll_phys_offset[i]);
1985  }
1986  }
1987 
1988  Array<OneD, NekDouble> tmp;
1989  for (int i = 0; i < m_collections.size(); ++i)
1990  {
1991  m_collections[i].ApplyOperator
1992  (Collections::eHelmholtz,
1993  inarray + m_coll_coeff_offset[i],
1994  tmp = outarray + m_coll_coeff_offset[i]);
1995  }
1996  }
1997  else
1998  {
1999  Array<OneD,NekDouble> tmp_outarray;
2000  for(int i= 0; i < (*m_exp).size(); ++i)
2001  {
2002  // need to be initialised with zero size for non
2003  // variable coefficient case
2004  StdRegions::VarCoeffMap varcoeffs;
2005 
2006  if(nvarcoeffs>0)
2007  {
2008  for (auto &x : gkey.GetVarCoeffs())
2009  {
2010  varcoeffs[x.first] = x.second + m_phys_offset[i];
2011  }
2012  }
2013 
2014  StdRegions::StdMatrixKey mkey(gkey.GetMatrixType(),
2015  (*m_exp)[i]->DetShapeType(),
2016  *((*m_exp)[i]),
2017  gkey.GetConstFactors(),varcoeffs);
2018 
2019  (*m_exp)[i]->GeneralMatrixOp(inarray + m_coeff_offset[i],
2020  tmp_outarray = outarray+
2021  m_coeff_offset[i],
2022  mkey);
2023  }
2024  }
2025  }
2026 
2027  /**
2028  * Retrieves local matrices from each expansion in the expansion list
2029  * and combines them together to generate a global matrix system.
2030  * @param mkey Matrix key for the matrix to be generated.
2031  * @param locToGloMap Local to global mapping.
2032  * @returns Shared pointer to the generated global matrix.
2033  */
2034  GlobalMatrixSharedPtr ExpList::GenGlobalMatrix(
2035  const GlobalMatrixKey &mkey,
2036  const AssemblyMapCGSharedPtr &locToGloMap)
2037  {
2038  int i,j,n,gid1,gid2,cntdim1,cntdim2;
2039  NekDouble sign1,sign2;
2040  DNekScalMatSharedPtr loc_mat;
2041 
2042  unsigned int glob_rows = 0;
2043  unsigned int glob_cols = 0;
2044  unsigned int loc_rows = 0;
2045  unsigned int loc_cols = 0;
2046 
2047  bool assembleFirstDim = false;
2048  bool assembleSecondDim = false;
2049 
2050  switch(mkey.GetMatrixType())
2051  {
2052  case StdRegions::eBwdTrans:
2053  {
2054  glob_rows = m_npoints;
2055  glob_cols = locToGloMap->GetNumGlobalCoeffs();
2056 
2057  assembleFirstDim = false;
2058  assembleSecondDim = true;
2059  }
2060  break;
2061  case StdRegions::eIProductWRTBase:
2062  {
2063  glob_rows = locToGloMap->GetNumGlobalCoeffs();
2064  glob_cols = m_npoints;
2065 
2066  assembleFirstDim = true;
2067  assembleSecondDim = false;
2068  }
2069  break;
2070  case StdRegions::eMass:
2071  case StdRegions::eHelmholtz:
2072  case StdRegions::eLaplacian:
2073  case StdRegions::eHybridDGHelmBndLam:
2074  {
2075  glob_rows = locToGloMap->GetNumGlobalCoeffs();
2076  glob_cols = locToGloMap->GetNumGlobalCoeffs();
2077 
2078  assembleFirstDim = true;
2079  assembleSecondDim = true;
2080  }
2081  break;
2082  default:
2083  {
2084  NEKERROR(ErrorUtil::efatal,
2085  "Global Matrix creation not defined for this "
2086  "type of matrix");
2087  }
2088  }
2089 
2090  COOMatType spcoomat;
2091  CoordType coord;
2092 
2093  int nvarcoeffs = mkey.GetNVarCoeffs();
2094  int eid;
2095 
2096  // fill global matrix
2097  for(n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
2098  {
2099  // need to be initialised with zero size for non variable coefficient case
2100  StdRegions::VarCoeffMap varcoeffs;
2101 
2102  eid = n;
2103  if(nvarcoeffs>0)
2104  {
2105  for (auto &x : mkey.GetVarCoeffs())
2106  {
2107  varcoeffs[x.first] = x.second + m_phys_offset[eid];
2108  }
2109  }
2110 
2111  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
2112  (*m_exp)[eid]->DetShapeType(),
2113  *((*m_exp)[eid]),
2114  mkey.GetConstFactors(),varcoeffs);
2115 
2116  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])->GetLocMatrix(matkey);
2117 
2118  loc_rows = loc_mat->GetRows();
2119  loc_cols = loc_mat->GetColumns();
2120 
2121  for(i = 0; i < loc_rows; ++i)
2122  {
2123  if(assembleFirstDim)
2124  {
2125  gid1 = locToGloMap->GetLocalToGlobalMap (cntdim1 + i);
2126  sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
2127  }
2128  else
2129  {
2130  gid1 = cntdim1 + i;
2131  sign1 = 1.0;
2132  }
2133 
2134  for(j = 0; j < loc_cols; ++j)
2135  {
2136  if(assembleSecondDim)
2137  {
2138  gid2 = locToGloMap
2139  ->GetLocalToGlobalMap(cntdim2 + j);
2140  sign2 = locToGloMap
2141  ->GetLocalToGlobalSign(cntdim2 + j);
2142  }
2143  else
2144  {
2145  gid2 = cntdim2 + j;
2146  sign2 = 1.0;
2147  }
2148 
2149  // sparse matrix fill
2150  coord = make_pair(gid1,gid2);
2151  if( spcoomat.count(coord) == 0 )
2152  {
2153  spcoomat[coord] = sign1*sign2*(*loc_mat)(i,j);
2154  }
2155  else
2156  {
2157  spcoomat[coord] += sign1*sign2*(*loc_mat)(i,j);
2158  }
2159  }
2160  }
2161  cntdim1 += loc_rows;
2162  cntdim2 += loc_cols;
2163  }
2164 
2165  return MemoryManager<GlobalMatrix>
2166  ::AllocateSharedPtr(m_session,glob_rows,glob_cols,spcoomat);
2167  }
2168 
2169 
2170  DNekMatSharedPtr ExpList::GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
2171  {
2172  int i,j,n,gid1,gid2,loc_lda,eid;
2173  NekDouble sign1,sign2,value;
2174  DNekScalMatSharedPtr loc_mat;
2175 
2176  int totDofs = locToGloMap->GetNumGlobalCoeffs();
2177  int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
2178 
2179  unsigned int rows = totDofs - NumDirBCs;
2180  unsigned int cols = totDofs - NumDirBCs;
2181  NekDouble zero = 0.0;
2182 
2183  DNekMatSharedPtr Gmat;
2184  int bwidth = locToGloMap->GetFullSystemBandWidth();
2185 
2186  int nvarcoeffs = mkey.GetNVarCoeffs();
2187  MatrixStorage matStorage;
2188 
2189  map<int, RobinBCInfoSharedPtr> RobinBCInfo = GetRobinBCInfo();
2190 
2191  switch(mkey.GetMatrixType())
2192  {
2193  // case for all symmetric matices
2194  case StdRegions::eHelmholtz:
2195  case StdRegions::eLaplacian:
2196  if( (2*(bwidth+1)) < rows)
2197  {
2198  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC_BANDED;
2199  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage,bwidth,bwidth);
2200  }
2201  else
2202  {
2203  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
2204  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
2205  }
2206 
2207  break;
2208  default: // Assume general matrix - currently only set up for full invert
2209  {
2210  matStorage = eFULL;
2211  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
2212  }
2213  }
2214 
2215  // fill global symmetric matrix
2216  for(n = 0; n < (*m_exp).size(); ++n)
2217  {
2218  // need to be initialised with zero size for non variable coefficient case
2219  StdRegions::VarCoeffMap varcoeffs;
2220 
2221  eid = n;
2222  if(nvarcoeffs>0)
2223  {
2224  for (auto &x : mkey.GetVarCoeffs())
2225  {
2226  varcoeffs[x.first] = x.second + m_phys_offset[eid];
2227  }
2228  }
2229 
2230  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
2231  (*m_exp)[eid]->DetShapeType(),
2232  *((*m_exp)[eid]),
2233  mkey.GetConstFactors(),varcoeffs);
2234 
2235  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])->GetLocMatrix(matkey);
2236 
2237 
2238  if(RobinBCInfo.count(n) != 0) // add robin mass matrix
2239  {
2240  RobinBCInfoSharedPtr rBC;
2241 
2242  // declare local matrix from scaled matrix.
2243  int rows = loc_mat->GetRows();
2244  int cols = loc_mat->GetColumns();
2245  const NekDouble *dat = loc_mat->GetRawPtr();
2246  DNekMatSharedPtr new_mat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,dat);
2247  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
2248 
2249  // add local matrix contribution
2250  for(rBC = RobinBCInfo.find(n)->second;rBC; rBC = rBC->next)
2251  {
2252  (*m_exp)[n]->AddRobinMassMatrix(rBC->m_robinID,rBC->m_robinPrimitiveCoeffs,new_mat);
2253  }
2254 
2255  NekDouble one = 1.0;
2256  // redeclare loc_mat to point to new_mat plus the scalar.
2257  loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,new_mat);
2258  }
2259 
2260  loc_lda = loc_mat->GetColumns();
2261 
2262  for(i = 0; i < loc_lda; ++i)
2263  {
2264  gid1 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + i) - NumDirBCs;
2265  sign1 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + i);
2266  if(gid1 >= 0)
2267  {
2268  for(j = 0; j < loc_lda; ++j)
2269  {
2270  gid2 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + j) - NumDirBCs;
2271  sign2 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + j);
2272  if(gid2 >= 0)
2273  {
2274  // When global matrix is symmetric,
2275  // only add the value for the upper
2276  // triangular part in order to avoid
2277  // entries to be entered twice
2278  if((matStorage == eFULL)||(gid2 >= gid1))
2279  {
2280  value = Gmat->GetValue(gid1,gid2) + sign1*sign2*(*loc_mat)(i,j);
2281  Gmat->SetValue(gid1,gid2,value);
2282  }
2283  }
2284  }
2285  }
2286  }
2287  }
2288 
2289  return Gmat;
2290  }
2291 
2292 
2293  /**
2294  * Consider a linear system
2295  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{f}\f$ to be solved. Dependent
2296  * on the solution method, this function constructs
2297  * - <b>The full linear system</b><BR>
2298  * A call to the function #GenGlobalLinSysFullDirect
2299  * - <b>The statically condensed linear system</b><BR>
2300  * A call to the function #GenGlobalLinSysStaticCond
2301  *
2302  * @param mkey A key which uniquely defines the global
2303  * matrix to be constructed.
2304  * @param locToGloMap Contains the mapping array and required
2305  * information for the transformation from
2306  * local to global degrees of freedom.
2307  * @return (A shared pointer to) the global linear system in
2308  * required format.
2309  */
2310  GlobalLinSysSharedPtr ExpList::GenGlobalLinSys(
2311  const GlobalLinSysKey &mkey,
2312  const AssemblyMapCGSharedPtr &locToGloMap)
2313  {
2314  GlobalLinSysSharedPtr returnlinsys;
2315  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
2316 
2317  MultiRegions::GlobalSysSolnType vType = mkey.GetGlobalSysSolnType();
2318 
2319  if (vType >= eSIZE_GlobalSysSolnType)
2320  {
2321  NEKERROR(ErrorUtil::efatal,"Matrix solution type not defined");
2322  }
2323  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
2324 
2325  return GetGlobalLinSysFactory().CreateInstance( vSolnType, mkey,
2326  vExpList, locToGloMap);
2327  }
2328 
2329  GlobalLinSysSharedPtr ExpList::GenGlobalBndLinSys(
2330  const GlobalLinSysKey &mkey,
2331  const AssemblyMapSharedPtr &locToGloMap)
2332  {
2333  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
2334  const map<int,RobinBCInfoSharedPtr> vRobinBCInfo = GetRobinBCInfo();
2335 
2336  MultiRegions::GlobalSysSolnType vType = mkey.GetGlobalSysSolnType();
2337 
2338  if (vType >= eSIZE_GlobalSysSolnType)
2339  {
2340  NEKERROR(ErrorUtil::efatal,"Matrix solution type not defined");
2341  }
2342  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
2343 
2344  return GetGlobalLinSysFactory().CreateInstance(vSolnType,mkey,
2345  vExpList,locToGloMap);
2346  }
2347 
2348 
2349  /**
2350  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of
2351  * an expansion, this function evaluates the spectral/hp
2352  * expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
2353  * quadrature points \f$\boldsymbol{x}_i\f$. The operation
2354  * is evaluated locally by the elemental function
2355  * StdRegions#StdExpansion#BwdTrans.
2356  *
2357  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
2358  * containing the local coefficients
2359  * \f$\hat{u}_n^e\f$.
2360  * @param outarray The resulting physical values at the
2361  * quadrature points
2362  * \f$u^{\delta}(\boldsymbol{x}_i)\f$
2363  * will be stored in this array of size
2364  * \f$Q_{\mathrm{tot}}\f$.
2365  */
2366  void ExpList::v_BwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
2367  Array<OneD, NekDouble> &outarray)
2368  {
2369  LibUtilities::Timer timer;
2370 
2371  if(m_expType == e0D)
2372  {
2373  Vmath::Vcopy(m_ncoeffs,inarray,1,outarray,1);
2374  }
2375  else
2376  {
2377  // initialise if required
2378  if(m_collections.size()&&m_collectionsDoInit[Collections::eBwdTrans])
2379  {
2380  for (int i = 0; i < m_collections.size(); ++i)
2381  {
2382  m_collections[i].Initialise(Collections::eBwdTrans);
2383  }
2384  m_collectionsDoInit[Collections::eBwdTrans] = false;
2385  }
2386 
2387  LIKWID_MARKER_START("v_BwdTrans_IterPerExp");
2388  // timer.Start();
2389 
2390  Array<OneD, NekDouble> tmp;
2391  for (int i = 0; i < m_collections.size(); ++i)
2392  {
2393  m_collections[i].ApplyOperator(Collections::eBwdTrans,
2394  inarray + m_coll_coeff_offset[i],
2395  tmp = outarray + m_coll_phys_offset[i]);
2396  }
2397  }
2398 
2399  timer.Stop();
2400  LIKWID_MARKER_STOP("v_BwdTrans_IterPerExp");
2401 
2402  // Elapsed time
2403  timer.AccumulateRegion("Collections:BwdTrans");
2404  }
2405 
2406  LocalRegions::ExpansionSharedPtr& ExpList::GetExp(
2407  const Array<OneD, const NekDouble> &gloCoord)
2408  {
2409  return GetExp(GetExpIndex(gloCoord));
2410  }
2411 
2412 
2413  /**
2414  * @todo need a smarter search here that first just looks at bounding
2415  * vertices - suggest first seeing if point is within 10% of
2416  * region defined by vertices. The do point search.
2417  */
2418  int ExpList::GetExpIndex(
2419  const Array<OneD, const NekDouble> &gloCoord,
2420  NekDouble tol,
2421  bool returnNearestElmt,
2422  int cachedId,
2423  NekDouble maxDistance)
2424  {
2425  Array<OneD, NekDouble> Lcoords(gloCoord.size());
2426 
2427  return GetExpIndex(gloCoord,Lcoords,tol,returnNearestElmt,cachedId,maxDistance);
2428  }
2429 
2430 
2431  int ExpList::GetExpIndex(
2432  const Array<OneD, const NekDouble> &gloCoords,
2433  Array<OneD, NekDouble> &locCoords,
2434  NekDouble tol,
2435  bool returnNearestElmt,
2436  int cachedId,
2437  NekDouble maxDistance)
2438  {
2439  if (GetNumElmts() == 0)
2440  {
2441  return -1;
2442  }
2443 
2444  if (m_elmtToExpId.size() == 0)
2445  {
2446  // Loop in reverse order so that in case where using a
2447  // Homogeneous expansion it sets geometry ids to first part of
2448  // m_exp list. Otherwise will set to second (complex) expansion
2449  for(int i = (*m_exp).size()-1; i >= 0; --i)
2450  {
2451  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2452  }
2453  }
2454 
2455  NekDouble nearpt = 1e6;
2456  NekDouble nearpt_min = 1e6;
2457  int min_id = -1;
2458  Array<OneD, NekDouble> savLocCoords(locCoords.size());
2459 
2460  if(cachedId >= 0 && cachedId < (*m_exp).size())
2461  {
2462  nearpt = 1e12;
2463  if((*m_exp)[cachedId]->GetGeom()->MinMaxCheck(gloCoords) &&
2464  (*m_exp)[cachedId]->GetGeom()->ContainsPoint(gloCoords,
2465  locCoords, tol, nearpt))
2466  {
2467  return cachedId;
2468  }
2469  else if(returnNearestElmt && (nearpt < nearpt_min))
2470  {
2471  // If it does not lie within, keep track of which element
2472  // is nearest.
2473  min_id = cachedId;
2474  nearpt_min = nearpt;
2475  Vmath::Vcopy(locCoords.size(),locCoords, 1, savLocCoords, 1);
2476  }
2477  }
2478 
2479  NekDouble x = (gloCoords.size() > 0 ? gloCoords[0] : 0.0);
2480  NekDouble y = (gloCoords.size() > 1 ? gloCoords[1] : 0.0);
2481  NekDouble z = (gloCoords.size() > 2 ? gloCoords[2] : 0.0);
2482  SpatialDomains::PointGeomSharedPtr p
2483  = MemoryManager<SpatialDomains::PointGeom>::AllocateSharedPtr(
2484  GetExp(0)->GetCoordim(), -1, x, y, z);
2485 
2486  // Get the list of elements whose bounding box contains the desired
2487  // point.
2488  std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
2489 
2490  // Check each element in turn to see if point lies within it.
2491  for (int i = 0; i < elmts.size(); ++i)
2492  {
2493  int id = m_elmtToExpId[elmts[i]];
2494  if(id == cachedId)
2495  {
2496  continue;
2497  }
2498  if ((*m_exp)[id]->GetGeom()->ContainsPoint(gloCoords,
2499  locCoords, tol, nearpt))
2500  {
2501  return id;
2502  }
2503  else if(returnNearestElmt && (nearpt < nearpt_min))
2504  {
2505  // If it does not lie within, keep track of which element
2506  // is nearest.
2507  min_id = id;
2508  nearpt_min = nearpt;
2509  Vmath::Vcopy(locCoords.size(),locCoords, 1, savLocCoords, 1);
2510  }
2511  }
2512 
2513  // If the calling function is with just the nearest element, return
2514  // that. Otherwise return -1 to indicate no matching elemenet found.
2515  if(returnNearestElmt && nearpt_min <= maxDistance)
2516  {
2517 
2518  std::string msg = "Failed to find point within element to "
2519  "tolerance of "
2520  + boost::lexical_cast<std::string>(tol)
2521  + " using local point ("
2522  + boost::lexical_cast<std::string>(locCoords[0]) +","
2523  + boost::lexical_cast<std::string>(locCoords[1]) +","
2524  + boost::lexical_cast<std::string>(locCoords[1])
2525  + ") in element: "
2526  + boost::lexical_cast<std::string>(min_id);
2527  WARNINGL1(false,msg.c_str());
2528 
2529  Vmath::Vcopy(locCoords.size(),savLocCoords, 1, locCoords, 1);
2530  return min_id;
2531  }
2532  else
2533  {
2534  return -1;
2535  }
2536  }
2537 
2538  /**
2539  * Given some coordinates, output the expansion field value at that
2540  * point
2541  */
2542  NekDouble ExpList::PhysEvaluate(
2543  const Array<OneD, const NekDouble> &coords,
2544  const Array<OneD, const NekDouble> &phys)
2545  {
2546  int dim = GetCoordim(0);
2547  ASSERTL0(dim == coords.size(),
2548  "Invalid coordinate dimension.");
2549 
2550  // Grab the element index corresponding to coords.
2551  Array<OneD, NekDouble> xi(dim);
2552  int elmtIdx = GetExpIndex(coords, xi);
2553  ASSERTL0(elmtIdx > 0, "Unable to find element containing point.");
2554 
2555  // Grab that element's physical storage.
2556  Array<OneD, NekDouble> elmtPhys = phys + m_phys_offset[elmtIdx];
2557 
2558  // Evaluate the element at the appropriate point.
2559  return (*m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
2560  }
2561 
2562  /**
2563  * Configures geometric info, such as tangent direction, on each
2564  * expansion.
2565  * @param graph2D Mesh
2566  */
2567  void ExpList::ApplyGeomInfo()
2568  {
2569 
2570  }
2571 
2572  /**
2573  * @brief Reset geometry information, metrics, matrix managers and
2574  * geometry information.
2575  *
2576  * This routine clears all matrix managers and resets all geometry
2577  * information, which allows the geometry information to be dynamically
2578  * updated as the solver is run.
2579  */
2580  void ExpList::v_Reset()
2581  {
2582  // Reset matrix managers.
2583  LibUtilities::NekManager<LocalRegions::MatrixKey,
2584  DNekScalMat, LocalRegions::MatrixKey::opLess>::ClearManager();
2585  LibUtilities::NekManager<LocalRegions::MatrixKey,
2586  DNekScalBlkMat, LocalRegions::MatrixKey::opLess>::ClearManager();
2587 
2588  // Loop over all elements and reset geometry information.
2589  for (int i = 0; i < m_exp->size(); ++i)
2590  {
2591  (*m_exp)[i]->GetGeom()->Reset(m_graph->GetCurvedEdges(),
2592  m_graph->GetCurvedFaces());
2593  }
2594 
2595  // Loop over all elements and rebuild geometric factors.
2596  for (int i = 0; i < m_exp->size(); ++i)
2597  {
2598  (*m_exp)[i]->Reset();
2599  }
2600  }
2601 
2602  /**
2603  * Write Tecplot Files Header
2604  * @param outfile Output file name.
2605  * @param var variables names
2606  */
2607  void ExpList::v_WriteTecplotHeader(std::ostream &outfile,
2608  std::string var)
2609  {
2610  if (GetNumElmts() == 0)
2611  {
2612  return;
2613  }
2614 
2615  int coordim = GetExp(0)->GetCoordim();
2616  char vars[3] = { 'x', 'y', 'z' };
2617 
2618  if (m_expType == e3DH1D)
2619  {
2620  coordim += 1;
2621  }
2622  else if (m_expType == e3DH2D)
2623  {
2624  coordim += 2;
2625  }
2626 
2627  outfile << "Variables = x";
2628  for (int i = 1; i < coordim; ++i)
2629  {
2630  outfile << ", " << vars[i];
2631  }
2632 
2633  if (var.size() > 0)
2634  {
2635  outfile << ", " << var;
2636  }
2637 
2638  outfile << std::endl << std::endl;
2639  }
2640 
2641  /**
2642  * Write Tecplot Files Zone
2643  * @param outfile Output file name.
2644  * @param expansion Expansion that is considered
2645  */
2646  void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
2647  {
2648  int i, j;
2649  int coordim = GetCoordim(0);
2650  int nPoints = GetTotPoints();
2651  int nBases = (*m_exp)[0]->GetNumBases();
2652  int numBlocks = 0;
2653 
2654  Array<OneD, Array<OneD, NekDouble> > coords(3);
2655 
2656  if (expansion == -1)
2657  {
2658  nPoints = GetTotPoints();
2659 
2660  coords[0] = Array<OneD, NekDouble>(nPoints);
2661  coords[1] = Array<OneD, NekDouble>(nPoints);
2662  coords[2] = Array<OneD, NekDouble>(nPoints);
2663 
2664  GetCoords(coords[0], coords[1], coords[2]);
2665 
2666  for (i = 0; i < m_exp->size(); ++i)
2667  {
2668  int numInt = 1;
2669 
2670  for (j = 0; j < nBases; ++j)
2671  {
2672  numInt *= (*m_exp)[i]->GetNumPoints(j)-1;
2673  }
2674 
2675  numBlocks += numInt;
2676  }
2677  }
2678  else
2679  {
2680  nPoints = (*m_exp)[expansion]->GetTotPoints();
2681 
2682  coords[0] = Array<OneD, NekDouble>(nPoints);
2683  coords[1] = Array<OneD, NekDouble>(nPoints);
2684  coords[2] = Array<OneD, NekDouble>(nPoints);
2685 
2686  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
2687 
2688  numBlocks = 1;
2689  for (j = 0; j < nBases; ++j)
2690  {
2691  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j)-1;
2692  }
2693  }
2694 
2695  if (m_expType == e3DH1D)
2696  {
2697  nBases += 1;
2698  coordim += 1;
2699  int nPlanes = GetZIDs().size();
2700  NekDouble tmp = numBlocks * (nPlanes-1.0) / nPlanes;
2701  numBlocks = (int)tmp;
2702  }
2703  else if (m_expType == e3DH2D)
2704  {
2705  nBases += 2;
2706  coordim += 1;
2707  }
2708 
2709  outfile << "Zone, N=" << nPoints << ", E="
2710  << numBlocks << ", F=FEBlock" ;
2711 
2712  switch(nBases)
2713  {
2714  case 2:
2715  outfile << ", ET=QUADRILATERAL" << std::endl;
2716  break;
2717  case 3:
2718  outfile << ", ET=BRICK" << std::endl;
2719  break;
2720  default:
2721  NEKERROR(ErrorUtil::efatal,
2722  "Not set up for this type of output");
2723  break;
2724  }
2725 
2726  // Write out coordinates
2727  for (j = 0; j < coordim; ++j)
2728  {
2729  for (i = 0; i < nPoints; ++i)
2730  {
2731  outfile << coords[j][i] << " ";
2732  if (i % 1000 == 0 && i)
2733  {
2734  outfile << std::endl;
2735  }
2736  }
2737  outfile << std::endl;
2738  }
2739  }
2740 
2741  void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile,
2742  int expansion)
2743  {
2744  int i,j,k,l;
2745  int nbase = (*m_exp)[0]->GetNumBases();
2746  int cnt = 0;
2747 
2748  std::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
2749 
2750  if (expansion != -1)
2751  {
2752  exp = std::shared_ptr<LocalRegions::ExpansionVector>(
2753  new LocalRegions::ExpansionVector(1));
2754  (*exp)[0] = (*m_exp)[expansion];
2755  }
2756 
2757  if (nbase == 2)
2758  {
2759  for(i = 0; i < (*exp).size(); ++i)
2760  {
2761  const int np0 = (*exp)[i]->GetNumPoints(0);
2762  const int np1 = (*exp)[i]->GetNumPoints(1);
2763 
2764  for(j = 1; j < np1; ++j)
2765  {
2766  for(k = 1; k < np0; ++k)
2767  {
2768  outfile << cnt + (j-1)*np0 + k << " ";
2769  outfile << cnt + (j-1)*np0 + k+1 << " ";
2770  outfile << cnt + j *np0 + k+1 << " ";
2771  outfile << cnt + j *np0 + k << endl;
2772  }
2773  }
2774 
2775  cnt += np0*np1;
2776  }
2777  }
2778  else if (nbase == 3)
2779  {
2780  for(i = 0; i < (*exp).size(); ++i)
2781  {
2782  const int np0 = (*exp)[i]->GetNumPoints(0);
2783  const int np1 = (*exp)[i]->GetNumPoints(1);
2784  const int np2 = (*exp)[i]->GetNumPoints(2);
2785  const int np01 = np0*np1;
2786 
2787  for(j = 1; j < np2; ++j)
2788  {
2789  for(k = 1; k < np1; ++k)
2790  {
2791  for(l = 1; l < np0; ++l)
2792  {
2793  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l << " ";
2794  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l+1 << " ";
2795  outfile << cnt + (j-1)*np01 + k *np0 + l+1 << " ";
2796  outfile << cnt + (j-1)*np01 + k *np0 + l << " ";
2797  outfile << cnt + j *np01 + (k-1)*np0 + l << " ";
2798  outfile << cnt + j *np01 + (k-1)*np0 + l+1 << " ";
2799  outfile << cnt + j *np01 + k *np0 + l+1 << " ";
2800  outfile << cnt + j *np01 + k *np0 + l << endl;
2801  }
2802  }
2803  }
2804  cnt += np0*np1*np2;
2805  }
2806  }
2807  else
2808  {
2809  NEKERROR(ErrorUtil::efatal,"Not set up for this dimension");
2810  }
2811  }
2812 
2813  /**
2814  * Write Tecplot Files Field
2815  * @param outfile Output file name.
2816  * @param expansion Expansion that is considered
2817  */
2818  void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
2819  {
2820  if (expansion == -1)
2821  {
2822  int totpoints = GetTotPoints();
2823  if(m_physState == false)
2824  {
2825  BwdTrans(m_coeffs,m_phys);
2826  }
2827 
2828  for(int i = 0; i < totpoints; ++i)
2829  {
2830  outfile << m_phys[i] << " ";
2831  if(i % 1000 == 0 && i)
2832  {
2833  outfile << std::endl;
2834  }
2835  }
2836  outfile << std::endl;
2837 
2838  }
2839  else
2840  {
2841  int nPoints = (*m_exp)[expansion]->GetTotPoints();
2842 
2843  for (int i = 0; i < nPoints; ++i)
2844  {
2845  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
2846  }
2847 
2848  outfile << std::endl;
2849  }
2850  }
2851 
2852  void ExpList::WriteVtkHeader(std::ostream &outfile)
2853  {
2854  outfile << "<?xml version=\"1.0\"?>" << endl;
2855  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
2856  << "byte_order=\"LittleEndian\">" << endl;
2857  outfile << " <UnstructuredGrid>" << endl;
2858  }
2859 
2860  void ExpList::WriteVtkFooter(std::ostream &outfile)
2861  {
2862  outfile << " </UnstructuredGrid>" << endl;
2863  outfile << "</VTKFile>" << endl;
2864  }
2865 
2866  void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
2867  {
2868  boost::ignore_unused(istrip);
2869  int i,j,k;
2870  int nbase = (*m_exp)[expansion]->GetNumBases();
2871  int ntot = (*m_exp)[expansion]->GetTotPoints();
2872  int nquad[3];
2873 
2874  int ntotminus = 1;
2875  for(i = 0; i < nbase; ++i)
2876  {
2877  nquad[i] = (*m_exp)[expansion]->GetNumPoints(i);
2878  ntotminus *= (nquad[i]-1);
2879  }
2880 
2881  Array<OneD,NekDouble> coords[3];
2882  coords[0] = Array<OneD,NekDouble>(ntot, 0.0);
2883  coords[1] = Array<OneD,NekDouble>(ntot, 0.0);
2884  coords[2] = Array<OneD,NekDouble>(ntot, 0.0);
2885  (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
2886 
2887  outfile << " <Piece NumberOfPoints=\""
2888  << ntot << "\" NumberOfCells=\""
2889  << ntotminus << "\">" << endl;
2890  outfile << " <Points>" << endl;
2891  outfile << " <DataArray type=\"Float64\" "
2892  << "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
2893  outfile << " ";
2894  for (i = 0; i < ntot; ++i)
2895  {
2896  for (j = 0; j < 3; ++j)
2897  {
2898  outfile << setprecision(8) << scientific
2899  << (float)coords[j][i] << " ";
2900  }
2901  outfile << endl;
2902  }
2903  outfile << endl;
2904  outfile << " </DataArray>" << endl;
2905  outfile << " </Points>" << endl;
2906  outfile << " <Cells>" << endl;
2907  outfile << " <DataArray type=\"Int32\" "
2908  << "Name=\"connectivity\" format=\"ascii\">" << endl;
2909 
2910  int ns = 0; // pow(2,dim) for later usage
2911  string ostr;
2912  switch(m_expType)
2913  {
2914  case e1D:
2915  {
2916  ns = 2;
2917  ostr = "3 ";
2918  for (i = 0; i < nquad[0]-1; ++i)
2919  {
2920  outfile << i << " " << i+1 << endl;
2921  }
2922  }
2923  break;
2924  case e2D:
2925  {
2926  ns = 4;
2927  ostr = "9 ";
2928  for (i = 0; i < nquad[0]-1; ++i)
2929  {
2930  for (j = 0; j < nquad[1]-1; ++j)
2931  {
2932  outfile << j*nquad[0] + i << " "
2933  << j*nquad[0] + i + 1 << " "
2934  << (j+1)*nquad[0] + i + 1 << " "
2935  << (j+1)*nquad[0] + i << endl;
2936  }
2937  }
2938  }
2939  break;
2940  case e3D:
2941  {
2942  ns = 8;
2943  ostr = "12 ";
2944  for (i = 0; i < nquad[0]-1; ++i)
2945  {
2946  for (j = 0; j < nquad[1]-1; ++j)
2947  {
2948  for (k = 0; k < nquad[2]-1; ++k)
2949  {
2950  outfile << k*nquad[0]*nquad[1] + j*nquad[0] + i << " "
2951  << k*nquad[0]*nquad[1] + j*nquad[0] + i + 1 << " "
2952  << k*nquad[0]*nquad[1] + (j+1)*nquad[0] + i + 1 << " "
2953  << k*nquad[0]*nquad[1] + (j+1)*nquad[0] + i << " "
2954  << (k+1)*nquad[0]*nquad[1] + j*nquad[0] + i << " "
2955  << (k+1)*nquad[0]*nquad[1] + j*nquad[0] + i + 1 << " "
2956  << (k+1)*nquad[0]*nquad[1] + (j+1)*nquad[0] + i + 1 << " "
2957  << (k+1)*nquad[0]*nquad[1] + (j+1)*nquad[0] + i << " " << endl;
2958  }
2959  }
2960  }
2961  }
2962  break;
2963  default:
2964  break;
2965  }
2966 
2967 
2968  outfile << endl;
2969  outfile << " </DataArray>" << endl;
2970  outfile << " <DataArray type=\"Int32\" "
2971  << "Name=\"offsets\" format=\"ascii\">" << endl;
2972  for (i = 0; i < ntotminus; ++i)
2973  {
2974  outfile << i*ns+ns << " ";
2975  }
2976  outfile << endl;
2977  outfile << " </DataArray>" << endl;
2978  outfile << " <DataArray type=\"UInt8\" "
2979  << "Name=\"types\" format=\"ascii\">" << endl;
2980  for (i = 0; i < ntotminus; ++i)
2981  {
2982  outfile << ostr;
2983  }
2984  outfile << endl;
2985  outfile << " </DataArray>" << endl;
2986  outfile << " </Cells>" << endl;
2987  outfile << " <PointData>" << endl;
2988 
2989  }
2990 
2991  void ExpList::WriteVtkPieceFooter(std::ostream &outfile, int expansion)
2992  {
2993  boost::ignore_unused(expansion);
2994  outfile << " </PointData>" << endl;
2995  outfile << " </Piece>" << endl;
2996  }
2997 
2998  void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
2999  std::string var)
3000  {
3001  int i;
3002  int nq = (*m_exp)[expansion]->GetTotPoints();
3003 
3004  // printing the fields of that zone
3005  outfile << " <DataArray type=\"Float64\" Name=\""
3006  << var << "\">" << endl;
3007  outfile << " ";
3008  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
3009  for(i = 0; i < nq; ++i)
3010  {
3011  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
3012  }
3013  outfile << endl;
3014  outfile << " </DataArray>" << endl;
3015  }
3016 
3017  /**
3018  * Given a spectral/hp approximation
3019  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
3020  * (which should be contained in #m_phys), this function calculates the
3021  * \f$L_\infty\f$ error of this approximation with respect to an exact
3022  * solution. The local distribution of the quadrature points allows an
3023  * elemental evaluation of this operation through the functions
3024  * StdRegions#StdExpansion#Linf.
3025  *
3026  * The exact solution, also evaluated at the quadrature
3027  * points, should be contained in the variable #m_phys of
3028  * the ExpList object \a Sol.
3029  *
3030  * @param soln A 1D array, containing the discrete
3031  * evaluation of the exact solution at the
3032  * quadrature points in its array #m_phys.
3033  * @return The \f$L_\infty\f$ error of the approximation.
3034  */
3035  NekDouble ExpList::Linf(
3036  const Array<OneD, const NekDouble> &inarray,
3037  const Array<OneD, const NekDouble> &soln)
3038  {
3039  NekDouble err = 0.0;
3040 
3041  if (soln == NullNekDouble1DArray)
3042  {
3043  err = Vmath::Vmax(m_npoints, inarray, 1);
3044  }
3045  else
3046  {
3047  for (int i = 0; i < m_npoints; ++i)
3048  {
3049  err = max(err, abs(inarray[i] - soln[i]));
3050  }
3051  }
3052 
3053  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
3054 
3055  return err;
3056  }
3057 
3058  /**
3059  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
3060  * evaluated at the quadrature points (which should be contained in
3061  * #m_phys), this function calculates the \f$L_2\f$ error of this
3062  * approximation with respect to an exact solution. The local
3063  * distribution of the quadrature points allows an elemental evaluation
3064  * of this operation through the functions StdRegions#StdExpansion#L2.
3065  *
3066  * The exact solution, also evaluated at the quadrature points, should
3067  * be contained in the variable #m_phys of the ExpList object \a Sol.
3068  *
3069  * @param Sol An ExpList, containing the discrete
3070  * evaluation of the exact solution at the
3071  * quadrature points in its array #m_phys.
3072  * @return The \f$L_2\f$ error of the approximation.
3073  */
3074  NekDouble ExpList::v_L2(
3075  const Array<OneD, const NekDouble> &inarray,
3076  const Array<OneD, const NekDouble> &soln)
3077  {
3078  NekDouble err = 0.0, errl2;
3079  int i;
3080 
3081  if (soln == NullNekDouble1DArray)
3082  {
3083  for (i = 0; i < (*m_exp).size(); ++i)
3084  {
3085  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
3086  err += errl2*errl2;
3087  }
3088  }
3089  else
3090  {
3091  for (i = 0; i < (*m_exp).size(); ++i)
3092  {
3093  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
3094  soln + m_phys_offset[i]);
3095  err += errl2*errl2;
3096  }
3097  }
3098 
3099  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
3100 
3101  return sqrt(err);
3102  }
3103 
3104  /**
3105  * The integration is evaluated locally, that is
3106  * \f[\int
3107  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
3108  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
3109  * where the integration over the separate elements is done by the
3110  * function StdRegions#StdExpansion#Integral, which discretely
3111  * evaluates the integral using Gaussian quadrature.
3112  *
3113  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
3114  * containing the values of the function
3115  * \f$f(\boldsymbol{x})\f$ at the quadrature
3116  * points \f$\boldsymbol{x}_i\f$.
3117  * @return The value of the discretely evaluated integral
3118  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
3119  */
3120  NekDouble ExpList::v_Integral(const Array<OneD, const NekDouble> &inarray)
3121  {
3122  NekDouble sum = 0.0;
3123  int i = 0;
3124 
3125  for (i = 0; i < (*m_exp).size(); ++i)
3126  {
3127  sum += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
3128  }
3129  m_comm->GetRowComm()->AllReduce(sum, LibUtilities::ReduceSum);
3130 
3131  return sum;
3132  }
3133 
3134  NekDouble ExpList::v_VectorFlux(const Array<OneD, Array<OneD, NekDouble> > &inarray)
3135  {
3136  NekDouble flux = 0.0;
3137  int i = 0;
3138  int j;
3139 
3140  for (i = 0; i < (*m_exp).size(); ++i)
3141  {
3142  Array<OneD, Array<OneD, NekDouble> > tmp(inarray.size());
3143  for (j = 0; j < inarray.size(); ++j)
3144  {
3145  tmp[j] = Array<OneD, NekDouble>(inarray[j] + m_phys_offset[i]);
3146  }
3147  flux += (*m_exp)[i]->VectorFlux(tmp);
3148  }
3149 
3150  return flux;
3151  }
3152 
3153  Array<OneD, const NekDouble> ExpList::v_HomogeneousEnergy (void)
3154  {
3155  NEKERROR(ErrorUtil::efatal,
3156  "This method is not defined or valid for this class type");
3157  Array<OneD, NekDouble> NoEnergy(1,0.0);
3158  return NoEnergy;
3159  }
3160 
3161  LibUtilities::TranspositionSharedPtr ExpList::v_GetTransposition(void)
3162  {
3163  NEKERROR(ErrorUtil::efatal,
3164  "This method is not defined or valid for this class type");
3165  LibUtilities::TranspositionSharedPtr trans;
3166  return trans;
3167  }
3168 
3169  NekDouble ExpList::v_GetHomoLen(void)
3170  {
3171  NEKERROR(ErrorUtil::efatal,
3172  "This method is not defined or valid for this class type");
3173  NekDouble len = 0.0;
3174  return len;
3175  }
3176 
3177  void ExpList::v_SetHomoLen(const NekDouble lhom)
3178  {
3179  boost::ignore_unused(lhom);
3180  NEKERROR(ErrorUtil::efatal,
3181  "This method is not defined or valid for this class type");
3182  }
3183 
3184  Array<OneD, const unsigned int> ExpList::v_GetZIDs(void)
3185  {
3186  NEKERROR(ErrorUtil::efatal,
3187  "This method is not defined or valid for this class type");
3188  Array<OneD, unsigned int> NoModes(1);
3189  return NoModes;
3190  }
3191 
3192  Array<OneD, const unsigned int> ExpList::v_GetYIDs(void)
3193  {
3194  NEKERROR(ErrorUtil::efatal,
3195  "This method is not defined or valid for this class type");
3196  Array<OneD, unsigned int> NoModes(1);
3197  return NoModes;
3198  }
3199 
3200 
3201  void ExpList::v_ClearGlobalLinSysManager(void)
3202  {
3203  NEKERROR(ErrorUtil::efatal,
3204  "This method is not defined or valid for this class type");
3205  }
3206 
3207  void ExpList::ExtractFileBCs(
3208  const std::string &fileName,
3209  LibUtilities::CommSharedPtr comm,
3210  const std::string &varName,
3211  const std::shared_ptr<ExpList> locExpList)
3212  {
3213  string varString = fileName.substr(0, fileName.find_last_of("."));
3214  int j, k, len = varString.length();
3215  varString = varString.substr(len-1, len);
3216 
3217  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
3218  std::vector<std::vector<NekDouble> > FieldData;
3219 
3220  std::string ft = LibUtilities::FieldIO::GetFileType(fileName, comm);
3221  LibUtilities::FieldIOSharedPtr f = LibUtilities::GetFieldIOFactory()
3222  .CreateInstance(ft, comm, m_session->GetSharedFilesystem());
3223 
3224  f->Import(fileName, FieldDef, FieldData);
3225 
3226  bool found = false;
3227  for (j = 0; j < FieldDef.size(); ++j)
3228  {
3229  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
3230  {
3231  if (FieldDef[j]->m_fields[k] == varName)
3232  {
3233  // Copy FieldData into locExpList
3234  locExpList->ExtractDataToCoeffs(
3235  FieldDef[j], FieldData[j],
3236  FieldDef[j]->m_fields[k],
3237  locExpList->UpdateCoeffs());
3238  found = true;
3239  }
3240  }
3241  }
3242 
3243  ASSERTL0(found, "Could not find variable '"+varName+
3244  "' in file boundary condition "+fileName);
3245  locExpList->BwdTrans_IterPerExp(
3246  locExpList->GetCoeffs(),
3247  locExpList->UpdatePhys());
3248  }
3249 
3250  /**
3251  * Given a spectral/hp approximation
3252  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
3253  * (which should be contained in #m_phys), this function calculates the
3254  * \f$H^1_2\f$ error of this approximation with respect to an exact
3255  * solution. The local distribution of the quadrature points allows an
3256  * elemental evaluation of this operation through the functions
3257  * StdRegions#StdExpansion#H1.
3258  *
3259  * The exact solution, also evaluated at the quadrature points, should
3260  * be contained in the variable #m_phys of the ExpList object \a Sol.
3261  *
3262  * @param soln An 1D array, containing the discrete evaluation
3263  * of the exact solution at the quadrature points.
3264  *
3265  * @return The \f$H^1_2\f$ error of the approximation.
3266  */
3267  NekDouble ExpList::H1(
3268  const Array<OneD, const NekDouble> &inarray,
3269  const Array<OneD, const NekDouble> &soln)
3270  {
3271  NekDouble err = 0.0, errh1;
3272  int i;
3273 
3274  for (i = 0; i < (*m_exp).size(); ++i)
3275  {
3276  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
3277  soln + m_phys_offset[i]);
3278  err += errh1*errh1;
3279  }
3280 
3281  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
3282 
3283  return sqrt(err);
3284  }
3285 
3286  void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
3287  int NumHomoDir,
3288  Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
3289  std::vector<NekDouble> &HomoLen,
3290  bool homoStrips,
3291  std::vector<unsigned int> &HomoSIDs,
3292  std::vector<unsigned int> &HomoZIDs,
3293  std::vector<unsigned int> &HomoYIDs)
3294  {
3295  int startenum = (int) LibUtilities::eSegment;
3296  int endenum = (int) LibUtilities::eHexahedron;
3297  int s = 0;
3298  LibUtilities::ShapeType shape;
3299 
3300  ASSERTL1(NumHomoDir == HomoBasis.size(),"Homogeneous basis is not the same length as NumHomoDir");
3301  ASSERTL1(NumHomoDir == HomoLen.size(),"Homogeneous length vector is not the same length as NumHomDir");
3302 
3303  // count number of shapes
3304  switch((*m_exp)[0]->GetShapeDimension())
3305  {
3306  case 1:
3307  startenum = (int) LibUtilities::eSegment;
3308  endenum = (int) LibUtilities::eSegment;
3309  break;
3310  case 2:
3311  startenum = (int) LibUtilities::eTriangle;
3312  endenum = (int) LibUtilities::eQuadrilateral;
3313  break;
3314  case 3:
3315  startenum = (int) LibUtilities::eTetrahedron;
3316  endenum = (int) LibUtilities::eHexahedron;
3317  break;
3318  }
3319 
3320  for(s = startenum; s <= endenum; ++s)
3321  {
3322  std::vector<unsigned int> elementIDs;
3323  std::vector<LibUtilities::BasisType> basis;
3324  std::vector<unsigned int> numModes;
3325  std::vector<std::string> fields;
3326 
3327  bool first = true;
3328  bool UniOrder = true;
3329  int n;
3330 
3331  shape = (LibUtilities::ShapeType) s;
3332 
3333  for(int i = 0; i < (*m_exp).size(); ++i)
3334  {
3335  if((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
3336  {
3337  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
3338  if(first)
3339  {
3340  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3341  {
3342  basis.push_back((*m_exp)[i]->GetBasis(j)->GetBasisType());
3343  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
3344  }
3345 
3346  // add homogeneous direction details if defined
3347  for(n = 0 ; n < NumHomoDir; ++n)
3348  {
3349  basis.push_back(HomoBasis[n]->GetBasisType());
3350  numModes.push_back(HomoBasis[n]->GetNumModes());
3351  }
3352 
3353  first = false;
3354  }
3355  else
3356  {
3357  ASSERTL0((*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],"Routine is not set up for multiple bases definitions");
3358 
3359  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3360  {
3361  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
3362  if(numModes[j] != (*m_exp)[i]->GetBasis(j)->GetNumModes())
3363  {
3364  UniOrder = false;
3365  }
3366  }
3367  // add homogeneous direction details if defined
3368  for(n = 0 ; n < NumHomoDir; ++n)
3369  {
3370  numModes.push_back(HomoBasis[n]->GetNumModes());
3371  }
3372  }
3373  }
3374  }
3375 
3376 
3377  if(elementIDs.size() > 0)
3378  {
3379  LibUtilities::FieldDefinitionsSharedPtr fdef =
3380  MemoryManager<LibUtilities::FieldDefinitions>::
3381  AllocateSharedPtr(shape, elementIDs, basis,
3382  UniOrder, numModes,fields,
3383  NumHomoDir, HomoLen, homoStrips,
3384  HomoSIDs, HomoZIDs, HomoYIDs);
3385  fielddef.push_back(fdef);
3386  }
3387  }
3388  }
3389 
3390  //
3391  // Virtual functions
3392  //
3393  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::v_GetFieldDefinitions()
3394  {
3395  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
3396  v_GetFieldDefinitions(returnval);
3397  return returnval;
3398  }
3399 
3400  void ExpList::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
3401  {
3402  GeneralGetFieldDefinitions(fielddef);
3403  }
3404 
3405  //Append the element data listed in elements
3406  //fielddef->m_ElementIDs onto fielddata
3407  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
3408  {
3409  v_AppendFieldData(fielddef,fielddata,m_coeffs);
3410  }
3411 
3412  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
3413  {
3414  int i;
3415  // Determine mapping from element ids to location in
3416  // expansion list
3417  // Determine mapping from element ids to location in
3418  // expansion list
3419  map<int, int> ElmtID_to_ExpID;
3420  for(i = 0; i < (*m_exp).size(); ++i)
3421  {
3422  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3423  }
3424 
3425  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
3426  {
3427  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
3428  int datalen = (*m_exp)[eid]->GetNcoeffs();
3429  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]],&coeffs[m_coeff_offset[eid]]+datalen);
3430  }
3431 
3432  }
3433 
3434  /// Extract the data in fielddata into the coeffs
3435  void ExpList::ExtractDataToCoeffs(
3436  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
3437  std::vector<NekDouble> &fielddata,
3438  std::string &field,
3439  Array<OneD, NekDouble> &coeffs)
3440  {
3441  v_ExtractDataToCoeffs(fielddef,fielddata,field,coeffs);
3442  }
3443 
3444  void ExpList::ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
3445  {
3446  v_ExtractCoeffsToCoeffs(fromExpList,fromCoeffs,toCoeffs);
3447  }
3448 
3449  /**
3450  * @brief Extract data from raw field data into expansion list.
3451  *
3452  * @param fielddef Field definitions.
3453  * @param fielddata Data for associated field.
3454  * @param field Field variable name.
3455  * @param coeffs Resulting coefficient array.
3456  */
3457  void ExpList::v_ExtractDataToCoeffs(
3458  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
3459  std::vector<NekDouble> &fielddata,
3460  std::string &field,
3461  Array<OneD, NekDouble> &coeffs)
3462  {
3463  int i, expId;
3464  int offset = 0;
3465  int modes_offset = 0;
3466  int datalen = fielddata.size()/fielddef->m_fields.size();
3467 
3468  // Find data location according to field definition
3469  for(i = 0; i < fielddef->m_fields.size(); ++i)
3470  {
3471  if(fielddef->m_fields[i] == field)
3472  {
3473  break;
3474  }
3475  offset += datalen;
3476  }
3477 
3478  ASSERTL0(i != fielddef->m_fields.size(),
3479  "Field (" + field + ") not found in file.");
3480 
3481  if (m_elmtToExpId.size() == 0)
3482  {
3483  // Loop in reverse order so that in case where using a
3484  // Homogeneous expansion it sets geometry ids to first part of
3485  // m_exp list. Otherwise will set to second (complex) expansion
3486  for(i = (*m_exp).size()-1; i >= 0; --i)
3487  {
3488  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3489  }
3490  }
3491 
3492  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
3493  {
3494  // Reset modes_offset in the case where all expansions of
3495  // the same order.
3496  if (fielddef->m_uniOrder == true)
3497  {
3498  modes_offset = 0;
3499  }
3500 
3501  datalen = LibUtilities::GetNumberOfCoefficients(
3502  fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
3503 
3504  const int elmtId = fielddef->m_elementIDs[i];
3505  auto eIt = m_elmtToExpId.find(elmtId);
3506 
3507  if (eIt == m_elmtToExpId.end())
3508  {
3509  offset += datalen;
3510  modes_offset += (*m_exp)[0]->GetNumBases();
3511  continue;
3512  }
3513 
3514  expId = eIt->second;
3515 
3516  bool sameBasis = true;
3517  for (int j = 0; j < fielddef->m_basis.size(); ++j)
3518  {
3519  if (fielddef->m_basis[j] != (*m_exp)[expId]->GetBasisType(j))
3520  {
3521  sameBasis = false;
3522  break;
3523  }
3524  }
3525 
3526  if (datalen == (*m_exp)[expId]->GetNcoeffs() && sameBasis)
3527  {
3528  Vmath::Vcopy(datalen, &fielddata[offset], 1,
3529  &coeffs[m_coeff_offset[expId]], 1);
3530  }
3531  else
3532  {
3533  (*m_exp)[expId]->ExtractDataToCoeffs(
3534  &fielddata[offset], fielddef->m_numModes,
3535  modes_offset, &coeffs[m_coeff_offset[expId]],
3536  fielddef->m_basis);
3537  }
3538 
3539  offset += datalen;
3540  modes_offset += (*m_exp)[0]->GetNumBases();
3541  }
3542 
3543  return;
3544  }
3545 
3546  void ExpList::v_ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
3547  {
3548  int i;
3549  int offset = 0;
3550 
3551  for(i = 0; i < (*m_exp).size(); ++i)
3552  {
3553  std::vector<unsigned int> nummodes;
3554  vector<LibUtilities::BasisType> basisTypes;
3555  for(int j= 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
3556  {
3557  nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
3558  basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
3559  }
3560 
3561  (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
3562  &toCoeffs[m_coeff_offset[i]],
3563  basisTypes);
3564 
3565  offset += fromExpList->GetExp(i)->GetNcoeffs();
3566  }
3567  }
3568 
3569  /**
3570  * Get the weight value on boundaries
3571  */
3572  void ExpList::GetBwdWeight(
3573  Array<OneD, NekDouble> &weightAver,
3574  Array<OneD, NekDouble> &weightJump)
3575  {
3576  size_t nTracePts = weightAver.size();
3577  // average for interior traces
3578  for(int i = 0; i < nTracePts; ++i)
3579  {
3580  weightAver[i] = 0.5;
3581  weightJump[i] = 1.0;
3582  }
3583  FillBwdWithBwdWeight(weightAver, weightJump);
3584  }
3585 
3586  void ExpList::v_GetMovingFrames(
3587  const SpatialDomains::GeomMMF MMFdir,
3588  const Array<OneD, const NekDouble> &CircCentre,
3589  Array<OneD, Array<OneD, NekDouble> > &outarray)
3590  {
3591  int npts;
3592 
3593  int MFdim = 3;
3594  int nq = outarray[0].size()/MFdim;
3595 
3596  // Assume whole array is of same coordinate dimension
3597  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
3598 
3599  Array<OneD, Array<OneD, NekDouble> > MFloc(MFdim*coordim);
3600  // Process each expansion.
3601  for(int i = 0; i < m_exp->size(); ++i)
3602  {
3603  npts = (*m_exp)[i]->GetTotPoints();
3604 
3605  for (int j=0; j< MFdim*coordim; ++j)
3606  {
3607  MFloc[j] = Array<OneD, NekDouble>(npts, 0.0);
3608  }
3609 
3610  // MF from LOCALREGIONS
3611  (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
3612  (*m_exp)[i]->GetPointsKeys(),
3613  MMFdir,
3614  CircCentre,
3615  MFloc );
3616 
3617  // Get the physical data offset for this expansion.
3618  for (int j = 0; j < MFdim; ++j)
3619  {
3620  for (int k = 0; k < coordim; ++k)
3621  {
3622  Vmath::Vcopy(npts,
3623  &MFloc[j*coordim+k][0], 1,
3624  &outarray[j][k*nq+m_phys_offset[i]], 1);
3625  }
3626  }
3627  }
3628 
3629  }
3630 
3631 
3632 
3633  /**
3634  * @brief Generate vector v such that v[i] = scalar1 if i is in the
3635  * element < ElementID. Otherwise, v[i] = scalar2.
3636  *
3637  */
3638  void ExpList::GenerateElementVector(
3639  const int ElementID,
3640  const NekDouble scalar1,
3641  const NekDouble scalar2,
3642  Array<OneD, NekDouble> &outarray)
3643  {
3644  int npoints_e;
3645  NekDouble coeff;
3646 
3647  Array<OneD, NekDouble> outarray_e;
3648 
3649  for(int i = 0 ; i < (*m_exp).size(); ++i)
3650  {
3651  npoints_e = (*m_exp)[i]->GetTotPoints();
3652 
3653  if(i <= ElementID)
3654  {
3655  coeff = scalar1;
3656  }
3657  else
3658  {
3659  coeff = scalar2;
3660  }
3661 
3662  outarray_e = Array<OneD, NekDouble>(npoints_e, coeff);
3663  Vmath::Vcopy(npoints_e, &outarray_e[0], 1,
3664  &outarray[m_phys_offset[i]], 1);
3665  }
3666  }
3667 
3668  const Array<OneD,const std::shared_ptr<ExpList> >
3669  &ExpList::v_GetBndCondExpansions(void)
3670  {
3671  NEKERROR(ErrorUtil::efatal,
3672  "This method is not defined or valid for this class type");
3673  static Array<OneD,const std::shared_ptr<ExpList> > result;
3674  return result;
3675  }
3676 
3677  std::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
3678  {
3679  boost::ignore_unused(i);
3680  NEKERROR(ErrorUtil::efatal,
3681  "This method is not defined or valid for this class type");
3682  static std::shared_ptr<ExpList> result;
3683  return result;
3684  }
3685 
3686  /**
3687  * Upwind the left and right states given by the Arrays Fwd and Bwd
3688  * using the vector quantity Vec and ouput the upwinded value in the
3689  * array upwind.
3690  *
3691  * @param Vec Velocity field.
3692  * @param Fwd Left state.
3693  * @param Bwd Right state.
3694  * @param Upwind Output vector.
3695  */
3696  void ExpList::v_Upwind(
3697  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
3698  const Array<OneD, const NekDouble> &Fwd,
3699  const Array<OneD, const NekDouble> &Bwd,
3700  Array<OneD, NekDouble> &Upwind)
3701  {
3702  switch(m_expType)
3703  {
3704  case e1D:
3705  {
3706  int i,j,k,e_npoints,offset;
3707  Array<OneD,NekDouble> normals;
3708  NekDouble Vn;
3709 
3710  // Assume whole array is of same coordimate dimension
3711  int coordim = GetCoordim(0);
3712 
3713  ASSERTL1(Vec.size() >= coordim,
3714  "Input vector does not have sufficient dimensions to "
3715  "match coordim");
3716 
3717  // Process each expansion
3718  for(i = 0; i < m_exp->size(); ++i)
3719  {
3720  // Get the number of points in the expansion and the normals.
3721  e_npoints = (*m_exp)[i]->GetNumPoints(0);
3722  normals = (*m_exp)[i]->GetPhysNormals();
3723 
3724  // Get the physical data offset of the expansion in m_phys.
3725  offset = m_phys_offset[i];
3726 
3727  // Compute each data point.
3728  for(j = 0; j < e_npoints; ++j)
3729  {
3730  // Calculate normal velocity.
3731  Vn = 0.0;
3732  for(k = 0; k < coordim; ++k)
3733  {
3734  Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
3735  }
3736 
3737  // Upwind based on direction of normal velocity.
3738  if(Vn > 0.0)
3739  {
3740  Upwind[offset + j] = Fwd[offset + j];
3741  }
3742  else
3743  {
3744  Upwind[offset + j] = Bwd[offset + j];
3745  }
3746  }
3747  }
3748  }
3749  break;
3750  default:
3751  NEKERROR(ErrorUtil::efatal,
3752  "This method is not defined or valid for this class type");
3753  break;
3754  }
3755  }
3756 
3757  /**
3758  * One-dimensional upwind.
3759  * \see ExpList::Upwind(
3760  * const Array<OneD, const Array<OneD, NekDouble> >,
3761  * const Array<OneD, const NekDouble>,
3762  * const Array<OneD, const NekDouble>,
3763  * Array<OneD, NekDouble>, int)
3764  *
3765  * @param Vn Velocity field.
3766  * @param Fwd Left state.
3767  * @param Bwd Right state.
3768  * @param Upwind Output vector.
3769  */
3770  void ExpList::v_Upwind(
3771  const Array<OneD, const NekDouble> &Vn,
3772  const Array<OneD, const NekDouble> &Fwd,
3773  const Array<OneD, const NekDouble> &Bwd,
3774  Array<OneD, NekDouble> &Upwind)
3775  {
3776  ASSERTL1(Vn.size() >= m_npoints,"Vn is not of sufficient length");
3777  ASSERTL1(Fwd.size() >= m_npoints,"Fwd is not of sufficient length");
3778  ASSERTL1(Bwd.size() >= m_npoints,"Bwd is not of sufficient length");
3779  ASSERTL1(Upwind.size() >= m_npoints,
3780  "Upwind is not of sufficient length");
3781 
3782  // Process each point in the expansion.
3783  for(int j = 0; j < m_npoints; ++j)
3784  {
3785  // Upwind based on one-dimensional velocity.
3786  if(Vn[j] > 0.0)
3787  {
3788  Upwind[j] = Fwd[j];
3789  }
3790  else
3791  {
3792  Upwind[j] = Bwd[j];
3793  }
3794  }
3795  }
3796 
3797  std::shared_ptr<ExpList> &ExpList::v_GetTrace()
3798  {
3799  NEKERROR(ErrorUtil::efatal,
3800  "This method is not defined or valid for this class type");
3801  static std::shared_ptr<ExpList> returnVal;
3802  return returnVal;
3803  }
3804 
3805  std::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
3806  {
3807  NEKERROR(ErrorUtil::efatal,
3808  "This method is not defined or valid for this class type");
3809  static std::shared_ptr<AssemblyMapDG> result;
3810  return result;
3811  }
3812 
3813  const Array<OneD, const int> &ExpList::v_GetTraceBndMap()
3814  {
3815  return GetTraceMap()->GetBndCondIDToGlobalTraceID();
3816  }
3817 
3818  /**
3819  * @brief Helper function to re-align face to a given orientation.
3820  */
3821  void AlignFace(const StdRegions::Orientation orient,
3822  const int nquad1,
3823  const int nquad2,
3824  const Array<OneD, const NekDouble> &in,
3825  Array<OneD, NekDouble> &out)
3826  {
3827  // Copy transpose.
3828  if (orient == StdRegions::eDir1FwdDir2_Dir2FwdDir1 ||
3829  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
3830  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
3831  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
3832  {
3833  for (int i = 0; i < nquad2; ++i)
3834  {
3835  for (int j = 0; j < nquad1; ++j)
3836  {
3837  out[i*nquad1 + j] = in[j*nquad2 + i];
3838  }
3839  }
3840  }
3841  else
3842  {
3843  for (int i = 0; i < nquad2; ++i)
3844  {
3845  for (int j = 0; j < nquad1; ++j)
3846  {
3847  out[i*nquad1 + j] = in[i*nquad1 + j];
3848  }
3849  }
3850  }
3851 
3852  if (orient == StdRegions::eDir1BwdDir1_Dir2FwdDir2 ||
3853  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
3854  orient == StdRegions::eDir1BwdDir2_Dir2FwdDir1 ||
3855  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
3856  {
3857  // Reverse x direction
3858  for (int i = 0; i < nquad2; ++i)
3859  {
3860  for (int j = 0; j < nquad1/2; ++j)
3861  {
3862  swap(out[i*nquad1 + j],
3863  out[i*nquad1 + nquad1-j-1]);
3864  }
3865  }
3866  }
3867 
3868  if (orient == StdRegions::eDir1FwdDir1_Dir2BwdDir2 ||
3869  orient == StdRegions::eDir1BwdDir1_Dir2BwdDir2 ||
3870  orient == StdRegions::eDir1FwdDir2_Dir2BwdDir1 ||
3871  orient == StdRegions::eDir1BwdDir2_Dir2BwdDir1)
3872  {
3873  // Reverse y direction
3874  for (int j = 0; j < nquad1; ++j)
3875  {
3876  for (int i = 0; i < nquad2/2; ++i)
3877  {
3878  swap(out[i*nquad1 + j],
3879  out[(nquad2-i-1)*nquad1 + j]);
3880  }
3881  }
3882  }
3883  }
3884 
3885 
3886  /**
3887  * For each local element, copy the normals stored in the element list
3888  * into the array \a normals.
3889  * @param normals Multidimensional array in which to copy normals
3890  * to. Must have dimension equal to or larger than
3891  * the spatial dimension of the elements.
3892  */
3893  void ExpList::v_GetNormals(
3894  Array<OneD, Array<OneD, NekDouble> > &normals)
3895  {
3896  int i,j,k,e_npoints,offset;
3897  Array<OneD,Array<OneD,NekDouble> > locnormals;
3898 
3899  // Assume whole array is of same coordinate dimension
3900  int coordim = GetCoordim(0);
3901 
3902  ASSERTL1(normals.size() >= coordim,
3903  "Output vector does not have sufficient dimensions to "
3904  "match coordim");
3905 
3906  switch(m_expType)
3907  {
3908  case e0D:
3909  {
3910  // Process each expansion.
3911  for(i = 0; i < m_exp->size(); ++i)
3912  {
3913  LocalRegions::ExpansionSharedPtr loc_exp = (*m_exp)[i];
3914 
3915  LocalRegions::ExpansionSharedPtr loc_elmt =
3916  loc_exp->GetLeftAdjacentElementExp();
3917 
3918  // Get the number of points and normals for this expansion.
3919  e_npoints = 1;
3920  locnormals = loc_elmt->GetTraceNormal(loc_exp->
3921  GetLeftAdjacentElementTrace());
3922 
3923  // Get the physical data offset for this expansion.
3924  offset = m_phys_offset[i];
3925 
3926  // Process each point in the expansion.
3927  for(j = 0; j < e_npoints; ++j)
3928  {
3929  // Process each spatial dimension and copy the
3930  // values into the output array.
3931  for(k = 0; k < coordim; ++k)
3932  {
3933  normals[k][offset] = locnormals[k][0];
3934  }
3935  }
3936  }
3937  }
3938  break;
3939  case e1D:
3940  {
3941  SpatialDomains::Geometry1DSharedPtr segGeom;
3942  Array<OneD,Array<OneD,NekDouble> > locnormals2;
3943  Array<OneD,Array<OneD,NekDouble> > Norms;
3944 
3945  for (i = 0; i < m_exp->size(); ++i)
3946  {
3947  LocalRegions::ExpansionSharedPtr loc_exp =(*m_exp)[i];
3948 
3949  LocalRegions::ExpansionSharedPtr loc_elmt =
3950  loc_exp->GetLeftAdjacentElementExp();
3951 
3952  int edgeNumber = loc_exp->GetLeftAdjacentElementTrace();
3953 
3954  // Get the number of points and normals for this expansion.
3955  e_npoints = (*m_exp)[i]->GetNumPoints(0);
3956 
3957  locnormals = loc_elmt->GetTraceNormal(edgeNumber);
3958  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
3959  int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
3960 
3961  if (e_nmodes != loc_nmodes)
3962  {
3963  if (loc_exp->GetRightAdjacentElementTrace() >= 0)
3964  {
3965  LocalRegions::ExpansionSharedPtr loc_elmt =
3966  loc_exp->GetRightAdjacentElementExp();
3967 
3968  int EdgeNumber = loc_exp->
3969  GetRightAdjacentElementTrace();
3970 
3971  // Serial case: right element is connected so we can
3972  // just grab that normal.
3973  locnormals = loc_elmt->GetTraceNormal(EdgeNumber);
3974 
3975  offset = m_phys_offset[i];
3976 
3977  // Process each point in the expansion.
3978  for (j = 0; j < e_npoints; ++j)
3979  {
3980  // Process each spatial dimension and
3981  // copy the values into the output
3982  // array.
3983  for (k = 0; k < coordim; ++k)
3984  {
3985  normals[k][offset + j] = -locnormals[k][j];
3986  }
3987  }
3988  }
3989  else
3990  {
3991  // Parallel case: need to interpolate normal.
3992  Array<OneD, Array<OneD, NekDouble> > normal(coordim);
3993 
3994  for (int p = 0; p < coordim; ++p)
3995  {
3996  normal[p] = Array<OneD, NekDouble>(e_npoints,0.0);
3997  LibUtilities::PointsKey to_key =
3998  loc_exp->GetBasis(0)->GetPointsKey();
3999  LibUtilities::PointsKey from_key =
4000  loc_elmt->GetBasis(0)->GetPointsKey();
4001  LibUtilities::Interp1D(from_key,
4002  locnormals[p],
4003  to_key,
4004  normal[p]);
4005  }
4006 
4007  offset = m_phys_offset[i];
4008 
4009  // Process each point in the expansion.
4010  for (j = 0; j < e_npoints; ++j)
4011  {
4012  // Process each spatial dimension and copy the values
4013  // into the output array.
4014  for (k = 0; k < coordim; ++k)
4015  {
4016  normals[k][offset + j] = normal[k][j];
4017  }
4018  }
4019  }
4020  }
4021  else
4022  {
4023  // Get the physical data offset for this expansion.
4024  offset = m_phys_offset[i];
4025 
4026  // Process each point in the expansion.
4027  for (j = 0; j < e_npoints; ++j)
4028  {
4029  // Process each spatial dimension and copy the values
4030  // into the output array.
4031  for (k = 0; k < coordim; ++k)
4032  {
4033  normals[k][offset + j] = locnormals[k][j];
4034  }
4035  }
4036  }
4037  }
4038  }
4039  break;
4040  case e2D:
4041  {
4042  Array<OneD, NekDouble> tmp;
4043 
4044  // Process each expansion.
4045  for (i = 0; i < m_exp->size(); ++i)
4046  {
4047  LocalRegions::ExpansionSharedPtr traceExp = (*m_exp)[i];
4048  LocalRegions::ExpansionSharedPtr exp3D =
4049  traceExp->GetLeftAdjacentElementExp();
4050 
4051  // Get the number of points and normals for this expansion.
4052  int faceNum = traceExp->GetLeftAdjacentElementTrace();
4053  int offset = m_phys_offset[i];
4054 
4055  const Array<OneD, const Array<OneD, NekDouble> > &locNormals
4056  = exp3D->GetTraceNormal(faceNum);
4057 
4058  // Project normals from 3D element onto the same
4059  // orientation as the trace expansion.
4060  StdRegions::Orientation orient = exp3D->
4061  GetTraceOrient(faceNum);
4062 
4063  int fromid0,fromid1;
4064 
4065  if(orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
4066  {
4067  fromid0 = 0;
4068  fromid1 = 1;
4069  }
4070  else
4071  {
4072  fromid0 = 1;
4073  fromid1 = 0;
4074  }
4075 
4076  LibUtilities::BasisKey faceBasis0
4077  = exp3D->GetTraceBasisKey(faceNum, fromid0);
4078  LibUtilities::BasisKey faceBasis1
4079  = exp3D->GetTraceBasisKey(faceNum, fromid1);
4080  LibUtilities::BasisKey traceBasis0
4081  = traceExp->GetBasis(0)->GetBasisKey();
4082  LibUtilities::BasisKey traceBasis1
4083  = traceExp->GetBasis(1)->GetBasisKey();
4084 
4085  const int faceNq0 = faceBasis0.GetNumPoints();
4086  const int faceNq1 = faceBasis1.GetNumPoints();
4087 
4088  Array<OneD, int> faceids;
4089  exp3D->ReOrientTracePhysMap(orient,faceids,faceNq0,faceNq1);
4090 
4091  Array<OneD, NekDouble> traceNormals(faceNq0 * faceNq1);
4092 
4093  for (j = 0; j < coordim; ++j)
4094  {
4095  Vmath::Scatr(faceNq0*faceNq1,locNormals[j],faceids,
4096  traceNormals);
4097 
4098  LibUtilities::Interp2D(faceBasis0.GetPointsKey(),
4099  faceBasis1.GetPointsKey(),
4100  traceNormals,
4101  traceBasis0.GetPointsKey(),
4102  traceBasis1.GetPointsKey(),
4103  tmp = normals[j]+offset);
4104  }
4105  }
4106  }
4107  break;
4108  default:
4109  {
4110  NEKERROR(ErrorUtil::efatal,
4111  "This method is not defined or valid for this class type");
4112  }
4113  }
4114  }
4115 
4116  void ExpList::GetElmtNormalLength(
4117  Array<OneD, NekDouble> &lengthsFwd,
4118  Array<OneD, NekDouble> &lengthsBwd)
4119  {
4120  int e_npoints;
4121 
4122  Array<OneD, NekDouble> locLeng;
4123  Array<OneD, NekDouble> lengintp;
4124  Array<OneD, NekDouble> lengAdd;
4125  Array<OneD, int > LRbndnumbs(2);
4126  Array<OneD, Array<OneD,NekDouble> > lengLR(2);
4127  lengLR[0] = lengthsFwd;
4128  lengLR[1] = lengthsBwd;
4129  Array<OneD, LocalRegions::ExpansionSharedPtr> LRelmts(2);
4130  LocalRegions::ExpansionSharedPtr loc_elmt;
4131  LocalRegions::ExpansionSharedPtr loc_exp;
4132  int e_npoints0 = -1;
4133  if(m_expType == e1D)
4134  {
4135  for (int i = 0; i < m_exp->size(); ++i)
4136  {
4137  loc_exp = (*m_exp)[i];
4138  int offset = m_phys_offset[i];
4139 
4140  int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
4141  e_npoints = (*m_exp)[i]->GetNumPoints(0);
4142  if ( e_npoints0 < e_npoints)
4143  {
4144  lengintp = Array<OneD, NekDouble>{size_t(e_npoints),0.0};
4145  e_npoints0 = e_npoints;
4146  }
4147 
4148  LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4149  LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4150 
4151  LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4152  LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4153  for (int nlr = 0; nlr < 2; ++nlr)
4154  {
4155  Vmath::Zero(e_npoints0, lengintp, 1);
4156  lengAdd = lengintp;
4157  int bndNumber = LRbndnumbs[nlr];
4158  loc_elmt = LRelmts[nlr];
4159  if (bndNumber >= 0)
4160  {
4161  locLeng = loc_elmt->GetElmtBndNormDirElmtLen(
4162  bndNumber);
4163  lengAdd = locLeng;
4164 
4165  int loc_nmodes = loc_elmt->GetBasis(0)->
4166  GetNumModes();
4167  if (e_nmodes != loc_nmodes)
4168  {
4169  // Parallel case: need to interpolate.
4170  LibUtilities::PointsKey to_key =
4171  loc_exp->GetBasis(0)->GetPointsKey();
4172  LibUtilities::PointsKey from_key =
4173  loc_elmt->GetBasis(0)->GetPointsKey();
4174  LibUtilities::Interp1D(from_key, locLeng,
4175  to_key, lengintp);
4176  lengAdd = lengintp;
4177  }
4178  }
4179  for (int j = 0; j < e_npoints; ++j)
4180  {
4181  lengLR[nlr][offset + j] = lengAdd[j];
4182  }
4183  }
4184  }
4185  }
4186  else if (m_expType == e2D)
4187  {
4188  for (int i = 0; i < m_exp->size(); ++i)
4189  {
4190  loc_exp = (*m_exp)[i];
4191  int offset = m_phys_offset[i];
4192 
4193  LibUtilities::BasisKey traceBasis0
4194  = loc_exp->GetBasis(0)->GetBasisKey();
4195  LibUtilities::BasisKey traceBasis1
4196  = loc_exp->GetBasis(1)->GetBasisKey();
4197  const int TraceNq0 = traceBasis0.GetNumPoints();
4198  const int TraceNq1 = traceBasis1.GetNumPoints();
4199  e_npoints = TraceNq0*TraceNq1;
4200  if (e_npoints0 < e_npoints)
4201  {
4202  lengintp = Array<OneD,NekDouble>{size_t(e_npoints),
4203  0.0};
4204  e_npoints0 = e_npoints;
4205  }
4206 
4207  LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4208  LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4209 
4210  LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4211  LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4212  for (int nlr = 0; nlr < 2; ++nlr)
4213  {
4214  Vmath::Zero(e_npoints0, lengintp, 1);
4215  int bndNumber = LRbndnumbs[nlr];
4216  loc_elmt = LRelmts[nlr];
4217  if (bndNumber >= 0)
4218  {
4219  locLeng = loc_elmt->GetElmtBndNormDirElmtLen(
4220  bndNumber);
4221  // Project normals from 3D element onto the
4222  // same orientation as the trace expansion.
4223  StdRegions::Orientation orient = loc_elmt->
4224  GetTraceOrient(bndNumber);
4225 
4226  int fromid0,fromid1;
4227  if (orient < StdRegions::eDir1FwdDir2_Dir2FwdDir1)
4228  {
4229  fromid0 = 0;
4230  fromid1 = 1;
4231  }
4232  else
4233  {
4234  fromid0 = 1;
4235  fromid1 = 0;
4236  }
4237 
4238  LibUtilities::BasisKey faceBasis0
4239  = loc_elmt->GetTraceBasisKey(bndNumber, fromid0);
4240  LibUtilities::BasisKey faceBasis1
4241  = loc_elmt->GetTraceBasisKey(bndNumber, fromid1);
4242  const int faceNq0 = faceBasis0.GetNumPoints();
4243  const int faceNq1 = faceBasis1.GetNumPoints();
4244  Array<OneD, NekDouble> alignedLeng(faceNq0*faceNq1);
4245 
4246  AlignFace(orient, faceNq0, faceNq1,
4247  locLeng, alignedLeng);
4248  LibUtilities::Interp2D(faceBasis0.GetPointsKey(),
4249  faceBasis1.GetPointsKey(),
4250  alignedLeng,
4251  traceBasis0.GetPointsKey(),
4252  traceBasis1.GetPointsKey(),
4253  lengintp);
4254  }
4255  for (int j = 0; j < e_npoints; ++j)
4256  {
4257  lengLR[nlr][offset + j] = lengintp[j];
4258  }
4259  }
4260  }
4261  }
4262  }
4263 
4264  void ExpList::v_AddTraceIntegral(
4265  const Array<OneD, const NekDouble> &Fx,
4266  const Array<OneD, const NekDouble> &Fy,
4267  Array<OneD, NekDouble> &outarray)
4268  {
4269  boost::ignore_unused(Fx, Fy, outarray);
4270  NEKERROR(ErrorUtil::efatal,
4271  "This method is not defined or valid for this class type");
4272  }
4273 
4274  void ExpList::v_AddTraceIntegral(
4275  const Array<OneD, const NekDouble> &Fn,
4276  Array<OneD, NekDouble> &outarray)
4277  {
4278  boost::ignore_unused(Fn, outarray);
4279  NEKERROR(ErrorUtil::efatal,
4280  "This method is not defined or valid for this class type");
4281  }
4282 
4283  void ExpList::v_AddFwdBwdTraceIntegral(
4284  const Array<OneD, const NekDouble> &Fwd,
4285  const Array<OneD, const NekDouble> &Bwd,
4286  Array<OneD, NekDouble> &outarray)
4287  {
4288  boost::ignore_unused(Fwd, Bwd, outarray);
4289  NEKERROR(ErrorUtil::efatal,
4290  "This method is not defined or valid for this class type");
4291  }
4292 
4293  void ExpList::v_GetFwdBwdTracePhys(Array<OneD,NekDouble> &Fwd,
4294  Array<OneD,NekDouble> &Bwd)
4295  {
4296  boost::ignore_unused(Fwd, Bwd);
4297  NEKERROR(ErrorUtil::efatal,
4298  "This method is not defined or valid for this class type");
4299  }
4300 
4301  void ExpList::v_GetFwdBwdTracePhys(
4302  const Array<OneD,const NekDouble> &field,
4303  Array<OneD,NekDouble> &Fwd,
4304  Array<OneD,NekDouble> &Bwd,
4305  bool FillBnd,
4306  bool PutFwdInBwdOnBCs,
4307  bool DoExchange)
4308  {
4309  boost::ignore_unused(field, Fwd, Bwd, FillBnd,
4310  PutFwdInBwdOnBCs, DoExchange);
4311  NEKERROR(ErrorUtil::efatal,
4312  "This method is not defined or valid for this class type");
4313  }
4314 
4315  void ExpList::v_FillBwdWithBoundCond(
4316  const Array<OneD,NekDouble> &Fwd,
4317  Array<OneD,NekDouble> &Bwd,
4318  bool PutFwdInBwdOnBCs)
4319  {
4320  boost::ignore_unused(Fwd, Bwd, PutFwdInBwdOnBCs);
4321  NEKERROR(ErrorUtil::efatal,
4322  "This method is not defined or valid for this class type");
4323  }
4324 
4325  void ExpList::v_AddTraceQuadPhysToField(
4326  const Array<OneD, const NekDouble> &Fwd,
4327  const Array<OneD, const NekDouble> &Bwd,
4328  Array<OneD, NekDouble> &field)
4329  {
4330  boost::ignore_unused(field, Fwd, Bwd);
4331  NEKERROR(ErrorUtil::efatal,
4332  "v_AddTraceQuadPhysToField is not defined for this class type");
4333  }
4334 
4335  void ExpList::v_AddTraceQuadPhysToOffDiag(
4336  const Array<OneD, const NekDouble> &Fwd,
4337  const Array<OneD, const NekDouble> &Bwd,
4338  Array<OneD, NekDouble> &field)
4339  {
4340  boost::ignore_unused(field, Fwd, Bwd);
4341  NEKERROR(ErrorUtil::efatal,
4342  "v_AddTraceQuadPhysToOffDiag is not defined for this class");
4343  }
4344 
4345  void ExpList::v_GetLocTraceFromTracePts(
4346  const Array<OneD, const NekDouble> &Fwd,
4347  const Array<OneD, const NekDouble> &Bwd,
4348  Array<OneD, NekDouble> &locTraceFwd,
4349  Array<OneD, NekDouble> &locTraceBwd)
4350  {
4351  boost::ignore_unused(Fwd,Bwd,locTraceFwd,locTraceBwd);
4352  NEKERROR(ErrorUtil::efatal,
4353  "v_GetLocTraceFromTracePts is not defined for this class");
4354  }
4355 
4356  const Array<OneD,const NekDouble>
4357  &ExpList::v_GetBndCondBwdWeight()
4358  {
4359  NEKERROR(ErrorUtil::efatal,
4360  "v_GetBndCondBwdWeight is not defined for this class type");
4361  static Array<OneD, NekDouble> tmp;
4362  return tmp;
4363  }
4364 
4365  void ExpList::v_SetBndCondBwdWeight(
4366  const int index,
4367  const NekDouble value)
4368  {
4369  boost::ignore_unused(index, value);
4370  NEKERROR(ErrorUtil::efatal,
4371  "v_setBndCondBwdWeight is not defined for this class type");
4372  }
4373 
4374  const vector<bool> &ExpList::v_GetLeftAdjacentFaces(void) const
4375  {
4376  NEKERROR(ErrorUtil::efatal,
4377  "This method is not defined or valid for this class type");
4378  static vector<bool> tmp;
4379  return tmp;
4380  }
4381 
4382 
4383  void ExpList::v_ExtractTracePhys(Array<OneD,NekDouble> &outarray)
4384  {
4385  boost::ignore_unused(outarray);
4386  NEKERROR(ErrorUtil::efatal,
4387  "This method is not defined or valid for this class type");
4388  }
4389 
4390  void ExpList::v_ExtractTracePhys(
4391  const Array<OneD, const NekDouble> &inarray,
4392  Array<OneD,NekDouble> &outarray)
4393  {
4394  boost::ignore_unused(inarray, outarray);
4395  NEKERROR(ErrorUtil::efatal,
4396  "This method is not defined or valid for this class type");
4397  }
4398 
4399  void ExpList::v_MultiplyByInvMassMatrix(
4400  const Array<OneD,const NekDouble> &inarray,
4401  Array<OneD, NekDouble> &outarray)
4402  {
4403  boost::ignore_unused(inarray, outarray);
4404  NEKERROR(ErrorUtil::efatal,
4405  "This method is not defined or valid for this class type");
4406  }
4407 
4408  void ExpList::v_HelmSolve(
4409  const Array<OneD, const NekDouble> &inarray,
4410  Array<OneD, NekDouble> &outarray,
4411  const StdRegions::ConstFactorMap &factors,
4412  const StdRegions::VarCoeffMap &varcoeff,
4413  const MultiRegions::VarFactorsMap &varfactors,
4414  const Array<OneD, const NekDouble> &dirForcing,
4415  const bool PhysSpaceForcing)
4416  {
4417  boost::ignore_unused(inarray, outarray, factors, varcoeff,
4418  varfactors, dirForcing, PhysSpaceForcing);
4419  NEKERROR(ErrorUtil::efatal, "HelmSolve not implemented.");
4420  }
4421 
4422  void ExpList::v_LinearAdvectionDiffusionReactionSolve(
4423  const Array<OneD, Array<OneD, NekDouble> > &velocity,
4424  const Array<OneD, const NekDouble> &inarray,
4425  Array<OneD, NekDouble> &outarray,
4426  const NekDouble lambda,
4427  const Array<OneD, const NekDouble>& dirForcing)
4428  {
4429  boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4430  NEKERROR(ErrorUtil::efatal,
4431  "This method is not defined or valid for this class type");
4432  }
4433 
4434  void ExpList::v_LinearAdvectionReactionSolve(
4435  const Array<OneD, Array<OneD, NekDouble> > &velocity,
4436  const Array<OneD, const NekDouble> &inarray,
4437  Array<OneD, NekDouble> &outarray,
4438  const NekDouble lambda,
4439  const Array<OneD, const NekDouble>& dirForcing)
4440  {
4441  boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4442  NEKERROR(ErrorUtil::efatal,
4443  "This method is not defined or valid for this class type");
4444  }
4445 
4446  void ExpList::v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
4447  Array<OneD, NekDouble> &outarray,
4448  bool Shuff,
4449  bool UnShuff)
4450  {
4451  boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4452  NEKERROR(ErrorUtil::efatal,
4453  "This method is not defined or valid for this class type");
4454  }
4455 
4456  void ExpList::v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
4457  Array<OneD, NekDouble> &outarray,
4458  bool Shuff,
4459  bool UnShuff)
4460  {
4461  boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4462  NEKERROR(ErrorUtil::efatal,
4463  "This method is not defined or valid for this class type");
4464  }
4465 
4466  void ExpList::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,
4467  const Array<OneD, NekDouble> &inarray2,
4468  Array<OneD, NekDouble> &outarray)
4469  {
4470  boost::ignore_unused(inarray1, inarray2, outarray);
4471  NEKERROR(ErrorUtil::efatal,
4472  "This method is not defined or valid for this class type");
4473  }
4474 
4475  void ExpList::v_DealiasedDotProd(
4476  const Array<OneD, Array<OneD, NekDouble> > &inarray1,
4477  const Array<OneD, Array<OneD, NekDouble> > &inarray2,
4478  Array<OneD, Array<OneD, NekDouble> > &outarray)
4479  {
4480  boost::ignore_unused(inarray1, inarray2, outarray);
4481  NEKERROR(ErrorUtil::efatal,
4482  "This method is not defined or valid for this class type");
4483  }
4484 
4485  void ExpList::v_GetBCValues(Array<OneD, NekDouble> &BndVals,
4486  const Array<OneD, NekDouble> &TotField,
4487  int BndID)
4488  {
4489  boost::ignore_unused(BndVals, TotField, BndID);
4490  NEKERROR(ErrorUtil::efatal,
4491  "This method is not defined or valid for this class type");
4492  }
4493 
4494  void ExpList::v_NormVectorIProductWRTBase(Array<OneD, const NekDouble> &V1,
4495  Array<OneD, const NekDouble> &V2,
4496  Array<OneD, NekDouble> &outarray,
4497  int BndID)
4498  {
4499  boost::ignore_unused(V1, V2, outarray, BndID);
4500  NEKERROR(ErrorUtil::efatal,
4501  "This method is not defined or valid for this class type");
4502  }
4503 
4504  void ExpList::v_NormVectorIProductWRTBase(
4505  Array<OneD, Array<OneD, NekDouble> > &V,
4506  Array<OneD, NekDouble> &outarray)
4507  {
4508  Array< OneD, NekDouble> tmp;
4509  switch (GetCoordim(0))
4510  {
4511  case 1:
4512  {
4513  for(int i = 0; i < GetExpSize(); ++i)
4514  {
4515  (*m_exp)[i]->NormVectorIProductWRTBase(
4516  V[0] + GetPhys_Offset(i),
4517  tmp = outarray + GetCoeff_Offset(i));
4518  }
4519  }
4520  break;
4521  case 2:
4522  {
4523  for(int i = 0; i < GetExpSize(); ++i)
4524  {
4525  (*m_exp)[i]->NormVectorIProductWRTBase(
4526  V[0] + GetPhys_Offset(i),
4527  V[1] + GetPhys_Offset(i),
4528  tmp = outarray + GetCoeff_Offset(i));
4529  }
4530  }
4531  break;
4532  case 3:
4533  {
4534  for(int i = 0; i < GetExpSize(); ++i)
4535  {
4536  (*m_exp)[i]->NormVectorIProductWRTBase(
4537  V[0] + GetPhys_Offset(i),
4538  V[1] + GetPhys_Offset(i),
4539  V[2] + GetPhys_Offset(i),
4540  tmp = outarray + GetCoeff_Offset(i));
4541  }
4542  }
4543  break;
4544  default:
4545  NEKERROR(ErrorUtil::efatal,"Dimension not supported");
4546  break;
4547  }
4548  }
4549 
4550  void ExpList::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
4551  {
4552  boost::ignore_unused(outarray);
4553  NEKERROR(ErrorUtil::efatal,
4554  "This method is not defined or valid for this class type");
4555  }
4556 
4557  /**
4558  */
4559  void ExpList::v_FillBndCondFromField()
4560  {
4561  NEKERROR(ErrorUtil::efatal,
4562  "This method is not defined or valid for this class type");
4563  }
4564 
4565  /**
4566  */
4567  void ExpList::v_FillBndCondFromField(const int nreg)
4568  {
4569  boost::ignore_unused(nreg);
4570  NEKERROR(ErrorUtil::efatal,
4571  "This method is not defined or valid for this class type");
4572  }
4573 
4574  void ExpList::v_LocalToGlobal(bool useComm)
4575  {
4576  boost::ignore_unused(useComm);
4577  NEKERROR(ErrorUtil::efatal,
4578  "This method is not defined or valid for this class type");
4579  }
4580 
4581 
4582  void ExpList::v_LocalToGlobal(const Array<OneD, const NekDouble> &inarray,
4583  Array<OneD,NekDouble> &outarray,
4584  bool useComm)
4585  {
4586  boost::ignore_unused(inarray, outarray, useComm);
4587  NEKERROR(ErrorUtil::efatal,
4588  "This method is not defined or valid for this class type");
4589  }
4590 
4591 
4592  void ExpList::v_GlobalToLocal(void)
4593  {
4594  NEKERROR(ErrorUtil::efatal,
4595  "This method is not defined or valid for this class type");
4596  }
4597 
4598 
4599  void ExpList::v_GlobalToLocal(const Array<OneD, const NekDouble> &inarray,
4600  Array<OneD,NekDouble> &outarray)
4601  {
4602  boost::ignore_unused(inarray, outarray);
4603  NEKERROR(ErrorUtil::efatal,
4604  "This method is not defined or valid for this class type");
4605  }
4606 
4607 
4608  void ExpList::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
4609  Array<OneD, NekDouble> &outarray)
4610  {
4611  v_BwdTrans_IterPerExp(inarray,outarray);
4612  }
4613 
4614  void ExpList::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
4615  Array<OneD, NekDouble> &outarray)
4616  {
4617  v_FwdTrans_IterPerExp(inarray,outarray);
4618  }
4619 
4620  void ExpList::v_IProductWRTBase(
4621  const Array<OneD, const NekDouble> &inarray,
4622  Array<OneD, NekDouble> &outarray)
4623  {
4624  // initialise if required
4625  if(m_collectionsDoInit[Collections::eIProductWRTBase])
4626  {
4627  for (int i = 0; i < m_collections.size(); ++i)
4628  {
4629  m_collections[i].Initialise(Collections::eIProductWRTBase);
4630  }
4631  m_collectionsDoInit[Collections::eIProductWRTBase] = false;
4632  }
4633 
4634  Array<OneD,NekDouble> tmp;
4635  for (int i = 0; i < m_collections.size(); ++i)
4636  {
4637 
4638  m_collections[i].ApplyOperator(Collections::eIProductWRTBase,
4639  inarray + m_coll_phys_offset[i],
4640  tmp = outarray + m_coll_coeff_offset[i]);
4641  }
4642  }
4643 
4644  void ExpList::v_GeneralMatrixOp(
4645  const GlobalMatrixKey &gkey,
4646  const Array<OneD,const NekDouble> &inarray,
4647  Array<OneD, NekDouble> &outarray)
4648  {
4649  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
4650  }
4651 
4652  /**
4653  * The operation is evaluated locally by the elemental
4654  * function StdRegions#StdExpansion#GetCoords.
4655  *
4656  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
4657  * will be stored in this array.
4658  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
4659  * will be stored in this array.
4660  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
4661  * will be stored in this array.
4662  */
4663  void ExpList::v_GetCoords(Array<OneD, NekDouble> &coord_0,
4664  Array<OneD, NekDouble> &coord_1,
4665  Array<OneD, NekDouble> &coord_2)
4666  {
4667  if (GetNumElmts() == 0)
4668  {
4669  return;
4670  }
4671 
4672  int i;
4673  Array<OneD, NekDouble> e_coord_0;
4674  Array<OneD, NekDouble> e_coord_1;
4675  Array<OneD, NekDouble> e_coord_2;
4676 
4677  switch(GetExp(0)->GetCoordim())
4678  {
4679  case 1:
4680  for(i= 0; i < (*m_exp).size(); ++i)
4681  {
4682  e_coord_0 = coord_0 + m_phys_offset[i];
4683  (*m_exp)[i]->GetCoords(e_coord_0);
4684  }
4685  break;
4686  case 2:
4687  ASSERTL0(coord_1.size() != 0,
4688  "output coord_1 is not defined");
4689 
4690  for(i= 0; i < (*m_exp).size(); ++i)
4691  {
4692  e_coord_0 = coord_0 + m_phys_offset[i];
4693  e_coord_1 = coord_1 + m_phys_offset[i];
4694  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1);
4695  }
4696  break;
4697  case 3:
4698  ASSERTL0(coord_1.size() != 0,
4699  "output coord_1 is not defined");
4700  ASSERTL0(coord_2.size() != 0,
4701  "output coord_2 is not defined");
4702 
4703  for(i= 0; i < (*m_exp).size(); ++i)
4704  {
4705  e_coord_0 = coord_0 + m_phys_offset[i];
4706  e_coord_1 = coord_1 + m_phys_offset[i];
4707  e_coord_2 = coord_2 + m_phys_offset[i];
4708  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1,e_coord_2);
4709  }
4710  break;
4711  }
4712  }
4713 
4714  /**
4715  */
4716  void ExpList::v_SetUpPhysNormals()
4717  {
4718  for (int i = 0; i < m_exp->size(); ++i)
4719  {
4720  for (int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
4721  {
4722  (*m_exp)[i]->ComputeTraceNormal(j);
4723  }
4724  }
4725  }
4726 
4727  /**
4728  */
4729  void ExpList::v_GetBndElmtExpansion(int i,
4730  std::shared_ptr<ExpList> &result,
4731  const bool DeclareCoeffPhysArrays)
4732  {
4733  boost::ignore_unused(i, result, DeclareCoeffPhysArrays);
4734  NEKERROR(ErrorUtil::efatal,
4735  "This method is not defined or valid for this class type");
4736  }
4737 
4738  /**
4739  */
4740  void ExpList::v_ExtractElmtToBndPhys(const int i,
4741  const Array<OneD, NekDouble> &element,
4742  Array<OneD, NekDouble> &boundary)
4743  {
4744  int n, cnt;
4745  Array<OneD, NekDouble> tmp1, tmp2;
4746  LocalRegions::ExpansionSharedPtr elmt;
4747 
4748  Array<OneD, int> ElmtID,EdgeID;
4749  GetBoundaryToElmtMap(ElmtID,EdgeID);
4750 
4751  // Initialise result
4752  boundary = Array<OneD, NekDouble>
4753  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
4754 
4755  // Skip other boundary regions
4756  for (cnt = n = 0; n < i; ++n)
4757  {
4758  cnt += GetBndCondExpansions()[n]->GetExpSize();
4759  }
4760 
4761  int offsetBnd;
4762  int offsetElmt = 0;
4763  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4764  {
4765  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4766 
4767  elmt = GetExp(ElmtID[cnt+n]);
4768  elmt->GetTracePhysVals(EdgeID[cnt+n],
4769  GetBndCondExpansions()[i]->GetExp(n),
4770  tmp1 = element + offsetElmt,
4771  tmp2 = boundary + offsetBnd);
4772 
4773  offsetElmt += elmt->GetTotPoints();
4774  }
4775  }
4776 
4777  /**
4778  */
4779  void ExpList::v_ExtractPhysToBndElmt(const int i,
4780  const Array<OneD, const NekDouble> &phys,
4781  Array<OneD, NekDouble> &bndElmt)
4782  {
4783  int n, cnt, nq;
4784 
4785  Array<OneD, int> ElmtID,EdgeID;
4786  GetBoundaryToElmtMap(ElmtID,EdgeID);
4787 
4788  // Skip other boundary regions
4789  for (cnt = n = 0; n < i; ++n)
4790  {
4791  cnt += GetBndCondExpansions()[n]->GetExpSize();
4792  }
4793 
4794  // Count number of points
4795  int npoints = 0;
4796  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4797  {
4798  npoints += GetExp(ElmtID[cnt+n])->GetTotPoints();
4799  }
4800 
4801  // Initialise result
4802  bndElmt = Array<OneD, NekDouble> (npoints, 0.0);
4803 
4804  // Extract data
4805  int offsetPhys;
4806  int offsetElmt = 0;
4807  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4808  {
4809  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
4810  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
4811  Vmath::Vcopy(nq, &phys[offsetPhys], 1,
4812  &bndElmt[offsetElmt], 1);
4813  offsetElmt += nq;
4814  }
4815  }
4816 
4817  /**
4818  */
4819  void ExpList::v_ExtractPhysToBnd(int i,
4820  const Array<OneD, const NekDouble> &phys,
4821  Array<OneD, NekDouble> &bnd)
4822  {
4823  int n, cnt;
4824  Array<OneD, NekDouble> tmp1;
4825  LocalRegions::ExpansionSharedPtr elmt;
4826 
4827  Array<OneD, int> ElmtID,EdgeID;
4828  GetBoundaryToElmtMap(ElmtID,EdgeID);
4829 
4830  // Initialise result
4831  bnd = Array<OneD, NekDouble>
4832  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
4833 
4834  // Skip other boundary regions
4835  for (cnt = n = 0; n < i; ++n)
4836  {
4837  cnt += GetBndCondExpansions()[n]->GetExpSize();
4838  }
4839 
4840  int offsetBnd;
4841  int offsetPhys;
4842  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4843  {
4844  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
4845  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4846 
4847  elmt = GetExp(ElmtID[cnt+n]);
4848  elmt->GetTracePhysVals(EdgeID[cnt+n],
4849  GetBndCondExpansions()[i]->GetExp(n),
4850  phys + offsetPhys,
4851  tmp1 = bnd + offsetBnd);
4852  }
4853  }
4854 
4855  /**
4856  */
4857  void ExpList::v_GetBoundaryNormals(int i,
4858  Array<OneD, Array<OneD, NekDouble> > &normals)
4859  {
4860  int j, n, cnt, nq;
4861  int coordim = GetCoordim(0);
4862  Array<OneD, NekDouble> tmp;
4863  LocalRegions::ExpansionSharedPtr elmt;
4864 
4865  Array<OneD, int> ElmtID,EdgeID;
4866  GetBoundaryToElmtMap(ElmtID,EdgeID);
4867 
4868  // Initialise result
4869  normals = Array<OneD, Array<OneD, NekDouble> > (coordim);
4870  for (j = 0; j < coordim; ++j)
4871  {
4872  normals[j] = Array<OneD, NekDouble> (
4873  GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
4874  }
4875 
4876  // Skip other boundary regions
4877  for (cnt = n = 0; n < i; ++n)
4878  {
4879  cnt += GetBndCondExpansions()[n]->GetExpSize();
4880  }
4881 
4882  int offset;
4883  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
4884  {
4885  offset = GetBndCondExpansions()[i]->GetPhys_Offset(n);
4886  nq = GetBndCondExpansions()[i]->GetExp(n)->GetTotPoints();
4887 
4888  elmt = GetExp(ElmtID[cnt+n]);
4889  const Array<OneD, const Array<OneD, NekDouble> > normalsElmt
4890  = elmt->GetTraceNormal(EdgeID[cnt+n]);
4891  // Copy to result
4892  for (j = 0; j < coordim; ++j)
4893  {
4894  Vmath::Vcopy(nq, normalsElmt[j], 1,
4895  tmp = normals[j] + offset, 1);
4896  }
4897  }
4898  }
4899 
4900  /**
4901  */
4902  void ExpList::v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
4903  Array<OneD,int> &EdgeID)
4904  {
4905  boost::ignore_unused(ElmtID, EdgeID);
4906  NEKERROR(ErrorUtil::efatal,
4907  "This method is not defined or valid for this class type");
4908  }
4909 
4910  void ExpList::v_FillBwdWithBwdWeight(
4911  Array<OneD, NekDouble> &weightave,
4912  Array<OneD, NekDouble> &weightjmp)
4913  {
4914  boost::ignore_unused(weightave, weightjmp);
4915  NEKERROR(ErrorUtil::efatal, "v_FillBwdWithBwdWeight not defined");
4916  }
4917 
4918  void ExpList::v_PeriodicBwdCopy(
4919  const Array<OneD, const NekDouble> &Fwd,
4920  Array<OneD, NekDouble> &Bwd)
4921  {
4922  boost::ignore_unused(Fwd, Bwd);
4923  NEKERROR(ErrorUtil::efatal, "v_PeriodicBwdCopy not defined");
4924  }
4925 
4926  /**
4927  */
4928  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>
4929  &ExpList::v_GetBndConditions(void)
4930  {
4931  NEKERROR(ErrorUtil::efatal,
4932  "This method is not defined or valid for this class type");
4933  static Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
4934  result;
4935  return result;
4936  }
4937 
4938  /**
4939  */
4940  Array<OneD,SpatialDomains::BoundaryConditionShPtr> &ExpList::v_UpdateBndConditions()
4941  {
4942  NEKERROR(ErrorUtil::efatal,
4943  "This method is not defined or valid for this class type");
4944  static Array<OneD, SpatialDomains::BoundaryConditionShPtr> result;
4945  return result;
4946  }
4947 
4948  /**
4949  */
4950  void ExpList::v_EvaluateBoundaryConditions(
4951  const NekDouble time,
4952  const std::string varName,
4953  const NekDouble x2_in,
4954  const NekDouble x3_in)
4955  {
4956  boost::ignore_unused(time, varName, x2_in, x3_in);
4957  NEKERROR(ErrorUtil::efatal,
4958  "This method is not defined or valid for this class type");
4959  }
4960 
4961  /**
4962  */
4963  map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
4964  {
4965  NEKERROR(ErrorUtil::efatal,
4966  "This method is not defined or valid for this class type");
4967  static map<int,RobinBCInfoSharedPtr> result;
4968  return result;
4969  }
4970 
4971  /**
4972  */
4973  void ExpList::v_GetPeriodicEntities(
4974  PeriodicMap &periodicVerts,
4975  PeriodicMap &periodicEdges,
4976  PeriodicMap &periodicFaces)
4977  {
4978  boost::ignore_unused(periodicVerts, periodicEdges, periodicFaces);
4979  NEKERROR(ErrorUtil::efatal,
4980  "This method is not defined or valid for this class type");
4981  }
4982 
4983  SpatialDomains::BoundaryConditionShPtr ExpList::GetBoundaryCondition(
4984  const SpatialDomains::BoundaryConditionCollection& collection,
4985  unsigned int regionId,
4986  const std::string& variable)
4987  {
4988  auto collectionIter = collection.find(regionId);
4989  ASSERTL1(collectionIter != collection.end(),
4990  "Unable to locate collection " +
4991  boost::lexical_cast<string>(regionId));
4992 
4993  const SpatialDomains::BoundaryConditionMapShPtr bndCondMap
4994  = (*collectionIter).second;
4995  auto conditionMapIter = bndCondMap->find(variable);
4996  ASSERTL1(conditionMapIter != bndCondMap->end(),
4997  "Unable to locate condition map.");
4998 
4999  const SpatialDomains::BoundaryConditionShPtr boundaryCondition
5000  = (*conditionMapIter).second;
5001 
5002  return boundaryCondition;
5003  }
5004 
5005  ExpListSharedPtr &ExpList::v_GetPlane(int n)
5006  {
5007  boost::ignore_unused(n);
5008  NEKERROR(ErrorUtil::efatal,
5009  "This method is not defined or valid for this class type");
5010  return NullExpListSharedPtr;
5011  }
5012 
5013  /**
5014  * @brief Construct collections of elements containing a single element
5015  * type and polynomial order from the list of expansions.
5016  */
5017  void ExpList::CreateCollections(Collections::ImplementationType ImpType)
5018  {
5019  map<LibUtilities::ShapeType,
5020  vector<std::pair<LocalRegions::ExpansionSharedPtr,int> > >
5021  collections;
5022 
5023  //Set up initialisation flags
5024  m_collectionsDoInit = std::vector<bool>(Collections::SIZE_OperatorType,true);
5025 
5026  // Figure out optimisation parameters if provided in
5027  // session file or default given
5028  Collections::CollectionOptimisation colOpt(m_session, ImpType);
5029  ImpType = colOpt.GetDefaultImplementationType();
5030 
5031  bool autotuning = colOpt.IsUsingAutotuning();
5032  bool verbose = (m_session->DefinesCmdLineArgument("verbose")) &&
5033  (m_comm->GetRank() == 0);
5034  int collmax = (colOpt.GetMaxCollectionSize() > 0
5035  ? colOpt.GetMaxCollectionSize()
5036  : 2*m_exp->size());
5037 
5038  // clear vectors in case previously called
5039  m_collections.clear();
5040  m_coll_coeff_offset.clear();
5041  m_coll_phys_offset.clear();
5042 
5043  // Loop over expansions, and create collections for each element type
5044  for (int i = 0; i < m_exp->size(); ++i)
5045  {
5046  collections[(*m_exp)[i]->DetShapeType()].push_back(
5047  std::pair<LocalRegions::ExpansionSharedPtr,int>((*m_exp)[i],i));
5048  }
5049 
5050  for (auto &it : collections)
5051  {
5052  LocalRegions::ExpansionSharedPtr exp = it.second[0].first;
5053 
5054  Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(exp);
5055  vector<StdRegions::StdExpansionSharedPtr> collExp;
5056 
5057  int prevCoeffOffset = m_coeff_offset[it.second[0].second];
5058  int prevPhysOffset = m_phys_offset [it.second[0].second];
5059  int collcnt;
5060 
5061  m_coll_coeff_offset.push_back(prevCoeffOffset);
5062  m_coll_phys_offset .push_back(prevPhysOffset);
5063 
5064  if(it.second.size() == 1) // single element case
5065  {
5066  collExp.push_back(it.second[0].first);
5067 
5068  // if no Imp Type provided and No
5069  // settign in xml file. reset
5070  // impTypes using timings
5071  if(autotuning)
5072  {
5073  impTypes = colOpt.SetWithTimings(collExp,
5074  impTypes, verbose);
5075  }
5076 
5077  Collections::Collection tmp(collExp, impTypes);
5078  m_collections.push_back(tmp);
5079  }
5080  else
5081  {
5082  // set up first geometry
5083  collExp.push_back(it.second[0].first);
5084  int prevnCoeff = it.second[0].first->GetNcoeffs();
5085  int prevnPhys = it.second[0].first->GetTotPoints();
5086  bool prevDeformed = it.second[0].first->GetMetricInfo()->GetGtype()
5087  == SpatialDomains::eDeformed;
5088  collcnt = 1;
5089 
5090  for (int i = 1; i < it.second.size(); ++i)
5091  {
5092  int nCoeffs = it.second[i].first->GetNcoeffs();
5093  int nPhys = it.second[i].first->GetTotPoints();
5094  bool Deformed = it.second[i].first->GetMetricInfo()->GetGtype()
5095  == SpatialDomains::eDeformed;
5096  int coeffOffset = m_coeff_offset[it.second[i].second];
5097  int physOffset = m_phys_offset [it.second[i].second];
5098 
5099  // check to see if next elmt is different or
5100  // collmax reached and if so end collection
5101  // and start new one
5102  if(prevCoeffOffset + nCoeffs != coeffOffset ||
5103  prevnCoeff != nCoeffs ||
5104  prevPhysOffset + nPhys != physOffset ||
5105  prevDeformed != Deformed ||
5106  prevnPhys != nPhys || collcnt >= collmax)
5107  {
5108 
5109  // if no Imp Type provided and No
5110  // settign in xml file. reset
5111  // impTypes using timings
5112  if(autotuning)
5113  {
5114  impTypes = colOpt.SetWithTimings(collExp,
5115  impTypes,
5116  verbose);
5117  }
5118 
5119  Collections::Collection tmp(collExp, impTypes);
5120  m_collections.push_back(tmp);
5121 
5122  // start new geom list
5123  collExp.clear();
5124 
5125  m_coll_coeff_offset.push_back(coeffOffset);
5126  m_coll_phys_offset .push_back(physOffset);
5127  collExp.push_back(it.second[i].first);
5128  collcnt = 1;
5129  }
5130  else // add to list of collections
5131  {
5132  collExp.push_back(it.second[i].first);
5133  collcnt++;
5134  }
5135 
5136  // if end of list finish up collection
5137  if (i == it.second.size() - 1)
5138  {
5139  // if no Imp Type provided and No
5140  // settign in xml file.
5141  if(autotuning)
5142  {
5143  impTypes = colOpt.SetWithTimings(collExp,
5144  impTypes,verbose);
5145  }
5146 
5147  Collections::Collection tmp(collExp, impTypes);
5148  m_collections.push_back(tmp);
5149 
5150  collExp.clear();
5151  collcnt = 0;
5152 
5153  }
5154 
5155  prevCoeffOffset = coeffOffset;
5156  prevPhysOffset = physOffset;
5157  prevDeformed = Deformed;
5158  prevnCoeff = nCoeffs;
5159  prevnPhys = nPhys;
5160  }
5161  }
5162  }
5163  }
5164 
5165  void ExpList::ClearGlobalLinSysManager(void)
5166  {
5167  v_ClearGlobalLinSysManager();
5168  }
5169 
5170  void ExpList::v_PhysInterp1DScaled(
5171  const NekDouble scale,
5172  const Array<OneD, NekDouble> &inarray,
5173  Array<OneD, NekDouble> &outarray)
5174  {
5175  int cnt,cnt1;
5176 
5177  cnt = cnt1 = 0;
5178 
5179  switch(m_expType)
5180  {
5181  case e2D:
5182  {
5183  for(int i = 0; i < GetExpSize(); ++i)
5184  {
5185  // get new points key
5186  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5187  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5188  int npt0 = (int) pt0*scale;
5189  int npt1 = (int) pt1*scale;
5190 
5191  LibUtilities::PointsKey newPointsKey0(npt0,
5192  (*m_exp)[i]->GetPointsType(0));
5193  LibUtilities::PointsKey newPointsKey1(npt1,
5194  (*m_exp)[i]->GetPointsType(1));
5195 
5196  // Interpolate points;
5197  LibUtilities::Interp2D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5198  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5199  &inarray[cnt],newPointsKey0,
5200  newPointsKey1,&outarray[cnt1]);
5201 
5202  cnt += pt0*pt1;
5203  cnt1 += npt0*npt1;
5204  }
5205  }
5206  break;
5207  case e3D:
5208  {
5209  for(int i = 0; i < GetExpSize(); ++i)
5210  {
5211  // get new points key
5212  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5213  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5214  int pt2 = (*m_exp)[i]->GetNumPoints(2);
5215  int npt0 = (int) pt0*scale;
5216  int npt1 = (int) pt1*scale;
5217  int npt2 = (int) pt2*scale;
5218 
5219  LibUtilities::PointsKey newPointsKey0(npt0,(*m_exp)[i]->GetPointsType(0));
5220  LibUtilities::PointsKey newPointsKey1(npt1,(*m_exp)[i]->GetPointsType(1));
5221  LibUtilities::PointsKey newPointsKey2(npt2,(*m_exp)[i]->GetPointsType(2));
5222 
5223  // Interpolate points;
5224  LibUtilities::Interp3D((*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5225  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5226  (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
5227  &inarray[cnt], newPointsKey0,
5228  newPointsKey1, newPointsKey2,
5229  &outarray[cnt1]);
5230 
5231  cnt += pt0*pt1*pt2;
5232  cnt1 += npt0*npt1*npt2;
5233  }
5234  }
5235  break;
5236  default:
5237  {
5238  NEKERROR(ErrorUtil::efatal,"This expansion is not set");
5239  }
5240  break;
5241  }
5242  }
5243 
5244  void ExpList::v_AddTraceIntegralToOffDiag(
5245  const Array<OneD, const NekDouble> &FwdFlux,
5246  const Array<OneD, const NekDouble> &BwdFlux,
5247  Array<OneD, NekDouble> &outarray)
5248  {
5249  boost::ignore_unused(FwdFlux, BwdFlux, outarray);
5250  NEKERROR(ErrorUtil::efatal,"AddTraceIntegralToOffDiag not defined");
5251  }
5252 
5253  void ExpList::GetMatIpwrtDeriveBase(
5254  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
5255  const int nDirctn,
5256  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5257  {
5258  int nTotElmt = (*m_exp).size();
5259  int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5260  int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5261 
5262 
5263  Array<OneD, NekDouble> tmpCoef(nElmtCoef,0.0);
5264  Array<OneD, NekDouble> tmpPhys(nElmtPnt,0.0);
5265 
5266  for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
5267  {
5268  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5269  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5270 
5271  if (tmpPhys.size()!=nElmtPnt||tmpCoef.size()!=nElmtCoef)
5272  {
5273  tmpPhys = Array<OneD, NekDouble>(nElmtPnt,0.0);
5274  tmpCoef = Array<OneD, NekDouble>(nElmtCoef,0.0);
5275  }
5276 
5277  for(int ncl = 0; ncl < nElmtPnt; ncl++)
5278  {
5279  tmpPhys[ncl] = inarray[nelmt][ncl];
5280 
5281  (*m_exp)[nelmt]->IProductWRTDerivBase(nDirctn,tmpPhys,
5282  tmpCoef);
5283 
5284  for(int nrw = 0; nrw < nElmtCoef; nrw++)
5285  {
5286  (*mtxPerVar[nelmt])(nrw,ncl) = tmpCoef[nrw];
5287  }
5288  // to maintain all the other columes are zero.
5289  tmpPhys[ncl] = 0.0;
5290  }
5291  }
5292  }
5293 
5294  void ExpList::GetMatIpwrtDeriveBase(
5295  const TensorOfArray3D<NekDouble> &inarray,
5296  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5297  {
5298  int nTotElmt = (*m_exp).size();
5299 
5300  int nspacedim = m_graph->GetSpaceDimension();
5301  Array<OneD, Array<OneD, NekDouble> > projectedpnts(nspacedim);
5302  Array<OneD, Array<OneD, NekDouble> > tmppnts(nspacedim);
5303  Array<OneD, DNekMatSharedPtr> ArrayStdMat(nspacedim);
5304  Array<OneD, Array<OneD, NekDouble> > ArrayStdMat_data(nspacedim);
5305 
5306  Array<OneD, NekDouble > clmnArray,clmnStdMatArray;
5307 
5308  LibUtilities::ShapeType ElmtTypePrevious =
5309  LibUtilities::eNoShapeType;
5310  int nElmtPntPrevious = 0;
5311  int nElmtCoefPrevious = 0;
5312 
5313  int nElmtPnt,nElmtCoef;
5314  for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
5315  {
5316  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5317  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5318  LibUtilities::ShapeType ElmtTypeNow =
5319  (*m_exp)[nelmt]->DetShapeType();
5320 
5321  if (nElmtPntPrevious!=nElmtPnt||nElmtCoefPrevious!=nElmtCoef||
5322  (ElmtTypeNow!=ElmtTypePrevious))
5323  {
5324  if(nElmtPntPrevious!=nElmtPnt)
5325  {
5326  for(int ndir =0;ndir<nspacedim;ndir++)
5327  {
5328  projectedpnts[ndir] =
5329  Array<OneD, NekDouble>(nElmtPnt,0.0);
5330  tmppnts[ndir] =
5331  Array<OneD, NekDouble>(nElmtPnt,0.0);
5332  }
5333  }
5334  StdRegions::StdExpansionSharedPtr stdExp;
5335  stdExp = (*m_exp)[nelmt]->GetStdExp();
5336  StdRegions::StdMatrixKey matkey(StdRegions::eDerivBase0,
5337  stdExp->DetShapeType(), *stdExp);
5338 
5339  ArrayStdMat[0] = stdExp->GetStdMatrix(matkey);
5340  ArrayStdMat_data[0] = ArrayStdMat[0]->GetPtr();
5341 
5342  if(nspacedim>1)
5343  {
5344  StdRegions::StdMatrixKey matkey(StdRegions::eDerivBase1,
5345  stdExp->DetShapeType(), *stdExp);
5346 
5347  ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
5348  ArrayStdMat_data[1] = ArrayStdMat[1]->GetPtr();
5349 
5350  if(nspacedim>2)
5351  {
5352  StdRegions::StdMatrixKey matkey(
5353  StdRegions::eDerivBase2,
5354  stdExp->DetShapeType(), *stdExp);
5355 
5356  ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
5357  ArrayStdMat_data[2] = ArrayStdMat[2]->GetPtr();
5358  }
5359  }
5360 
5361  ElmtTypePrevious = ElmtTypeNow;
5362  nElmtPntPrevious = nElmtPnt;
5363  nElmtCoefPrevious = nElmtCoef;
5364  }
5365  else
5366  {
5367  for(int ndir =0;ndir<nspacedim;ndir++)
5368  {
5369  Vmath::Zero(nElmtPnt,projectedpnts[ndir],1);
5370  }
5371  }
5372 
5373  for(int ndir =0;ndir<nspacedim;ndir++)
5374  {
5375  (*m_exp)[nelmt]->AlignVectorToCollapsedDir(
5376  ndir,inarray[ndir][nelmt],tmppnts);
5377  for(int n =0;n<nspacedim;n++)
5378  {
5379  Vmath::Vadd(nElmtPnt,tmppnts[n],1,projectedpnts[n],1,
5380  projectedpnts[n],1);
5381  }
5382  }
5383 
5384  for(int ndir =0;ndir<nspacedim;ndir++)
5385  {
5386  // weight with metric
5387  (*m_exp)[nelmt]->MultiplyByQuadratureMetric(
5388  projectedpnts[ndir],projectedpnts[ndir]);
5389  Array<OneD, NekDouble> MatDataArray =
5390  mtxPerVar[nelmt]->GetPtr();
5391 
5392  for(int np=0;np<nElmtPnt;np++)
5393  {
5394  NekDouble factor = projectedpnts[ndir][np];
5395  clmnArray = MatDataArray + np*nElmtCoef;
5396  clmnStdMatArray = ArrayStdMat_data[ndir] + np*nElmtCoef;
5397  Vmath::Svtvp(nElmtCoef,factor,clmnStdMatArray,1,
5398  clmnArray,1,clmnArray,1);
5399  }
5400  }
5401  }
5402  }
5403 
5404  // TODO: Reduce cost by getting Bwd Matrix directly
5405  void ExpList::GetDiagMatIpwrtBase(
5406  const Array<OneD, const Array<OneD, NekDouble> > &inarray,
5407  Array<OneD, DNekMatSharedPtr> &mtxPerVar)
5408  {
5409  LibUtilities::ShapeType ElmtTypePrevious =
5410  LibUtilities::eNoShapeType;
5411  int nElmtPntPrevious = 0;
5412  int nElmtCoefPrevious = 0;
5413  int nTotElmt = (*m_exp).size();
5414  int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5415  int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5416 
5417  Array<OneD, NekDouble> tmpPhys;
5418  Array<OneD, NekDouble > clmnArray,clmnStdMatArray;
5419  Array<OneD, NekDouble > stdMat_data;
5420 
5421  for(int nelmt = 0; nelmt < nTotElmt; nelmt++)
5422  {
5423  nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5424  nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5425  LibUtilities::ShapeType ElmtTypeNow =
5426  (*m_exp)[nelmt]->DetShapeType();
5427 
5428  if (nElmtPntPrevious!=nElmtPnt||nElmtCoefPrevious!=nElmtCoef||
5429  (ElmtTypeNow!=ElmtTypePrevious))
5430  {
5431  StdRegions::StdExpansionSharedPtr stdExp;
5432  stdExp = (*m_exp)[nelmt]->GetStdExp();
5433  StdRegions::StdMatrixKey matkey(StdRegions::eBwdMat,
5434  stdExp->DetShapeType(), *stdExp);
5435 
5436  DNekMatSharedPtr BwdMat = stdExp->GetStdMatrix(matkey);
5437  stdMat_data = BwdMat->GetPtr();
5438 
5439  if (nElmtPntPrevious!=nElmtPnt)
5440  {
5441  tmpPhys = Array<OneD, NekDouble>(nElmtPnt,0.0);
5442  }
5443 
5444  ElmtTypePrevious = ElmtTypeNow;
5445  nElmtPntPrevious = nElmtPnt;
5446  nElmtCoefPrevious = nElmtCoef;
5447  }
5448 
5449  (*m_exp)[nelmt]->MultiplyByQuadratureMetric(inarray[nelmt],
5450  tmpPhys); // weight with metric
5451 
5452  Array<OneD, NekDouble> MatDataArray =
5453  mtxPerVar[nelmt]->GetPtr();
5454 
5455  for(int np=0;np<nElmtPnt;np++)
5456  {
5457  NekDouble factor = tmpPhys[np];
5458  clmnArray = MatDataArray + np*nElmtCoef;
5459  clmnStdMatArray = stdMat_data + np*nElmtCoef;
5460  Vmath::Smul(nElmtCoef,factor,clmnStdMatArray,1,clmnArray,1);
5461  }
5462  }
5463  }
5464 
5465  /**
5466  * @brief inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian,
5467  * with dimension NtotalTrace*TraceCoef*TracePhys.
5468  * return Elemental Jacobian matrix with dimension NtotalElement*ElementCoef*ElementPhys.
5469  *
5470  *
5471  * @param field is a NekDouble array which contains the 2D data
5472  * from which we wish to extract the backward and
5473  * forward orientated trace/edge arrays.
5474  * @param Fwd The resulting forwards space.
5475  * @param Bwd The resulting backwards space.
5476  */
5477  void ExpList::AddTraceJacToElmtJac(
5478  const Array<OneD, const DNekMatSharedPtr> &FwdMat,
5479  const Array<OneD, const DNekMatSharedPtr> &BwdMat,
5480  Array<OneD, DNekMatSharedPtr> &fieldMat)
5481  {
5482  MultiRegions::ExpListSharedPtr tracelist = GetTrace();
5483  std::shared_ptr<LocalRegions::ExpansionVector> traceExp=
5484  tracelist->GetExp();
5485  int ntotTrace = (*traceExp).size();
5486  int nTracePnt,nTraceCoef;
5487 
5488  std::shared_ptr<LocalRegions::ExpansionVector> fieldExp= GetExp();
5489  int nElmtCoef ;
5490 
5491  const MultiRegions::LocTraceToTraceMapSharedPtr locTraceToTraceMap =
5492  GetLocTraceToTraceMap();
5493  const Array<OneD, const Array<OneD, int >> LRAdjExpid =
5494  locTraceToTraceMap->GetLeftRightAdjacentExpId();
5495  const Array<OneD, const Array<OneD, bool>> LRAdjflag =
5496  locTraceToTraceMap->GetLeftRightAdjacentExpFlag();
5497 
5498  const Array<OneD, const Array<OneD, Array<OneD, int > > > elmtLRMap
5499  = locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffMap();
5500  const Array<OneD, const Array<OneD, Array<OneD, int > > > elmtLRSign
5501  = locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffSign();
5502  DNekMatSharedPtr ElmtMat;
5503  Array<OneD, NekDouble > ElmtMat_data;
5504  // int nclAdjExp;
5505  int nrwAdjExp;
5506  int MatIndex,nPnts;
5507  NekDouble sign = 1.0;
5508 
5509  int nTracePntsTtl = tracelist->GetTotPoints();
5510  int nlocTracePts = locTraceToTraceMap->GetNLocTracePts();
5511  int nlocTracePtsFwd = locTraceToTraceMap->GetNFwdLocTracePts();
5512  int nlocTracePtsBwd = nlocTracePts-nlocTracePtsFwd;
5513 
5514  Array<OneD, int > nlocTracePtsLR(2);
5515  nlocTracePtsLR[0] = nlocTracePtsFwd;
5516  nlocTracePtsLR[1] = nlocTracePtsBwd;
5517 
5518  size_t nFwdBwdNonZero = 0;
5519  Array<OneD, int> tmpIndex{2, -1};
5520  for (int i = 0; i < 2; ++i)
5521  {
5522  if (nlocTracePtsLR[i] > 0)
5523  {
5524  tmpIndex[nFwdBwdNonZero] = i;
5525  nFwdBwdNonZero++;
5526  }
5527  }
5528 
5529  Array<OneD, int> nlocTracePtsNonZeroIndex{nFwdBwdNonZero};
5530  for (int i = 0; i < nFwdBwdNonZero; ++i)
5531  {
5532  nlocTracePtsNonZeroIndex[i] = tmpIndex[i];
5533  }
5534 
5535  Array<OneD, NekDouble> TraceFwdPhy(nTracePntsTtl);
5536  Array<OneD, NekDouble> TraceBwdPhy(nTracePntsTtl);
5537  Array<OneD, Array<OneD, NekDouble> > tmplocTrace(2);
5538  for (int k = 0; k < 2; ++k)
5539  {
5540  tmplocTrace[k] = NullNekDouble1DArray;
5541  }
5542 
5543  for (int k = 0; k < nFwdBwdNonZero; ++k)
5544  {
5545  size_t i = nlocTracePtsNonZeroIndex[k];
5546  tmplocTrace[i] = Array<OneD, NekDouble> (nlocTracePtsLR[i]);
5547  }
5548 
5549  int nNumbElmt = fieldMat.size();
5550  Array<OneD, Array<OneD, NekDouble> > ElmtMatDataArray(nNumbElmt);
5551  Array<OneD, int> ElmtCoefArray(nNumbElmt);
5552  for(int i=0;i<nNumbElmt;i++)
5553  {
5554  ElmtMatDataArray[i] = fieldMat[i]->GetPtr();
5555  ElmtCoefArray[i] = GetNcoeffs(i);
5556  }
5557 
5558  int nTraceCoefMax = 0;
5559  int nTraceCoefMin = std::numeric_limits<int>::max();
5560  Array<OneD, int> TraceCoefArray(ntotTrace);
5561  Array<OneD, int> TracePntArray(ntotTrace);
5562  Array<OneD, int> TraceOffArray(ntotTrace);
5563  Array<OneD, Array<OneD, NekDouble> > FwdMatData(ntotTrace);
5564  Array<OneD, Array<OneD, NekDouble> > BwdMatData(ntotTrace);
5565  for(int nt = 0; nt < ntotTrace; nt++)
5566  {
5567  nTraceCoef = (*traceExp)[nt]->GetNcoeffs();
5568  nTracePnt = tracelist->GetTotPoints(nt);
5569  int noffset = tracelist->GetPhys_Offset(nt);
5570  TraceCoefArray[nt] = nTraceCoef;
5571  TracePntArray[nt] = nTracePnt;
5572  TraceOffArray[nt] = noffset;
5573  FwdMatData[nt] = FwdMat[nt]->GetPtr();
5574  BwdMatData[nt] = BwdMat[nt]->GetPtr();
5575  if(nTraceCoef>nTraceCoefMax)
5576  {
5577  nTraceCoefMax = nTraceCoef;
5578  }
5579  if(nTraceCoef<nTraceCoefMin)
5580  {
5581  nTraceCoefMin = nTraceCoef;
5582  }
5583  }
5584  WARNINGL1(nTraceCoefMax==nTraceCoefMin,
5585  "nTraceCoefMax!=nTraceCoefMin: Effeciency may be low ");
5586 
5587  int traceID, nfieldPnts, ElmtId, noffset;
5588  const Array<OneD, const Array<OneD, int > > LocTracephysToTraceIDMap
5589  = locTraceToTraceMap->GetLocTracephysToTraceIDMap();
5590  const Array<OneD, const int > fieldToLocTraceMap =
5591  locTraceToTraceMap->GetfieldToLocTraceMap();
5592  Array<OneD, Array<OneD, int > > fieldToLocTraceMapLR(2);
5593  noffset = 0;
5594  for (int k = 0; k < nFwdBwdNonZero; ++k)
5595  {
5596  size_t i = nlocTracePtsNonZeroIndex[k];
5597  fieldToLocTraceMapLR[i] = Array<OneD, int> (nlocTracePtsLR[i]);
5598  Vmath::Vcopy(nlocTracePtsLR[i],
5599  &fieldToLocTraceMap[0]+noffset,1,
5600  &fieldToLocTraceMapLR[i][0],1);
5601  noffset += nlocTracePtsLR[i];
5602  }
5603 
5604  Array<OneD, Array<OneD, int > > MatIndexArray(2);
5605  for (int k = 0; k < nFwdBwdNonZero; ++k)
5606  {
5607  size_t nlr = nlocTracePtsNonZeroIndex[k];
5608  MatIndexArray[nlr] = Array<OneD, int > (nlocTracePtsLR[nlr]);
5609  for(int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5610  {
5611  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5612  nTraceCoef = TraceCoefArray[traceID];
5613  ElmtId = LRAdjExpid[nlr][traceID];
5614  noffset = GetPhys_Offset(ElmtId);
5615  nElmtCoef = ElmtCoefArray[ElmtId];
5616  nfieldPnts = fieldToLocTraceMapLR[nlr][nloc];
5617  nPnts = nfieldPnts - noffset;
5618 
5619  MatIndexArray[nlr][nloc] = nPnts*nElmtCoef;
5620  }
5621  }
5622 
5623  for(int nc=0;nc<nTraceCoefMin;nc++)
5624  {
5625  for(int nt = 0; nt < ntotTrace; nt++)
5626  {
5627  nTraceCoef = TraceCoefArray[nt];
5628  nTracePnt = TracePntArray[nt] ;
5629  noffset = TraceOffArray[nt] ;
5630  Vmath::Vcopy(nTracePnt,
5631  &FwdMatData[nt][nc], nTraceCoef,
5632  &TraceFwdPhy[noffset],1);
5633  Vmath::Vcopy(nTracePnt,
5634  &BwdMatData[nt][nc],nTraceCoef,
5635  &TraceBwdPhy[noffset],1);
5636  }
5637 
5638  for (int k = 0; k < nFwdBwdNonZero; ++k)
5639  {
5640  size_t i = nlocTracePtsNonZeroIndex[k];
5641  Vmath::Zero(nlocTracePtsLR[i],tmplocTrace[i],1);
5642  }
5643 
5644  GetLocTraceFromTracePts(TraceFwdPhy,TraceBwdPhy,tmplocTrace[0],
5645  tmplocTrace[1]);
5646 
5647  for (int k = 0; k < nFwdBwdNonZero; ++k)
5648  {
5649  size_t nlr = nlocTracePtsNonZeroIndex[k];
5650  for(int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5651  {
5652  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5653  nTraceCoef = TraceCoefArray[traceID];
5654  ElmtId = LRAdjExpid[nlr][traceID];
5655  nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5656  sign = elmtLRSign[nlr][traceID][nc];
5657  MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5658 
5659  ElmtMatDataArray[ElmtId][MatIndex] -=
5660  sign*tmplocTrace[nlr][nloc];
5661  }
5662  }
5663  }
5664 
5665  for(int nc=nTraceCoefMin;nc<nTraceCoefMax;nc++)
5666  {
5667  for(int nt = 0; nt < ntotTrace; nt++)
5668  {
5669  nTraceCoef = TraceCoefArray[nt];
5670  nTracePnt = TracePntArray[nt] ;
5671  noffset = TraceOffArray[nt] ;
5672  if(nc<nTraceCoef)
5673  {
5674  Vmath::Vcopy(nTracePnt,
5675  &FwdMatData[nt][nc],nTraceCoef,
5676  &TraceFwdPhy[noffset],1);
5677  Vmath::Vcopy(nTracePnt,
5678  &BwdMatData[nt][nc],nTraceCoef,
5679  &TraceBwdPhy[noffset],1);
5680  }
5681  else
5682  {
5683  Vmath::Zero(nTracePnt,&TraceFwdPhy[noffset],1);
5684  Vmath::Zero(nTracePnt,&TraceBwdPhy[noffset],1);
5685  }
5686  }
5687 
5688  for (int k = 0; k < nFwdBwdNonZero; ++k)
5689  {
5690  size_t i = nlocTracePtsNonZeroIndex[k];
5691  Vmath::Zero(nlocTracePtsLR[i],tmplocTrace[i],1);
5692  }
5693  GetLocTraceFromTracePts(TraceFwdPhy,TraceBwdPhy,tmplocTrace[0],
5694  tmplocTrace[1]);
5695 
5696  for (int k = 0; k < nFwdBwdNonZero; ++k)
5697  {
5698  size_t nlr = nlocTracePtsNonZeroIndex[k];
5699  for(int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5700  {
5701  traceID = LocTracephysToTraceIDMap[nlr][nloc];
5702  nTraceCoef = TraceCoefArray[traceID];
5703  if(nc<nTraceCoef)
5704  {
5705  ElmtId = LRAdjExpid[nlr][traceID];
5706  nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5707  sign =-elmtLRSign[nlr][traceID][nc];
5708  MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5709 
5710  ElmtMatDataArray[ElmtId][MatIndex] +=
5711  sign*tmplocTrace[nlr][nloc];
5712  }
5713  }
5714  }
5715  }
5716  }
5717 
5718  void ExpList::AddRightIPTPhysDerivBase(
5719  const int dir,
5720  const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
5721  Array<OneD, DNekMatSharedPtr> ElmtJacCoef)
5722  {
5723  int nelmt;
5724  int nelmtcoef, nelmtpnts,nelmtcoef0, nelmtpnts0;
5725 
5726  nelmtcoef = GetNcoeffs(0);
5727  nelmtpnts = GetTotPoints(0);
5728 
5729  Array<OneD,NekDouble> innarray(nelmtpnts,0.0);
5730  Array<OneD,NekDouble> outarray(nelmtcoef,0.0);
5731 
5732  Array<OneD,NekDouble> MatQ_data;
5733  Array<OneD,NekDouble> MatC_data;
5734 
5735  DNekMatSharedPtr tmpMatQ,tmpMatC;
5736 
5737  nelmtcoef0 = nelmtcoef;
5738  nelmtpnts0 = nelmtpnts;
5739 
5740  for(nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5741  {
5742  nelmtcoef = GetNcoeffs(nelmt);
5743  nelmtpnts = GetTotPoints(nelmt);
5744 
5745  tmpMatQ = ElmtJacQuad[nelmt];
5746  tmpMatC = ElmtJacCoef[nelmt];
5747 
5748  MatQ_data = tmpMatQ->GetPtr();
5749  MatC_data = tmpMatC->GetPtr();
5750 
5751  if(nelmtcoef!=nelmtcoef0)
5752  {
5753  outarray = Array<OneD,NekDouble> (nelmtcoef,0.0);
5754  nelmtcoef0 = nelmtcoef;
5755  }
5756 
5757  if(nelmtpnts!=nelmtpnts0)
5758  {
5759  innarray = Array<OneD,NekDouble> (nelmtpnts,0.0);
5760  nelmtpnts0 = nelmtpnts;
5761  }
5762 
5763  for(int np=0; np<nelmtcoef;np++)
5764  {
5765  Vmath::Vcopy(nelmtpnts,&MatQ_data[0]+np,nelmtcoef,&innarray[0],1);
5766  (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray,innarray);
5767  (*m_exp)[nelmt]->IProductWRTDerivBase(dir,innarray,outarray);
5768 
5769  Vmath::Vadd(nelmtcoef,&outarray[0],1,&MatC_data[0]+np,nelmtcoef,&MatC_data[0]+np,nelmtcoef);
5770  }
5771  }
5772  }
5773 
5774 
5775  void ExpList::AddRightIPTBaseMatrix(
5776  const Array<OneD, const DNekMatSharedPtr> ElmtJacQuad,
5777  Array<OneD, DNekMatSharedPtr> ElmtJacCoef)
5778  {
5779  int nelmt;
5780  int nelmtcoef, nelmtpnts,nelmtcoef0, nelmtpnts0;
5781 
5782  nelmtcoef = GetNcoeffs(0);
5783  nelmtpnts = GetTotPoints(0);
5784 
5785  Array<OneD,NekDouble> innarray(nelmtpnts,0.0);
5786  Array<OneD,NekDouble> outarray(nelmtcoef,0.0);
5787 
5788  Array<OneD,NekDouble> MatQ_data;
5789  Array<OneD,NekDouble> MatC_data;
5790 
5791  DNekMatSharedPtr tmpMatQ,tmpMatC;
5792 
5793  nelmtcoef0 = nelmtcoef;
5794  nelmtpnts0 = nelmtpnts;
5795 
5796  for(nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5797  {
5798  nelmtcoef = GetNcoeffs(nelmt);
5799  nelmtpnts = GetTotPoints(nelmt);
5800 
5801  tmpMatQ = ElmtJacQuad[nelmt];
5802  tmpMatC = ElmtJacCoef[nelmt];
5803 
5804  MatQ_data = tmpMatQ->GetPtr();
5805  MatC_data = tmpMatC->GetPtr();
5806 
5807  if(nelmtcoef!=nelmtcoef0)
5808  {
5809  outarray = Array<OneD,NekDouble> (nelmtcoef,0.0);
5810  nelmtcoef0 = nelmtcoef;
5811  }
5812 
5813  if(nelmtpnts!=nelmtpnts0)
5814  {
5815  innarray = Array<OneD,NekDouble> (nelmtpnts,0.0);
5816  nelmtpnts0 = nelmtpnts;
5817  }
5818 
5819  for(int np=0; np<nelmtcoef;np++)
5820  {
5821  Vmath::Vcopy(nelmtpnts,&MatQ_data[0]+np,nelmtcoef,
5822  &innarray[0],1);
5823  (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray,innarray);
5824  (*m_exp)[nelmt]->IProductWRTBase(innarray,outarray);
5825 
5826  Vmath::Vadd(nelmtcoef,&outarray[0],1,
5827  &MatC_data[0]+np,nelmtcoef,
5828  &MatC_data[0]+np,nelmtcoef);
5829  }
5830  }
5831  }
5832 
5833  void ExpList::v_PhysGalerkinProjection1DScaled(
5834  const NekDouble scale,
5835  const Array<OneD, NekDouble> &inarray,
5836  Array<OneD, NekDouble> &outarray)
5837  {
5838  int cnt,cnt1;
5839 
5840  cnt = cnt1 = 0;
5841 
5842  switch(m_expType)
5843  {
5844  case e2D:
5845  {
5846  for(int i = 0; i < GetExpSize(); ++i)
5847  {
5848  // get new points key
5849  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5850  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5851  int npt0 = (int) pt0*scale;
5852  int npt1 = (int) pt1*scale;
5853 
5854  LibUtilities::PointsKey newPointsKey0(npt0,
5855  (*m_exp)[i]->GetPointsType(0));
5856  LibUtilities::PointsKey newPointsKey1(npt1,
5857  (*m_exp)[i]->GetPointsType(1));
5858 
5859  // Project points;
5860  LibUtilities::PhysGalerkinProject2D(newPointsKey0,
5861  newPointsKey1,
5862  &inarray[cnt],
5863  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5864  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5865  &outarray[cnt1]);
5866 
5867  cnt += npt0*npt1;
5868  cnt1 += pt0*pt1;
5869  }
5870  }
5871  break;
5872  case e3D:
5873  {
5874  for(int i = 0; i < GetExpSize(); ++i)
5875  {
5876  // get new points key
5877  int pt0 = (*m_exp)[i]->GetNumPoints(0);
5878  int pt1 = (*m_exp)[i]->GetNumPoints(1);
5879  int pt2 = (*m_exp)[i]->GetNumPoints(2);
5880  int npt0 = (int) pt0*scale;
5881  int npt1 = (int) pt1*scale;
5882  int npt2 = (int) pt2*scale;
5883 
5884  LibUtilities::PointsKey newPointsKey0(npt0,
5885  (*m_exp)[i]->GetPointsType(0));
5886  LibUtilities::PointsKey newPointsKey1(npt1,
5887  (*m_exp)[i]->GetPointsType(1));
5888  LibUtilities::PointsKey newPointsKey2(npt2,
5889  (*m_exp)[i]->GetPointsType(2));
5890 
5891  // Project points;
5892  LibUtilities::PhysGalerkinProject3D(newPointsKey0,
5893  newPointsKey1,
5894  newPointsKey2,
5895  &inarray[cnt],
5896  (*m_exp)[i]->GetBasis(0)->GetPointsKey(),
5897  (*m_exp)[i]->GetBasis(1)->GetPointsKey(),
5898  (*m_exp)[i]->GetBasis(2)->GetPointsKey(),
5899  &outarray[cnt1]);
5900 
5901  cnt += npt0*npt1*npt2;
5902  cnt1 += pt0*pt1*pt2;
5903  }
5904  }
5905  break;
5906  default:
5907  {
5908  NEKERROR(ErrorUtil::efatal,"not setup for this expansion");
5909  }
5910  break;
5911  }
5912  }
5913 
5914  const LocTraceToTraceMapSharedPtr
5915  &ExpList::v_GetLocTraceToTraceMap() const
5916  {
5917  NEKERROR(ErrorUtil::efatal, "v_GetLocTraceToTraceMap not coded");
5918  return NullLocTraceToTraceMapSharedPtr;
5919  }
5920 
5921  } //end of namespace
5922 } //end of namespace
5923 
WARNINGL1
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:251
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:250
LIKWID_MARKER_START
#define LIKWID_MARKER_START(regionTag)
Definition: Likwid.hpp:46
LIKWID_MARKER_STOP
#define LIKWID_MARKER_STOP(regionTag)
Definition: Likwid.hpp:47
sign
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Nektar::Collections::CollectionOptimisation::SetWithTimings
COLLECTIONS_EXPORT OperatorImpMap SetWithTimings(std::vector< StdRegions::StdExpansionSharedPtr > pGeom, OperatorImpMap &impTypes, bool verbose=true)
Definition: CollectionOptimisation.cpp:321
Nektar::Collections::CollectionOptimisation::GetOperatorImpMap
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Definition: CollectionOptimisation.cpp:283
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::BasisKey::GetNumPoints
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:133
Nektar::LibUtilities::BasisKey::GetBasisType
BasisType GetBasisType() const
Return type of expansion basis.
Definition: Basis.h:144
Nektar::LibUtilities::BasisKey::GetPointsKey
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:150
Nektar::LibUtilities::BasisKey::GetNumModes
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:86
Nektar::LibUtilities::FieldIO::GetFileType
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
Definition: FieldIO.cpp:97
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:145
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:60
Nektar::LibUtilities::Timer::AccumulateRegion
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
Definition: Timer.cpp:67
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:84
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:171
Nektar::MultiRegions::ExpList
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:107
Nektar::MultiRegions::ExpList::v_Reset
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:2580
Nektar::MultiRegions::ExpList::m_coeffs
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1252
Nektar::MultiRegions::ExpList::AddTraceJacToElmtJac
void AddTraceJacToElmtJac(const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian,...
Definition: ExpList.cpp:5477
Nektar::MultiRegions::ExpList::GetBndCondExpansions
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions()
Definition: ExpList.h:2470
Nektar::MultiRegions::ExpList::IProductWRTBase_IterPerExp
void IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
Definition: ExpList.h:1966
Nektar::MultiRegions::ExpList::GetTrace
std::shared_ptr< ExpList > & GetTrace()
Definition: ExpList.h:2525
Nektar::MultiRegions::ExpList::GetNcoeffs
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
Definition: ExpList.h:1800
Nektar::MultiRegions::ExpList::v_ExtractCoeffsToCoeffs
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Definition: ExpList.cpp:3546
Nektar::MultiRegions::ExpList::GetBlockMatrix
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:1946
Nektar::MultiRegions::ExpList::v_AddFwdBwdTraceIntegral
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4283
Nektar::MultiRegions::ExpList::ExpList
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
Definition: ExpList.cpp:104
Nektar::MultiRegions::ExpList::ExtractDataToCoeffs
void ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract the data in fielddata into the coeffs.
Definition: ExpList.cpp:3435
Nektar::MultiRegions::ExpList::CreateCollections
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
Definition: ExpList.cpp:5017
Nektar::MultiRegions::ExpList::v_FwdTrans_BndConstrained
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1745
Nektar::MultiRegions::ExpList::v_GetBoundaryToElmtMap
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.cpp:4902
Nektar::MultiRegions::ExpList::ExponentialFilter
void ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: ExpList.cpp:1670
Nektar::MultiRegions::ExpList::v_HomogeneousBwdTrans
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:4456
Nektar::MultiRegions::ExpList::GenGlobalMatrix
std::shared_ptr< GlobalMatrix > GenGlobalMatrix(const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Generates a global matrix from the given key and map.
Definition: ExpList.cpp:2034
Nektar::MultiRegions::ExpList::v_HomogeneousFwdTrans
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:4446
Nektar::MultiRegions::ExpList::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4620
Nektar::MultiRegions::ExpList::AddRightIPTPhysDerivBase
void AddRightIPTPhysDerivBase(const int dir, const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:5718
Nektar::MultiRegions::ExpList::v_ExtractTracePhys
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4383
Nektar::MultiRegions::ExpList::GetExpIndex
int GetExpIndex(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
This function returns the index of the local elemental expansion containing the arbitrary point given...
Definition: ExpList.cpp:2418
Nektar::MultiRegions::ExpList::v_WriteTecplotConnectivity
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2741
Nektar::MultiRegions::ExpList::m_blockMat
BlockMatrixMapShPtr m_blockMat
Definition: ExpList.h:1312
Nektar::MultiRegions::ExpList::v_GetYIDs
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
Definition: ExpList.cpp:3192
Nektar::MultiRegions::ExpList::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
multiply the metric jacobi and quadrature weights
Definition: ExpList.cpp:1220
Nektar::MultiRegions::ExpList::v_LinearAdvectionReactionSolve
virtual void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:4434
Nektar::MultiRegions::ExpList::GetBoundaryToElmtMap
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2709
Nektar::MultiRegions::ExpList::WriteVtkPieceFooter
void WriteVtkPieceFooter(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2991
Nektar::MultiRegions::ExpList::v_ExtractPhysToBnd
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
Definition: ExpList.cpp:4819
Nektar::MultiRegions::ExpList::GetBwdWeight
void GetBwdWeight(Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
Get the weight value for boundary conditions for boundary average and jump calculations.
Definition: ExpList.cpp:3572
Nektar::MultiRegions::ExpList::v_AddTraceIntegralToOffDiag
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5244
Nektar::MultiRegions::ExpList::InitialiseExpVector
void InitialiseExpVector(const SpatialDomains::ExpansionInfoMap &expmap)
Define a list of elements using the geometry and basis key information in expmap;.
Definition: ExpList.cpp:1022
Nektar::MultiRegions::ExpList::SetupCoeffPhys
void SetupCoeffPhys(bool DeclareCoeffPhysArrays=true, bool SetupOffsets=true)
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
Definition: ExpList.cpp:978
Nektar::MultiRegions::ExpList::GetBoundaryCondition
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:4983
Nektar::MultiRegions::ExpList::v_PhysDirectionalDeriv
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1638
Nektar::MultiRegions::ExpList::v_GetBCValues
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
Definition: ExpList.cpp:4485
Nektar::MultiRegions::ExpList::v_GetLocTraceFromTracePts
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.cpp:4345
Nektar::MultiRegions::ExpList::v_WriteTecplotHeader
virtual void v_WriteTecplotHeader(std::ostream &outfile, std::string var="")
Definition: ExpList.cpp:2607
Nektar::MultiRegions::ExpList::v_GetPlane
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
Definition: ExpList.cpp:5005
Nektar::MultiRegions::ExpList::GetCoeff_Offset
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n.
Definition: ExpList.h:2431
Nektar::MultiRegions::ExpList::GetExpSize
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2401
Nektar::MultiRegions::ExpList::v_AddTraceQuadPhysToOffDiag
virtual void v_AddTraceQuadPhysToOffDiag(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:4335
Nektar::MultiRegions::ExpList::v_ClearGlobalLinSysManager
virtual void v_ClearGlobalLinSysManager(void)
Definition: ExpList.cpp:3201
Nektar::MultiRegions::ExpList::GeneralGetFieldDefinitions
void GeneralGetFieldDefinitions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, int NumHomoDir=0, Array< OneD, LibUtilities::BasisSharedPtr > &HomoBasis=LibUtilities::NullBasisSharedPtr1DArray, std::vector< NekDouble > &HomoLen=LibUtilities::NullNekDoubleVector, bool homoStrips=false, std::vector< unsigned int > &HomoSIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoZIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoYIDs=LibUtilities::NullUnsignedIntVector)
Definition: ExpList.cpp:3286
Nektar::MultiRegions::ExpList::v_GetPeriodicEntities
virtual void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
Definition: ExpList.cpp:4973
Nektar::MultiRegions::ExpList::IProductWRTDerivBase
void IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to the derivative (in directio...
Definition: ExpList.cpp:1299
Nektar::MultiRegions::ExpList::v_CurlCurl
virtual void v_CurlCurl(Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
Definition: ExpList.cpp:1566
Nektar::MultiRegions::ExpList::v_LocalToGlobal
virtual void v_LocalToGlobal(bool UseComm)
Definition: ExpList.cpp:4574
Nektar::MultiRegions::ExpList::GetNumElmts
int GetNumElmts(void)
This function returns the number of elements in the expansion which may be different for a homogeoeno...
Definition: ExpList.h:743
Nektar::MultiRegions::ExpList::v_GetBndConditions
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions()
Definition: ExpList.cpp:4929
Nektar::MultiRegions::ExpList::~ExpList
virtual ~ExpList()
The default destructor.
Definition: ExpList.cpp:1187
Nektar::MultiRegions::ExpList::v_FillBndCondFromField
virtual void v_FillBndCondFromField()
Definition: ExpList.cpp:4559
Nektar::MultiRegions::ExpList::v_ImposeDirichletConditions
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4550
Nektar::MultiRegions::ExpList::v_EvaluateBoundaryConditions
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
Definition: ExpList.cpp:4950
Nektar::MultiRegions::ExpList::v_GetLeftAdjacentFaces
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
Definition: ExpList.cpp:4374
Nektar::MultiRegions::ExpList::m_coeff_offset
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1304
Nektar::MultiRegions::ExpList::GetTraceMap
std::shared_ptr< AssemblyMapDG > & GetTraceMap(void)
Definition: ExpList.h:2530
Nektar::MultiRegions::ExpList::v_GetMovingFrames
virtual void v_GetMovingFrames(const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ExpList.cpp:3586
Nektar::MultiRegions::ExpList::v_GetTraceMap
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap()
Definition: ExpList.cpp:3805
Nektar::MultiRegions::ExpList::m_coeffsToElmt
Array< OneD, std::pair< int, int > > m_coeffsToElmt
m_coeffs to elemental value map
Definition: ExpList.h:1310
Nektar::MultiRegions::ExpList::v_GetFwdBwdTracePhys
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:4293
Nektar::MultiRegions::ExpList::PhysEvaluate
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
Definition: ExpList.cpp:2542
Nektar::MultiRegions::ExpList::v_Integral
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Definition: ExpList.cpp:3120
Nektar::MultiRegions::ExpList::v_GetRobinBCInfo
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo(void)
Definition: ExpList.cpp:4963
Nektar::MultiRegions::ExpList::v_GetTraceBndMap
virtual const Array< OneD, const int > & v_GetTraceBndMap()
Definition: ExpList.cpp:3813
Nektar::MultiRegions::ExpList::v_GetBoundaryNormals
virtual void v_GetBoundaryNormals(int i, Array< OneD, Array< OneD, NekDouble > > &normals)
Definition: ExpList.cpp:4857
Nektar::MultiRegions::ExpList::v_LinearAdvectionDiffusionReactionSolve
virtual void v_LinearAdvectionDiffusionReactionSolve(const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:4422
Nektar::MultiRegions::ExpList::v_HomogeneousEnergy
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
Definition: ExpList.cpp:3153
Nektar::MultiRegions::ExpList::v_GetHomoLen
virtual NekDouble v_GetHomoLen(void)
Definition: ExpList.cpp:3169
Nektar::MultiRegions::ExpList::ExtractFileBCs
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:3207
Nektar::MultiRegions::ExpList::GetLocTraceToTraceMap
const LocTraceToTraceMapSharedPtr & GetLocTraceToTraceMap() const
Definition: ExpList.h:2484
Nektar::MultiRegions::ExpList::v_GetZIDs
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
Definition: ExpList.cpp:3184
Nektar::MultiRegions::ExpList::ExtractCoeffsToCoeffs
void ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Extract the data from fromField using fromExpList the coeffs using the basic ExpList Elemental expans...
Definition: ExpList.cpp:3444
Nektar::MultiRegions::ExpList::GetCoords
void GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
This function calculates the coordinates of all the elemental quadrature points .
Definition: ExpList.h:2089
Nektar::MultiRegions::ExpList::GenGlobalBndLinSys
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
Definition: ExpList.cpp:2329
Nektar::MultiRegions::ExpList::v_SetBndCondBwdWeight
virtual void v_SetBndCondBwdWeight(const int index, const NekDouble value)
Definition: ExpList.cpp:4365
Nektar::MultiRegions::ExpList::v_SmoothField
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:1767
Nektar::MultiRegions::ExpList::v_AppendFieldData
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
Definition: ExpList.cpp:3407
Nektar::MultiRegions::ExpList::m_collectionsDoInit
std::vector< bool > m_collectionsDoInit
Vector of bools to act as an initialise on first call flag.
Definition: ExpList.h:1295
Nektar::MultiRegions::ExpList::v_FwdTrans_IterPerExp
virtual void v_FwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1735
Nektar::MultiRegions::ExpList::m_physState
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1278
Nektar::MultiRegions::ExpList::AddRightIPTBaseMatrix
void AddRightIPTBaseMatrix(const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
Definition: ExpList.cpp:5775
Nektar::MultiRegions::ExpList::GetLocTraceFromTracePts
void GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
Definition: ExpList.h:2604
Nektar::MultiRegions::ExpList::v_GeneralMatrixOp
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4644
Nektar::MultiRegions::ExpList::v_WriteVtkPieceData
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
Definition: ExpList.cpp:2998
Nektar::MultiRegions::ExpList::v_WriteVtkPieceHeader
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Definition: ExpList.cpp:2866
Nektar::MultiRegions::ExpList::v_BwdTrans_IterPerExp
virtual void v_BwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2366
Nektar::MultiRegions::ExpList::GetShapeDimension
int GetShapeDimension()
This function returns the dimension of the shape of the element eid.
Definition: ExpList.h:2255
Nektar::MultiRegions::ExpList::H1
NekDouble H1(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:3267
Nektar::MultiRegions::ExpList::m_coll_phys_offset
std::vector< int > m_coll_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1301
Nektar::MultiRegions::ExpList::m_comm
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1220
Nektar::MultiRegions::ExpList::GenGlobalLinSys
std::shared_ptr< GlobalLinSys > GenGlobalLinSys(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
This operation constructs the global linear system of type mkey.
Definition: ExpList.cpp:2310
Nektar::MultiRegions::ExpList::WriteVtkHeader
void WriteVtkHeader(std::ostream &outfile)
Definition: ExpList.cpp:2852
Nektar::MultiRegions::ExpList::m_exp
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1290
Nektar::MultiRegions::ExpList::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
Definition: ExpList.cpp:4663
Nektar::MultiRegions::ExpList::m_ncoeffs
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1230
Nektar::MultiRegions::ExpList::GenerateElementVector
void GenerateElementVector(const int ElementID, const NekDouble scalar1, const NekDouble scalar2, Array< OneD, NekDouble > &outarray)
Generate vector v such that v[i] = scalar1 if i is in the element < ElementID. Otherwise,...
Definition: ExpList.cpp:3638
Nektar::MultiRegions::ExpList::GenGlobalMatrixFull
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Definition: ExpList.cpp:2170
Nektar::MultiRegions::ExpList::v_NormVectorIProductWRTBase
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
Definition: ExpList.cpp:4494
Nektar::MultiRegions::ExpList::v_ExtractElmtToBndPhys
virtual void v_ExtractElmtToBndPhys(const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
Definition: ExpList.cpp:4740
Nektar::MultiRegions::ExpList::v_ExtractPhysToBndElmt
virtual void v_ExtractPhysToBndElmt(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
Definition: ExpList.cpp:4779
Nektar::MultiRegions::ExpList::v_Upwind
virtual void v_Upwind(const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.cpp:3696
Nektar::MultiRegions::ExpList::v_AddTraceIntegral
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4264
Nektar::MultiRegions::ExpList::v_FillBwdWithBwdWeight
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Definition: ExpList.cpp:4910
Nektar::MultiRegions::ExpList::v_PhysGalerkinProjection1DScaled
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5833
Nektar::MultiRegions::ExpList::v_SetUpPhysNormals
virtual void v_SetUpPhysNormals()
Definition: ExpList.cpp:4716
Nektar::MultiRegions::ExpList::v_GlobalToLocal
virtual void v_GlobalToLocal(void)
Definition: ExpList.cpp:4592
Nektar::MultiRegions::ExpList::v_GetTransposition
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
Definition: ExpList.cpp:3161
Nektar::MultiRegions::ExpList::GeneralMatrixOp_IterPerExp
void GeneralMatrixOp_IterPerExp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1961
Nektar::MultiRegions::ExpList::GetExp
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2422
Nektar::MultiRegions::ExpList::v_HelmSolve
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Definition: ExpList.cpp:4408
Nektar::MultiRegions::ExpList::m_graph
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1226
Nektar::MultiRegions::ExpList::v_MultiplyByInvMassMatrix
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4399
Nektar::MultiRegions::ExpList::v_AddTraceQuadPhysToField
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:4325
Nektar::MultiRegions::ExpList::v_GetBndCondBwdWeight
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight()
Definition: ExpList.cpp:4357
Nektar::MultiRegions::ExpList::v_L2
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Definition: ExpList.cpp:3074
Nektar::MultiRegions::ExpList::WriteVtkFooter
void WriteVtkFooter(std::ostream &outfile)
Definition: ExpList.cpp:2860
Nektar::MultiRegions::ExpList::MultiplyByBlockMatrix
void MultiplyByBlockMatrix(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1199
Nektar::MultiRegions::ExpList::v_GetTrace
virtual std::shared_ptr< ExpList > & v_GetTrace()
Definition: ExpList.cpp:3797
Nektar::MultiRegions::ExpList::GetElmtNormalLength
void GetElmtNormalLength(Array< OneD, NekDouble > &lengthsFwd, Array< OneD, NekDouble > &lengthsBwd)
Get the length of elements in boundary normal direction.
Definition: ExpList.cpp:4116
Nektar::MultiRegions::ExpList::v_VectorFlux
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Definition: ExpList.cpp:3134
Nektar::MultiRegions::ExpList::v_GetFieldDefinitions
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
Definition: ExpList.cpp:3393
Nektar::MultiRegions::ExpList::m_coll_coeff_offset
std::vector< int > m_coll_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1298
Nektar::MultiRegions::ExpList::IProductWRTDirectionalDerivBase
void IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directional derivative along a given direction.
Definition: ExpList.cpp:1319
Nektar::MultiRegions::ExpList::FillBwdWithBwdWeight
void FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill Bwd with boundary conditions.
Definition: ExpList.h:2613
Nektar::MultiRegions::ExpList::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4608
Nektar::MultiRegions::ExpList::v_UpdateBndConditions
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
Definition: ExpList.cpp:4940
Nektar::MultiRegions::ExpList::GetMatIpwrtDeriveBase
void GetMatIpwrtDeriveBase(const Array< OneD, const Array< OneD, NekDouble > > &inarray, const int nDirctn, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:5253
Nektar::MultiRegions::ExpList::v_WriteTecplotField
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2818
Nektar::MultiRegions::ExpList::GetSharedThisPtr
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: ExpList.h:1129
Nektar::MultiRegions::ExpList::m_collections
Collections::CollectionVector m_collections
Definition: ExpList.h:1292
Nektar::MultiRegions::ExpList::v_PhysInterp1DScaled
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:5170
Nektar::MultiRegions::ExpList::ApplyGeomInfo
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:2567
Nektar::MultiRegions::ExpList::GetZIDs
Array< OneD, const unsigned int > GetZIDs(void)
This function returns a vector containing the wave numbers in z-direction associated with the 3D homo...
Definition: ExpList.h:677
Nektar::MultiRegions::ExpList::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: ExpList.cpp:1465
Nektar::MultiRegions::ExpList::m_elmtToExpId
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1321
Nektar::MultiRegions::ExpList::v_WriteTecplotZone
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:2646
Nektar::MultiRegions::ExpList::m_session
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1223
Nektar::MultiRegions::ExpList::v_DealiasedDotProd
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ExpList.cpp:4475
Nektar::MultiRegions::ExpList::v_PeriodicBwdCopy
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:4918
Nektar::MultiRegions::ExpList::GetRobinBCInfo
std::map< int, RobinBCInfoSharedPtr > GetRobinBCInfo()
Definition: ExpList.h:1048
Nektar::MultiRegions::ExpList::v_FillBwdWithBoundCond
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
Definition: ExpList.cpp:4315
Nektar::MultiRegions::ExpList::DivideByQuadratureMetric
void DivideByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Divided by the metric jacobi and quadrature weights.
Definition: ExpList.cpp:1237
Nektar::MultiRegions::ExpList::v_GetBndCondExpansions
virtual const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions(void)
Definition: ExpList.cpp:3669
Nektar::MultiRegions::ExpList::m_phys_offset
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1307
Nektar::MultiRegions::ExpList::v_SetHomoLen
virtual void v_SetHomoLen(const NekDouble lhom)
Definition: ExpList.cpp:3177
Nektar::MultiRegions::ExpList::Upwind
void Upwind(const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Definition: ExpList.h:2516
Nektar::MultiRegions::ExpList::v_GetNormals
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Populate normals with the normals of all expansions.
Definition: ExpList.cpp:3893
Nektar::MultiRegions::ExpList::GenBlockMatrix
const DNekScalBlkMatSharedPtr GenBlockMatrix(const GlobalMatrixKey &gkey)
This function assembles the block diagonal matrix of local matrices of the type mtype.
Definition: ExpList.cpp:1810
Nektar::MultiRegions::ExpList::v_IProductWRTBase_IterPerExp
virtual void v_IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1262
Nektar::MultiRegions::ExpList::PhysDeriv
void PhysDeriv(Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
Definition: ExpList.h:2132
Nektar::MultiRegions::ExpList::MultiplyByElmtInvMass
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally mulplies the coefficient space of Sin my the elemental inverse of the mass ...
Definition: ExpList.cpp:1696
Nektar::MultiRegions::ExpList::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.h:2015
Nektar::MultiRegions::ExpList::v_GetBndElmtExpansion
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Definition: ExpList.cpp:4729
Nektar::MultiRegions::ExpList::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4614
Nektar::MultiRegions::ExpList::v_UpdateBndCondExpansion
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
Definition: ExpList.cpp:3677
Nektar::MultiRegions::ExpList::GetPhys_Offset
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n.
Definition: ExpList.h:2439
Nektar::MultiRegions::ExpList::v_GetLocTraceToTraceMap
virtual const std::shared_ptr< LocTraceToTraceMap > & v_GetLocTraceToTraceMap(void) const
Definition: ExpList.cpp:5915
Nektar::MultiRegions::ExpList::v_ExtractDataToCoeffs
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
Definition: ExpList.cpp:3457
Nektar::MultiRegions::ExpList::m_expType
ExpansionType m_expType
Exapnsion type.
Definition: ExpList.h:1212
Nektar::MultiRegions::ExpList::m_phys
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1269
Nektar::MultiRegions::ExpList::GetExpType
ExpansionType GetExpType(void)
Returns the type of the expansion.
Definition: ExpList.cpp:1182
Nektar::MultiRegions::ExpList::Linf
NekDouble Linf(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global spectral/hp element approximation.
Definition: ExpList.cpp:3035
Nektar::MultiRegions::ExpList::v_DealiasedProd
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:4466
Nektar::MultiRegions::ExpList::GetTotPoints
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1850
Nektar::MultiRegions::ExpList::GetCoordim
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:2244
Nektar::MultiRegions::ExpList::GetDiagMatIpwrtBase
void GetDiagMatIpwrtBase(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
Definition: ExpList.cpp:5405
Nektar::MultiRegions::GlobalLinSysKey::GetGlobalSysSolnType
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Definition: GlobalLinSysKey.h:98
Nektar::MultiRegions::GlobalMatrixKey
Describes a matrix with ordering defined by a local to global map.
Definition: GlobalMatrixKey.h:48
Nektar::MultiRegions::GlobalMatrixKey::GetConstFactors
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
Definition: GlobalMatrixKey.h:149
Nektar::MultiRegions::GlobalMatrixKey::GetShapeType
LibUtilities::ShapeType GetShapeType() const
Return the expansion type associated with key.
Definition: GlobalMatrixKey.h:121
Nektar::MultiRegions::GlobalMatrixKey::GetVarCoeffs
const StdRegions::VarCoeffMap & GetVarCoeffs() const
Definition: GlobalMatrixKey.h:165
Nektar::MultiRegions::GlobalMatrixKey::GetMatrixType
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
Definition: GlobalMatrixKey.h:115
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:182
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::Collections::OperatorImpMap
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:108
Nektar::LibUtilities::PhysGalerkinProject3D
void PhysGalerkinProject3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
Definition: PhysGalerkinProject.cpp:177
Nektar::LibUtilities::FieldIOSharedPtr
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:306
Nektar::LibUtilities::NullBasisKey
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
Nektar::LibUtilities::Interp1D
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:53
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:120
Nektar::LibUtilities::Interp3D
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
Definition: Interp.cpp:185
Nektar::LibUtilities::FieldDefinitionsSharedPtr
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:179
Nektar::LibUtilities::Interp2D
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:115
Nektar::LibUtilities::GetFieldIOFactory
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:72
Nektar::LibUtilities::TranspositionSharedPtr
std::shared_ptr< Transposition > TranspositionSharedPtr
Definition: Transposition.h:165
Nektar::LibUtilities::eNodalTriElec
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
Definition: PointsType.h:69
Nektar::LibUtilities::CommSharedPtr
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:54
Nektar::LibUtilities::PhysGalerkinProject2D
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: PhysGalerkinProject.cpp:115
Nektar::LibUtilities::GetNumberOfCoefficients
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset)
Definition: ShapeType.hpp:313
Nektar::LibUtilities::eGauss_Lagrange
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:55
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:45
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:54
Nektar::LocalRegions::PointExpSharedPtr
std::shared_ptr< PointExp > PointExpSharedPtr
Definition: PointExp.h:129
Nektar::LocalRegions::ExpansionSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:68
Nektar::LocalRegions::Expansion2DSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
Definition: Expansion1D.h:47
Nektar::LocalRegions::Expansion0DSharedPtr
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
Definition: Expansion0D.h:48
Nektar::LocalRegions::ExpansionVector
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:70
Nektar::LocalRegions::Expansion1DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
Definition: Expansion1D.h:51
Nektar::LocalRegions::Expansion3DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
Definition: Expansion2D.h:49
Nektar::MultiRegions::DirCartesianMap
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:90
Nektar::MultiRegions::AlignFace
void AlignFace(const StdRegions::Orientation orient, const int nquad1, const int nquad2, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out)
Helper function to re-align face to a given orientation.
Definition: ExpList.cpp:3821
Nektar::MultiRegions::NullExpListSharedPtr
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1792
Nektar::MultiRegions::AssemblyMapCGSharedPtr
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:51
Nektar::MultiRegions::RobinBCInfoSharedPtr
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
Definition: MultiRegions.hpp:187
Nektar::MultiRegions::GlobalSysSolnTypeMap
const char *const GlobalSysSolnTypeMap[]
Definition: MultiRegions.hpp:87
Nektar::MultiRegions::GlobalLinSysSharedPtr
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
Nektar::MultiRegions::GlobalMatrixSharedPtr
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
Definition: GlobalMatrix.h:88
Nektar::MultiRegions::GetGlobalLinSysFactory
GlobalLinSysFactory & GetGlobalLinSysFactory()
Definition: GlobalLinSys.cpp:212
Nektar::MultiRegions::NullLocTraceToTraceMapSharedPtr
static LocTraceToTraceMapSharedPtr NullLocTraceToTraceMapSharedPtr
Definition: LocTraceToTraceMap.h:383
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::MultiRegions::LocTraceToTraceMapSharedPtr
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
Definition: LocTraceToTraceMap.h:381
Nektar::MultiRegions::VarFactorsMap
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
Definition: GlobalLinSysKey.h:46
Nektar::MultiRegions::AssemblyMapSharedPtr
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:52
Nektar::MultiRegions::PeriodicMap
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
Definition: MultiRegions.hpp:207
Nektar::MultiRegions::BlockMatrixMap
std::map< GlobalMatrixKey, DNekScalBlkMatSharedPtr > BlockMatrixMap
A map between global matrix keys and their associated block matrices.
Definition: ExpList.h:99
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::SpatialDomains::QuadGeomSharedPtr
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
Definition: HexGeom.h:46
Nektar::SpatialDomains::GeomMMF
GeomMMF
Principle direction for MMF.
Definition: SpatialDomains.hpp:74
Nektar::SpatialDomains::BoundaryConditionShPtr
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
Nektar::SpatialDomains::PrismGeomSharedPtr
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
Definition: PrismGeom.h:82
Nektar::SpatialDomains::CompositeMap
std::map< int, CompositeSharedPtr > CompositeMap
Definition: MeshGraph.h:138
Nektar::SpatialDomains::BoundaryConditionCollection
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:226
Nektar::SpatialDomains::eDeformed
@ eDeformed
Geometry is curved or has non-constant factors.
Definition: SpatialDomains.hpp:56
Nektar::SpatialDomains::SegGeomSharedPtr
std::shared_ptr< SegGeom > SegGeomSharedPtr
Definition: Geometry2D.h:62
Nektar::SpatialDomains::PyrGeomSharedPtr
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
Definition: PyrGeom.h:74
Nektar::SpatialDomains::ExpansionInfoShPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
Definition: MeshGraph.h:140
Nektar::SpatialDomains::TetGeomSharedPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
Nektar::SpatialDomains::MeshGraphSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
Nektar::SpatialDomains::BoundaryConditionMapShPtr
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:225
Nektar::SpatialDomains::HexGeomSharedPtr
std::shared_ptr< HexGeom > HexGeomSharedPtr
Definition: HexGeom.h:84
Nektar::SpatialDomains::PointGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
Nektar::SpatialDomains::Geometry2DSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
Nektar::SpatialDomains::TriGeomSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
Definition: TriGeom.h:58
Nektar::SpatialDomains::Geometry1DSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:63
Nektar::SpatialDomains::ExpansionInfoMap
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
Definition: MeshGraph.h:143
Nektar::StdRegions::StdExpansionSharedPtr
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
Definition: StdExpansion.h:1714
Nektar::StdRegions::VarCoeffMap
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:272
Nektar::StdRegions::ConstFactorMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:314
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::DNekScalMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
Definition: NekMatrixFwd.hpp:68
Nektar::DNekScalMatSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
Definition: NekMatrixFwd.hpp:69
Nektar::CoordType
std::pair< IndexType, IndexType > CoordType
Definition: SparseMatrixFwd.hpp:55
Nektar::DNekScalBlkMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
Definition: MatrixStorageType.h:51
Nektar::ePOSITIVE_DEFINITE_SYMMETRIC
@ ePOSITIVE_DEFINITE_SYMMETRIC
Definition: MatrixStorageType.h:48
Nektar::COOMatType
std::map< CoordType, NekDouble > COOMatType
Definition: SparseMatrixFwd.hpp:57
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:820
Nektar::DNekScalBlkMat
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:65
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:565
Vmath::Neg
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:461
Vmath::Vsum
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.cpp:846
Vmath::Scatr
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
Definition: Vmath.cpp:772
Vmath::Vadd
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Vmax
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
Definition: Vmath.cpp:892
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
Vmath::Vsub
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:372
tinysimd::abs
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:272
tinysimd::sqrt
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:267
Nektar::LocalRegions::MatrixKey::opLess
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:68
Nektar::MultiRegions::RobinBCInfo::next
std::shared_ptr< RobinBCInfo > next
Definition: MultiRegions.hpp:184
Generated on Wed Nov 24 2021 10:33:26 for Nektar++ by doxygen 1.9.1