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 <boost/core/ignore_unused.hpp>
36 
37 #include <MultiRegions/ExpList.h>
40 
41 #include <StdRegions/StdSegExp.h>
42 #include <StdRegions/StdTriExp.h>
44 #include <StdRegions/StdQuadExp.h>
45 #include <StdRegions/StdTetExp.h>
46 #include <StdRegions/StdPyrExp.h>
47 #include <StdRegions/StdPrismExp.h>
48 #include <StdRegions/StdHexExp.h>
49 
50 #include <LocalRegions/MatrixKey.h> // for MatrixKey
51 #include <LocalRegions/Expansion.h> // for Expansion
52 
53 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h> // for AssemblyMapCG, etc
54 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h> // for AssemblyMapCG, etc
55 #include <MultiRegions/GlobalLinSysKey.h> // for GlobalLinSysKey
56 #include <MultiRegions/GlobalMatrix.h> // for GlobalMatrix, etc
57 #include <MultiRegions/GlobalMatrixKey.h> // for GlobalMatrixKey
58 
62 
64 #include <Collections/Operator.h>
65 
66 using namespace std;
67 
68 namespace Nektar
69 {
70  namespace MultiRegions
71  {
72  /**
73  * @class ExpList
74  * All multi-elemental expansions \f$u^{\delta}(\boldsymbol{x})\f$ can
75  * be considered as the assembly of the various elemental contributions.
76  * On a discrete level, this yields,
77  * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
78  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
79  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
80  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
81  * As it is the lowest level class, it contains the definition of the
82  * common data and common routines to all multi-elemental expansions.
83  *
84  * The class stores a vector of expansions, \a m_exp, (each derived from
85  * StdRegions#StdExpansion) which define the constituent components of
86  * the domain. The coefficients from these expansions are concatenated
87  * in \a m_coeffs, while the expansion evaluated at the quadrature
88  * points is stored in \a m_phys.
89  */
90 
91  /**
92  * Creates an empty expansion list. The expansion list will typically be
93  * populated by a derived class (namely one of MultiRegions#ExpList1D,
94  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
95  */
96  ExpList::ExpList():
97  m_comm(),
98  m_session(),
99  m_graph(),
100  m_ncoeffs(0),
101  m_npoints(0),
102  m_coeffs(),
103  m_phys(),
104  m_physState(false),
105  m_exp(MemoryManager<LocalRegions::ExpansionVector>
106  ::AllocateSharedPtr()),
107  m_coeff_offset(),
108  m_phys_offset(),
109  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
110  m_WaveSpace(false)
111  {
113  }
114 
115 
116  /**
117  * Creates an empty expansion list. The expansion list will typically be
118  * populated by a derived class (namely one of MultiRegions#ExpList1D,
119  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
120  */
122  const LibUtilities::SessionReaderSharedPtr &pSession):
123  m_comm(pSession->GetComm()),
124  m_session(pSession),
125  m_graph(),
126  m_ncoeffs(0),
127  m_npoints(0),
128  m_coeffs(),
129  m_phys(),
130  m_physState(false),
131  m_exp(MemoryManager<LocalRegions::ExpansionVector>
132  ::AllocateSharedPtr()),
133  m_coeff_offset(),
134  m_phys_offset(),
135  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
136  m_WaveSpace(false)
137  {
139  }
140 
141 
142  /**
143  * Creates an empty expansion list. The expansion list will typically be
144  * populated by a derived class (namely one of MultiRegions#ExpList1D,
145  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
146  */
148  const LibUtilities::SessionReaderSharedPtr &pSession,
149  const SpatialDomains::MeshGraphSharedPtr &pGraph):
150  m_comm(pSession->GetComm()),
151  m_session(pSession),
152  m_graph(pGraph),
153  m_ncoeffs(0),
154  m_npoints(0),
155  m_coeffs(),
156  m_phys(),
157  m_physState(false),
158  m_exp(MemoryManager<LocalRegions::ExpansionVector>
159  ::AllocateSharedPtr()),
160  m_coeff_offset(),
161  m_phys_offset(),
162  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
163  m_WaveSpace(false)
164  {
166  }
167 
168  /**
169  * Copies the eIds elements from an existing expansion list.
170  * @param in Source expansion list.
171  * @param in elements that will be in the new exp list.
172  */
174  const std::vector<unsigned int> &eIDs,
175  const bool DeclareCoeffPhysArrays):
176  m_comm(in.m_comm),
177  m_session(in.m_session),
178  m_graph(in.m_graph),
179  m_ncoeffs(0),
180  m_npoints(0),
181  m_coeffs(),
182  m_phys(),
183  m_physState(false),
184  m_exp(MemoryManager<LocalRegions::ExpansionVector>
185  ::AllocateSharedPtr()),
186  m_coeff_offset(),
187  m_phys_offset(),
188  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
189  m_WaveSpace(false)
190  {
192 
193  for (int i=0; i < eIDs.size(); ++i)
194  {
195  (*m_exp).push_back( (*(in.m_exp))[eIDs[i]]);
196  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
197  m_npoints += (*m_exp)[i]->GetTotPoints();
198  }
199 
200  if(DeclareCoeffPhysArrays)
201  {
204  }
205  }
206 
207 
208  /**
209  * Copies an existing expansion list.
210  * @param in Source expansion list.
211  */
212  ExpList::ExpList(const ExpList &in, const bool DeclareCoeffPhysArrays):
213  std::enable_shared_from_this<ExpList>(in),
214  m_comm(in.m_comm),
215  m_session(in.m_session),
216  m_graph(in.m_graph),
217  m_ncoeffs(in.m_ncoeffs),
218  m_npoints(in.m_npoints),
219  m_physState(false),
220  m_exp(in.m_exp),
228  m_WaveSpace(false)
229  {
231 
232  if(DeclareCoeffPhysArrays)
233  {
236  }
237  }
238 
239 
240  /**
241  * Each expansion (local element) is processed in turn to
242  * determine the number of coefficients and physical data
243  * points it contributes to the domain. Twoe arrays,
244  * #m_coeff_offset are #m_phys_offset are also initialised and
245  * updated to store the data offsets of each element in the
246  * #m_coeffs and #m_phys arrays, and the element id that each
247  * consecutive block is associated respectively.
248  */
250  {
251  int i;
252 
253  // Set up offset information and array sizes
256 
257  m_ncoeffs = m_npoints = 0;
258 
259  for(i = 0; i < m_exp->size(); ++i)
260  {
262  m_phys_offset [i] = m_npoints;
263  m_ncoeffs += (*m_exp)[i]->GetNcoeffs();
264  m_npoints += (*m_exp)[i]->GetTotPoints();
265  }
266  }
267 
268  /**
269  *
270  */
272  {
273  return m_expType;
274  }
275 
276  /**
277  *
278  */
280  {
281  m_expType = Type;
282  }
283 
285  {
286  }
287 
288  /**
289  * The integration is evaluated locally, that is
290  * \f[\int
291  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
292  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
293  * where the integration over the separate elements is done by the
294  * function StdRegions#StdExpansion#Integral, which discretely
295  * evaluates the integral using Gaussian quadrature.
296  *
297  * Note that the array #m_phys should be filled with the values of the
298  * function \f$f(\boldsymbol{x})\f$ at the quadrature points
299  * \f$\boldsymbol{x}_i\f$.
300  *
301  * @return The value of the discretely evaluated integral
302  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
303  */
305  {
306  ASSERTL1(m_physState == true,
307  "local physical space is not true ");
308 
309  return PhysIntegral(m_phys);
310  }
311 
312 
313  /**
314  * The integration is evaluated locally, that is
315  * \f[\int
316  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
317  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
318  * where the integration over the separate elements is done by the
319  * function StdRegions#StdExpansion#Integral, which discretely
320  * evaluates the integral using Gaussian quadrature.
321  *
322  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
323  * containing the values of the function
324  * \f$f(\boldsymbol{x})\f$ at the quadrature
325  * points \f$\boldsymbol{x}_i\f$.
326  * @return The value of the discretely evaluated integral
327  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
328  */
330  const Array<OneD, const NekDouble> &inarray)
331  {
332  int i;
333  NekDouble sum = 0.0;
334 
335  for(i = 0; i < (*m_exp).size(); ++i)
336  {
337  sum += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
338  }
339 
340  return sum;
341  }
342 
343 
344  /**
345  * Retrieves the block matrix specified by \a bkey, and computes
346  * \f$ y=Mx \f$.
347  * @param gkey GlobalMatrixKey specifying the block matrix to
348  * use in the matrix-vector multiply.
349  * @param inarray Input vector \f$ x \f$.
350  * @param outarray Output vector \f$ y \f$.
351  */
353  const GlobalMatrixKey &gkey,
354  const Array<OneD,const NekDouble> &inarray,
355  Array<OneD, NekDouble> &outarray)
356  {
357  // Retrieve the block matrix using the given key.
358  const DNekScalBlkMatSharedPtr& blockmat = GetBlockMatrix(gkey);
359  int nrows = blockmat->GetRows();
360  int ncols = blockmat->GetColumns();
361 
362  // Create NekVectors from the given data arrays
363  NekVector<NekDouble> in (ncols,inarray, eWrapper);
364  NekVector< NekDouble> out(nrows,outarray,eWrapper);
365 
366  // Perform matrix-vector multiply.
367  out = (*blockmat)*in;
368  }
369 
370 
371  /**
372  * The operation is evaluated locally for every element by the function
373  * StdRegions#StdExpansion#IProductWRTBase.
374  *
375  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
376  * containing the values of the function
377  * \f$f(\boldsymbol{x})\f$ at the quadrature
378  * points \f$\boldsymbol{x}_i\f$.
379  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
380  * used to store the result.
381  */
383  const Array<OneD, const NekDouble> &inarray,
384  Array<OneD, NekDouble> &outarray)
385  {
387  for (int i = 0; i < m_collections.size(); ++i)
388  {
389 
391  inarray + m_coll_phys_offset[i],
392  tmp = outarray + m_coll_coeff_offset[i]);
393  }
394  }
395 
396  /**
397  * The operation is evaluated locally for every element by the function
398  * StdRegions#StdExpansion#IProductWRTDerivBase.
399  *
400  * @param dir {0,1} is the direction in which the
401  * derivative of the basis should be taken
402  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
403  * containing the values of the function
404  * \f$f(\boldsymbol{x})\f$ at the quadrature
405  * points \f$\boldsymbol{x}_i\f$.
406  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
407  * used to store the result.
408  */
409  void ExpList::IProductWRTDerivBase(const int dir,
410  const Array<OneD, const NekDouble> &inarray,
411  Array<OneD, NekDouble> &outarray)
412  {
413  int i;
414 
415  Array<OneD,NekDouble> e_outarray;
416 
417  for(i = 0; i < (*m_exp).size(); ++i)
418  {
419  (*m_exp)[i]->IProductWRTDerivBase(dir,inarray+m_phys_offset[i],
420  e_outarray = outarray+m_coeff_offset[i]);
421  }
422  }
423 
424 
425  /**
426  * @brief Directional derivative along a given direction
427  *
428  */
430  const Array<OneD, const NekDouble> &direction,
431  const Array<OneD, const NekDouble> &inarray,
432  Array<OneD, NekDouble> &outarray)
433  {
434  int npts_e;
435  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
436  int nq = direction.num_elements()/coordim;
437 
438  Array<OneD, NekDouble> e_outarray;
439  Array<OneD, NekDouble> e_MFdiv;
440 
441  Array<OneD, NekDouble> locdir;
442 
443  for(int i = 0; i < (*m_exp).size(); ++i)
444  {
445  npts_e = (*m_exp)[i]->GetTotPoints();
446  locdir = Array<OneD, NekDouble>(npts_e*coordim);
447 
448  for (int k = 0; k<coordim; ++k)
449  {
450  Vmath::Vcopy(npts_e, &direction[k*nq+m_phys_offset[i]], 1,
451  &locdir[k*npts_e], 1);
452  }
453 
454  (*m_exp)[i]->IProductWRTDirectionalDerivBase(
455  locdir,
456  inarray+m_phys_offset[i],
457  e_outarray = outarray+m_coeff_offset[i] );
458  }
459  }
460 
461 
462  /**
463  * The operation is evaluated locally for every element by the function
464  * StdRegions#StdExpansion#IProductWRTDerivBase.
465  *
466  * @param inarray An array of arrays of size \f$Q_{\mathrm{tot}}\f$
467  * containing the values of the function
468  * \f$f(\boldsymbol{x})\f$ at the quadrature
469  * points \f$\boldsymbol{x}_i\f$ in dir directions.
470  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
471  * used to store the result.
472  */
474  Array<OneD, NekDouble> &outarray)
475  {
476  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
477  // assume coord dimension defines the size of Deriv Base
478  int dim = GetCoordim(0);
479 
480  ASSERTL1(inarray.num_elements() >= dim,"inarray is not of sufficient dimension");
481 
482  switch(dim)
483  {
484  case 1:
485  for (int i = 0; i < m_collections.size(); ++i)
486  {
487  m_collections[i].ApplyOperator(
489  inarray[0] + m_coll_phys_offset[i],
490  tmp0 = outarray + m_coll_coeff_offset[i]);
491  }
492  break;
493  case 2:
494  for (int i = 0; i < m_collections.size(); ++i)
495  {
496  m_collections[i].ApplyOperator(
498  inarray[0] + m_coll_phys_offset[i],
499  tmp0 = inarray[1] + m_coll_phys_offset[i],
500  tmp1 = outarray + m_coll_coeff_offset[i]);
501  }
502  break;
503  case 3:
504  for (int i = 0; i < m_collections.size(); ++i)
505  {
506  m_collections[i].ApplyOperator(
508  inarray[0] + m_coll_phys_offset[i],
509  tmp0 = inarray[1] + m_coll_phys_offset[i],
510  tmp1 = inarray[2] + m_coll_phys_offset[i],
511  tmp2 = outarray + m_coll_coeff_offset[i]);
512  }
513  break;
514  default:
515  ASSERTL0(false,"Dimension of inarray not correct");
516  break;
517  }
518  }
519  /**
520  * Given a function \f$f(\boldsymbol{x})\f$ evaluated at
521  * the quadrature points, this function calculates the
522  * derivatives \f$\frac{d}{dx_1}\f$, \f$\frac{d}{dx_2}\f$
523  * and \f$\frac{d}{dx_3}\f$ of the function
524  * \f$f(\boldsymbol{x})\f$ at the same quadrature
525  * points. The local distribution of the quadrature points
526  * allows an elemental evaluation of the derivative. This
527  * is done by a call to the function
528  * StdRegions#StdExpansion#PhysDeriv.
529  *
530  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
531  * containing the values of the function
532  * \f$f(\boldsymbol{x})\f$ at the quadrature
533  * points \f$\boldsymbol{x}_i\f$.
534  * @param out_d0 The discrete evaluation of the
535  * derivative\f$\frac{d}{dx_1}\f$ will
536  * be stored in this array of size
537  * \f$Q_{\mathrm{tot}}\f$.
538  * @param out_d1 The discrete evaluation of the
539  * derivative\f$\frac{d}{dx_2}\f$ will be
540  * stored in this array of size
541  * \f$Q_{\mathrm{tot}}\f$. Note that if no
542  * memory is allocated for \a out_d1, the
543  * derivative \f$\frac{d}{dx_2}\f$ will not be
544  * calculated.
545  * @param out_d2 The discrete evaluation of the
546  * derivative\f$\frac{d}{dx_3}\f$ will be
547  * stored in this array of size
548  * \f$Q_{\mathrm{tot}}\f$. Note that if no
549  * memory is allocated for \a out_d2, the
550  * derivative \f$\frac{d}{dx_3}\f$ will not be
551  * calculated.
552  */
554  Array<OneD, NekDouble> &out_d0,
555  Array<OneD, NekDouble> &out_d1,
556  Array<OneD, NekDouble> &out_d2)
557  {
558  Array<OneD, NekDouble> e_out_d0;
559  Array<OneD, NekDouble> e_out_d1;
560  Array<OneD, NekDouble> e_out_d2;
561  int offset;
562  for (int i = 0; i < m_collections.size(); ++i)
563  {
564  offset = m_coll_phys_offset[i];
565  e_out_d0 = out_d0 + offset;
566  e_out_d1 = out_d1 + offset;
567  e_out_d2 = out_d2 + offset;
568 
569  m_collections[i].ApplyOperator(Collections::ePhysDeriv,
570  inarray + offset,
571  e_out_d0,e_out_d1, e_out_d2);
572 
573  }
574  }
575 
576  void ExpList::v_PhysDeriv(const int dir,
577  const Array<OneD, const NekDouble> &inarray,
578  Array<OneD, NekDouble> &out_d)
579  {
580  Direction edir = DirCartesianMap[dir];
581  v_PhysDeriv(edir, inarray,out_d);
582  }
583 
585  Array<OneD, NekDouble> &out_d)
586  {
587  int i;
588  if(edir==MultiRegions::eS)
589  {
590  Array<OneD, NekDouble> e_out_ds;
591  for(i=0; i<(*m_exp).size(); ++i)
592  {
593  e_out_ds = out_d + m_phys_offset[i];
594  (*m_exp)[i]->PhysDeriv_s(inarray+m_phys_offset[i],e_out_ds);
595  }
596  }
597  else if(edir==MultiRegions::eN)
598  {
599  Array<OneD, NekDouble > e_out_dn;
600  for(i=0; i<(*m_exp).size(); i++)
601  {
602  e_out_dn = out_d +m_phys_offset[i];
603  (*m_exp)[i]->PhysDeriv_n(inarray+m_phys_offset[i],e_out_dn);
604  }
605  }
606  else
607  {
608  // convert enum into int
609  int intdir= (int)edir;
610  Array<OneD, NekDouble> e_out_d;
611  int offset;
612  for (int i = 0; i < m_collections.size(); ++i)
613  {
614  offset = m_coll_phys_offset[i];
615  e_out_d = out_d + offset;
616 
617  m_collections[i].ApplyOperator(Collections::ePhysDeriv,
618  intdir,
619  inarray + offset,
620  e_out_d);
621  }
622  }
623  }
624 
628  {
629  int nq = GetTotPoints();
630  Array<OneD,NekDouble> Vx(nq);
631  Array<OneD,NekDouble> Uy(nq);
633 
634  bool halfMode = false;
635  if ( GetExpType() == e3DH1D)
636  {
637  m_session->MatchSolverInfo("ModeType", "HalfMode",
638  halfMode, false);
639  }
640 
641  switch(GetExpType())
642  {
643  case e2D:
644  {
645  PhysDeriv(xDir, Vel[yDir], Vx);
646  PhysDeriv(yDir, Vel[xDir], Uy);
647 
648 
649  Vmath::Vsub(nq, Vx, 1, Uy, 1, Dummy, 1);
650 
651  PhysDeriv(Dummy,Q[1],Q[0]);
652 
653  Vmath::Smul(nq, -1.0, Q[1], 1, Q[1], 1);
654  }
655  break;
656 
657  case e3D:
658  case e3DH1D:
659  case e3DH2D:
660  {
661  Array<OneD,NekDouble> Vz(nq);
662  Array<OneD,NekDouble> Uz(nq);
663  Array<OneD,NekDouble> Wx(nq);
664  Array<OneD,NekDouble> Wy(nq);
665 
666  PhysDeriv(Vel[xDir], Dummy, Uy, Uz);
667  PhysDeriv(Vel[yDir], Vx, Dummy, Vz);
668  PhysDeriv(Vel[zDir], Wx, Wy, Dummy);
669 
670  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
671  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
672  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
673 
674  PhysDeriv(Q[0], Dummy, Uy, Uz);
675  PhysDeriv(Q[1], Vx, Dummy, Vz);
676  PhysDeriv(Q[2], Wx, Wy, Dummy);
677 
678  // For halfmode, need to change the sign of z derivatives
679  if (halfMode)
680  {
681  Vmath::Neg(nq, Uz, 1);
682  Vmath::Neg(nq, Vz, 1);
683  }
684 
685  Vmath::Vsub(nq, Wy, 1, Vz, 1, Q[0], 1);
686  Vmath::Vsub(nq, Uz, 1, Wx, 1, Q[1], 1);
687  Vmath::Vsub(nq, Vx, 1, Uy, 1, Q[2], 1);
688  }
689  break;
690  default:
691  ASSERTL0(0,"Dimension not supported");
692  break;
693  }
694  }
695 
696 
698  const Array<OneD, const NekDouble> &direction,
699  const Array<OneD, const NekDouble> &inarray,
700  Array<OneD, NekDouble> &outarray)
701  {
702  int npts_e;
703  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
704  int nq = direction.num_elements() / coordim;
705 
706  Array<OneD, NekDouble> e_outarray;
707  Array<OneD, NekDouble> e_MFdiv;
708  Array<OneD, NekDouble> locdir;
709 
710  for(int i = 0; i < (*m_exp).size(); ++i)
711  {
712  npts_e = (*m_exp)[i]->GetTotPoints();
713  locdir = Array<OneD, NekDouble>(npts_e*coordim);
714 
715  for (int k = 0; k<coordim; ++k)
716  {
717  Vmath::Vcopy(npts_e, &direction[k*nq+m_phys_offset[i]], 1,
718  &locdir[k*npts_e], 1);
719  }
720 
721  (*m_exp)[i]->PhysDirectionalDeriv(
722  inarray + m_phys_offset[i],
723  locdir,
724  e_outarray = outarray+m_phys_offset[i] );
725  }
726  }
727 
728 
730  Array<OneD, NekDouble> &array,
731  const NekDouble alpha,
732  const NekDouble exponent,
733  const NekDouble cutoff)
734  {
735  Array<OneD,NekDouble> e_array;
736 
737  for(int i = 0; i < (*m_exp).size(); ++i)
738  {
739  (*m_exp)[i]->ExponentialFilter(
740  e_array = array+m_phys_offset[i],
741  alpha,
742  exponent,
743  cutoff);
744  }
745  }
746 
747  /**
748  * The coefficients of the function to be acted upon
749  * should be contained in the \param inarray. The
750  * resulting coefficients are stored in \param outarray
751  *
752  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
753  * containing the inner product.
754  */
756  const Array<OneD, const NekDouble> &inarray,
757  Array<OneD, NekDouble> &outarray)
758  {
760  const DNekScalBlkMatSharedPtr& InvMass = GetBlockMatrix(mkey);
761 
762  // Inverse mass matrix
764  if(inarray.get() == outarray.get())
765  {
766  NekVector<NekDouble> in(m_ncoeffs,inarray); // copy data
767  out = (*InvMass)*in;
768  }
769  else
770  {
772  out = (*InvMass)*in;
773  }
774  }
775 
776  /**
777  * Given a function \f$u(\boldsymbol{x})\f$ defined at the
778  * quadrature points, this function determines the
779  * transformed elemental coefficients \f$\hat{u}_n^e\f$
780  * employing a discrete elemental Galerkin projection from
781  * physical space to coefficient space. For each element,
782  * the operation is evaluated locally by the function
783  * StdRegions#StdExpansion#IproductWRTBase followed by a
784  * call to #MultiRegions#MultiplyByElmtInvMass.
785  *
786  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
787  * containing the values of the function
788  * \f$f(\boldsymbol{x})\f$ at the quadrature
789  * points \f$\boldsymbol{x}_i\f$.
790  * @param outarray The resulting coefficients
791  * \f$\hat{u}_n^e\f$ will be stored in this
792  * array of size \f$N_{\mathrm{eof}}\f$.
793  */
795  Array<OneD, NekDouble> &outarray)
796  {
798 
799  IProductWRTBase_IterPerExp(inarray,f);
800  MultiplyByElmtInvMass(f,outarray);
801 
802  }
803 
805  const Array<OneD, const NekDouble>& inarray,
806  Array<OneD, NekDouble> &outarray)
807  {
808  int i;
809 
810  Array<OneD,NekDouble> e_outarray;
811 
812  for(i= 0; i < (*m_exp).size(); ++i)
813  {
814  (*m_exp)[i]->FwdTrans_BndConstrained(inarray+m_phys_offset[i],
815  e_outarray = outarray+m_coeff_offset[i]);
816  }
817  }
818 
819  /**
820  * This function smooth a field after some calculaitons which have
821  * been done elementally.
822  *
823  * @param field An array containing the field in physical space
824  *
825  */
827  {
828  boost::ignore_unused(field);
829  // Do nothing unless the method is implemented in the appropriate
830  // class, i.e. ContField1D,ContField2D, etc.
831 
832  // So far it has been implemented just for ContField2D and
833  // ContField3DHomogeneous1D
834 
835  // Block in case users try the smoothing with a modal expansion.
836  // Maybe a different techique for the smoothing require
837  // implementation for modal basis.
838 
839  ASSERTL0((*m_exp)[0]->GetBasisType(0)
841  (*m_exp)[0]->GetBasisType(0)
843  "Smoothing is currently not allowed unless you are using "
844  "a nodal base for efficiency reasons. The implemented "
845  "smoothing technique requires the mass matrix inversion "
846  "which is trivial just for GLL_LAGRANGE_SEM and "
847  "GAUSS_LAGRANGE_SEMexpansions.");
848  }
849 
850 
851  /**
852  * This function assembles the block diagonal matrix
853  * \f$\underline{\boldsymbol{M}}^e\f$, which is the
854  * concatenation of the local matrices
855  * \f$\boldsymbol{M}^e\f$ of the type \a mtype, that is
856  *
857  * \f[
858  * \underline{\boldsymbol{M}}^e = \left[
859  * \begin{array}{cccc}
860  * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
861  * 0 & \boldsymbol{M}^2 & 0 & 0 \\
862  * 0 & 0 & \ddots & 0 \\
863  * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array}\right].\f]
864  *
865  * @param mtype the type of matrix to be assembled
866  * @param scalar an optional parameter
867  * @param constant an optional parameter
868  */
870  const GlobalMatrixKey &gkey)
871  {
872  int i,cnt1;
873  int n_exp = 0;
874  DNekScalMatSharedPtr loc_mat;
875  DNekScalBlkMatSharedPtr BlkMatrix;
876  map<int,int> elmt_id;
878 
879  if(ShapeType != LibUtilities::eNoShapeType)
880  {
881  for(i = 0 ; i < (*m_exp).size(); ++i)
882  {
883  if((*m_exp)[i]->DetShapeType()
884  == ShapeType)
885  {
886  elmt_id[n_exp++] = i;
887  }
888  }
889  }
890  else
891  {
892  n_exp = (*m_exp).size();
893  for(i = 0; i < n_exp; ++i)
894  {
895  elmt_id[i] = i;
896  }
897  }
898 
899  Array<OneD,unsigned int> nrows(n_exp);
900  Array<OneD,unsigned int> ncols(n_exp);
901 
902  switch(gkey.GetMatrixType())
903  {
905  {
906  // set up an array of integers for block matrix construction
907  for(i = 0; i < n_exp; ++i)
908  {
909  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
910  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
911  }
912  }
913  break;
915  {
916  // set up an array of integers for block matrix construction
917  for(i = 0; i < n_exp; ++i)
918  {
919  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
920  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
921  }
922  }
923  break;
924  case StdRegions::eMass:
929  {
930  // set up an array of integers for block matrix construction
931  for(i = 0; i < n_exp; ++i)
932  {
933  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
934  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
935  }
936  }
937  break;
938 
940  {
941  // set up an array of integers for block matrix construction
942  for(i = 0; i < n_exp; ++i)
943  {
944  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
945  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
946  }
947  }
948  break;
949 
951  {
952  // set up an array of integers for block matrix construction
953  for(i = 0; i < n_exp; ++i)
954  {
955  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
956  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
957  }
958  }
959  break;
960 
961  default:
962  {
964  "Global Matrix creation not defined for this type "
965  "of matrix");
966  }
967  }
968 
969  MatrixStorage blkmatStorage = eDIAGONAL;
971  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
972 
973  int nvarcoeffs = gkey.GetNVarCoeffs();
974  int eid;
975  Array<OneD, NekDouble> varcoeffs_wk;
976 
977  for(i = cnt1 = 0; i < n_exp; ++i)
978  {
979  // need to be initialised with zero size for non variable coefficient case
980  StdRegions::VarCoeffMap varcoeffs;
981 
982  eid = elmt_id[i];
983  if(nvarcoeffs>0)
984  {
985  for (auto &x : gkey.GetVarCoeffs())
986  {
987  varcoeffs[x.first] = x.second + m_phys_offset[eid];
988  }
989  }
990 
992  (*m_exp)[eid]->DetShapeType(),
993  *(*m_exp)[eid],
994  gkey.GetConstFactors(),
995  varcoeffs );
996 
997  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[elmt_id.find(i)->second])->GetLocMatrix(matkey);
998  BlkMatrix->SetBlock(i,i,loc_mat);
999  }
1000 
1001  return BlkMatrix;
1002  }
1003 
1005  const GlobalMatrixKey &gkey)
1006  {
1007  auto matrixIter = m_blockMat->find(gkey);
1008 
1009  if(matrixIter == m_blockMat->end())
1010  {
1011  return ((*m_blockMat)[gkey] = GenBlockMatrix(gkey));
1012  }
1013  else
1014  {
1015  return matrixIter->second;
1016  }
1017  }
1018 
1020  const GlobalMatrixKey &gkey,
1021  const Array<OneD,const NekDouble> &inarray,
1022  Array<OneD, NekDouble> &outarray)
1023  {
1024  const Array<OneD, const bool> doBlockMatOp
1025  = m_globalOptParam->DoBlockMatOp(gkey.GetMatrixType());
1026  const Array<OneD, const int> num_elmts
1027  = m_globalOptParam->GetShapeNumElements();
1028 
1029  Array<OneD,NekDouble> tmp_outarray;
1030  int cnt = 0;
1031  int eid;
1032  for(int n = 0; n < num_elmts.num_elements(); ++n)
1033  {
1034  if(doBlockMatOp[n])
1035  {
1036  const LibUtilities::ShapeType vType
1037  = m_globalOptParam->GetShapeList()[n];
1038  const MultiRegions::GlobalMatrixKey vKey(gkey, vType);
1039  if (cnt < m_coeff_offset.num_elements())
1040  {
1041  eid = cnt;
1042  MultiplyByBlockMatrix(vKey,inarray + m_coeff_offset[eid],
1043  tmp_outarray = outarray + m_coeff_offset[eid]);
1044  cnt += num_elmts[n];
1045  }
1046  }
1047  else
1048  {
1049  int i;
1050  int nvarcoeffs = gkey.GetNVarCoeffs();
1051 
1052  for(i= 0; i < num_elmts[n]; ++i)
1053  {
1054  // need to be initialised with zero size for non variable coefficient case
1055  StdRegions::VarCoeffMap varcoeffs;
1056 
1057  eid = cnt++;
1058  if(nvarcoeffs>0)
1059  {
1060  for (auto &x : gkey.GetVarCoeffs())
1061  {
1062  varcoeffs[x.first] = x.second + m_phys_offset[eid];
1063  }
1064  }
1065 
1067  (*m_exp)[eid]->DetShapeType(),
1068  *((*m_exp)[eid]),
1069  gkey.GetConstFactors(),varcoeffs);
1070 
1071  (*m_exp)[eid]->GeneralMatrixOp(inarray + m_coeff_offset[eid],
1072  tmp_outarray = outarray+m_coeff_offset[eid],
1073  mkey);
1074  }
1075  }
1076  }
1077  }
1078 
1079  /**
1080  * Retrieves local matrices from each expansion in the expansion list
1081  * and combines them together to generate a global matrix system.
1082  * @param mkey Matrix key for the matrix to be generated.
1083  * @param locToGloMap Local to global mapping.
1084  * @returns Shared pointer to the generated global matrix.
1085  */
1087  const GlobalMatrixKey &mkey,
1088  const AssemblyMapCGSharedPtr &locToGloMap)
1089  {
1090  int i,j,n,gid1,gid2,cntdim1,cntdim2;
1091  NekDouble sign1,sign2;
1092  DNekScalMatSharedPtr loc_mat;
1093 
1094  unsigned int glob_rows = 0;
1095  unsigned int glob_cols = 0;
1096  unsigned int loc_rows = 0;
1097  unsigned int loc_cols = 0;
1098 
1099  bool assembleFirstDim = false;
1100  bool assembleSecondDim = false;
1101 
1102  switch(mkey.GetMatrixType())
1103  {
1104  case StdRegions::eBwdTrans:
1105  {
1106  glob_rows = m_npoints;
1107  glob_cols = locToGloMap->GetNumGlobalCoeffs();
1108 
1109  assembleFirstDim = false;
1110  assembleSecondDim = true;
1111  }
1112  break;
1114  {
1115  glob_rows = locToGloMap->GetNumGlobalCoeffs();
1116  glob_cols = m_npoints;
1117 
1118  assembleFirstDim = true;
1119  assembleSecondDim = false;
1120  }
1121  break;
1122  case StdRegions::eMass:
1126  {
1127  glob_rows = locToGloMap->GetNumGlobalCoeffs();
1128  glob_cols = locToGloMap->GetNumGlobalCoeffs();
1129 
1130  assembleFirstDim = true;
1131  assembleSecondDim = true;
1132  }
1133  break;
1134  default:
1135  {
1137  "Global Matrix creation not defined for this type "
1138  "of matrix");
1139  }
1140  }
1141 
1142  COOMatType spcoomat;
1143  CoordType coord;
1144 
1145  int nvarcoeffs = mkey.GetNVarCoeffs();
1146  int eid;
1147 
1148  // fill global matrix
1149  for(n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
1150  {
1151  // need to be initialised with zero size for non variable coefficient case
1152  StdRegions::VarCoeffMap varcoeffs;
1153 
1154  eid = n;
1155  if(nvarcoeffs>0)
1156  {
1157  for (auto &x : mkey.GetVarCoeffs())
1158  {
1159  varcoeffs[x.first] = x.second + m_phys_offset[eid];
1160  }
1161  }
1162 
1163  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
1164  (*m_exp)[eid]->DetShapeType(),
1165  *((*m_exp)[eid]),
1166  mkey.GetConstFactors(),varcoeffs);
1167 
1168  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])->GetLocMatrix(matkey);
1169 
1170  loc_rows = loc_mat->GetRows();
1171  loc_cols = loc_mat->GetColumns();
1172 
1173  for(i = 0; i < loc_rows; ++i)
1174  {
1175  if(assembleFirstDim)
1176  {
1177  gid1 = locToGloMap->GetLocalToGlobalMap (cntdim1 + i);
1178  sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
1179  }
1180  else
1181  {
1182  gid1 = cntdim1 + i;
1183  sign1 = 1.0;
1184  }
1185 
1186  for(j = 0; j < loc_cols; ++j)
1187  {
1188  if(assembleSecondDim)
1189  {
1190  gid2 = locToGloMap
1191  ->GetLocalToGlobalMap(cntdim2 + j);
1192  sign2 = locToGloMap
1193  ->GetLocalToGlobalSign(cntdim2 + j);
1194  }
1195  else
1196  {
1197  gid2 = cntdim2 + j;
1198  sign2 = 1.0;
1199  }
1200 
1201  // sparse matrix fill
1202  coord = make_pair(gid1,gid2);
1203  if( spcoomat.count(coord) == 0 )
1204  {
1205  spcoomat[coord] = sign1*sign2*(*loc_mat)(i,j);
1206  }
1207  else
1208  {
1209  spcoomat[coord] += sign1*sign2*(*loc_mat)(i,j);
1210  }
1211  }
1212  }
1213  cntdim1 += loc_rows;
1214  cntdim2 += loc_cols;
1215  }
1216 
1218  ::AllocateSharedPtr(m_session,glob_rows,glob_cols,spcoomat);
1219  }
1220 
1221 
1223  {
1224  int i,j,n,gid1,gid2,loc_lda,eid;
1225  NekDouble sign1,sign2,value;
1226  DNekScalMatSharedPtr loc_mat;
1227 
1228  int totDofs = locToGloMap->GetNumGlobalCoeffs();
1229  int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
1230 
1231  unsigned int rows = totDofs - NumDirBCs;
1232  unsigned int cols = totDofs - NumDirBCs;
1233  NekDouble zero = 0.0;
1234 
1235  DNekMatSharedPtr Gmat;
1236  int bwidth = locToGloMap->GetFullSystemBandWidth();
1237 
1238  int nvarcoeffs = mkey.GetNVarCoeffs();
1239  MatrixStorage matStorage;
1240 
1241  map<int, RobinBCInfoSharedPtr> RobinBCInfo = GetRobinBCInfo();
1242 
1243  switch(mkey.GetMatrixType())
1244  {
1245  // case for all symmetric matices
1248  if( (2*(bwidth+1)) < rows)
1249  {
1251  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage,bwidth,bwidth);
1252  }
1253  else
1254  {
1255  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
1256  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
1257  }
1258 
1259  break;
1260  default: // Assume general matrix - currently only set up for full invert
1261  {
1262  matStorage = eFULL;
1263  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
1264  }
1265  }
1266 
1267  // fill global symmetric matrix
1268  for(n = 0; n < (*m_exp).size(); ++n)
1269  {
1270  // need to be initialised with zero size for non variable coefficient case
1271  StdRegions::VarCoeffMap varcoeffs;
1272 
1273  eid = n;
1274  if(nvarcoeffs>0)
1275  {
1276  for (auto &x : mkey.GetVarCoeffs())
1277  {
1278  varcoeffs[x.first] = x.second + m_phys_offset[eid];
1279  }
1280  }
1281 
1282  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
1283  (*m_exp)[eid]->DetShapeType(),
1284  *((*m_exp)[eid]),
1285  mkey.GetConstFactors(),varcoeffs);
1286 
1287  loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])->GetLocMatrix(matkey);
1288 
1289 
1290  if(RobinBCInfo.count(n) != 0) // add robin mass matrix
1291  {
1293 
1294  // declare local matrix from scaled matrix.
1295  int rows = loc_mat->GetRows();
1296  int cols = loc_mat->GetColumns();
1297  const NekDouble *dat = loc_mat->GetRawPtr();
1299  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
1300 
1301  // add local matrix contribution
1302  for(rBC = RobinBCInfo.find(n)->second;rBC; rBC = rBC->next)
1303  {
1304  (*m_exp)[n]->AddRobinMassMatrix(rBC->m_robinID,rBC->m_robinPrimitiveCoeffs,new_mat);
1305  }
1306 
1307  NekDouble one = 1.0;
1308  // redeclare loc_mat to point to new_mat plus the scalar.
1309  loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,new_mat);
1310  }
1311 
1312  loc_lda = loc_mat->GetColumns();
1313 
1314  for(i = 0; i < loc_lda; ++i)
1315  {
1316  gid1 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + i) - NumDirBCs;
1317  sign1 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + i);
1318  if(gid1 >= 0)
1319  {
1320  for(j = 0; j < loc_lda; ++j)
1321  {
1322  gid2 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + j) - NumDirBCs;
1323  sign2 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + j);
1324  if(gid2 >= 0)
1325  {
1326  // When global matrix is symmetric,
1327  // only add the value for the upper
1328  // triangular part in order to avoid
1329  // entries to be entered twice
1330  if((matStorage == eFULL)||(gid2 >= gid1))
1331  {
1332  value = Gmat->GetValue(gid1,gid2) + sign1*sign2*(*loc_mat)(i,j);
1333  Gmat->SetValue(gid1,gid2,value);
1334  }
1335  }
1336  }
1337  }
1338  }
1339  }
1340 
1341  return Gmat;
1342  }
1343 
1344 
1345  /**
1346  * Consider a linear system
1347  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{f}\f$ to be solved. Dependent
1348  * on the solution method, this function constructs
1349  * - <b>The full linear system</b><BR>
1350  * A call to the function #GenGlobalLinSysFullDirect
1351  * - <b>The statically condensed linear system</b><BR>
1352  * A call to the function #GenGlobalLinSysStaticCond
1353  *
1354  * @param mkey A key which uniquely defines the global
1355  * matrix to be constructed.
1356  * @param locToGloMap Contains the mapping array and required
1357  * information for the transformation from
1358  * local to global degrees of freedom.
1359  * @return (A shared pointer to) the global linear system in
1360  * required format.
1361  */
1363  const GlobalLinSysKey &mkey,
1364  const AssemblyMapCGSharedPtr &locToGloMap)
1365  {
1366  GlobalLinSysSharedPtr returnlinsys;
1367  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
1368 
1370 
1371  if (vType >= eSIZE_GlobalSysSolnType)
1372  {
1373  ASSERTL0(false,"Matrix solution type not defined");
1374  }
1375  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
1376 
1377  return GetGlobalLinSysFactory().CreateInstance( vSolnType, mkey,
1378  vExpList, locToGloMap);
1379  }
1380 
1382  const GlobalLinSysKey &mkey,
1383  const AssemblyMapSharedPtr &locToGloMap)
1384  {
1385  std::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
1386  const map<int,RobinBCInfoSharedPtr> vRobinBCInfo = GetRobinBCInfo();
1387 
1389 
1390  if (vType >= eSIZE_GlobalSysSolnType)
1391  {
1392  ASSERTL0(false,"Matrix solution type not defined");
1393  }
1394  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
1395 
1396  return GetGlobalLinSysFactory().CreateInstance(vSolnType,mkey,
1397  vExpList,locToGloMap);
1398  }
1399 
1400 
1401  /**
1402  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of
1403  * an expansion, this function evaluates the spectral/hp
1404  * expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
1405  * quadrature points \f$\boldsymbol{x}_i\f$. The operation
1406  * is evaluated locally by the elemental function
1407  * StdRegions#StdExpansion#BwdTrans.
1408  *
1409  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
1410  * containing the local coefficients
1411  * \f$\hat{u}_n^e\f$.
1412  * @param outarray The resulting physical values at the
1413  * quadrature points
1414  * \f$u^{\delta}(\boldsymbol{x}_i)\f$
1415  * will be stored in this array of size
1416  * \f$Q_{\mathrm{tot}}\f$.
1417  */
1419  Array<OneD, NekDouble> &outarray)
1420  {
1422  for (int i = 0; i < m_collections.size(); ++i)
1423  {
1424  m_collections[i].ApplyOperator(Collections::eBwdTrans,
1425  inarray + m_coll_coeff_offset[i],
1426  tmp = outarray + m_coll_phys_offset[i]);
1427  }
1428  }
1429 
1431  const Array<OneD, const NekDouble> &gloCoord)
1432  {
1433  Array<OneD, NekDouble> stdCoord(GetCoordim(0),0.0);
1434  for (int i = 0; i < (*m_exp).size(); ++i)
1435  {
1436  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoord))
1437  {
1438  return (*m_exp)[i];
1439  }
1440  }
1441  ASSERTL0(false, "Cannot find element for this point.");
1442  return (*m_exp)[0]; // avoid warnings
1443  }
1444 
1445 
1446  /**
1447  * @todo need a smarter search here that first just looks at bounding
1448  * vertices - suggest first seeing if point is within 10% of
1449  * region defined by vertices. The do point search.
1450  */
1452  const Array<OneD, const NekDouble> &gloCoord,
1453  NekDouble tol,
1454  bool returnNearestElmt)
1455  {
1456  Array<OneD, NekDouble> Lcoords(gloCoord.num_elements());
1457 
1458  return GetExpIndex(gloCoord,Lcoords,tol,returnNearestElmt);
1459  }
1460 
1461 
1463  const Array<OneD, const NekDouble> &gloCoords,
1464  Array<OneD, NekDouble> &locCoords,
1465  NekDouble tol,
1466  bool returnNearestElmt)
1467  {
1468  if (GetNumElmts() == 0)
1469  {
1470  return -1;
1471  }
1472 
1473  if (m_elmtToExpId.size() == 0)
1474  {
1475  // Loop in reverse order so that in case where using a
1476  // Homogeneous expansion it sets geometry ids to first part of
1477  // m_exp list. Otherwise will set to second (complex) expansion
1478  for(int i = (*m_exp).size()-1; i >= 0; --i)
1479  {
1480  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1481  }
1482  }
1483 
1484  NekDouble x = (gloCoords.num_elements() > 0 ? gloCoords[0] : 0.0);
1485  NekDouble y = (gloCoords.num_elements() > 1 ? gloCoords[1] : 0.0);
1486  NekDouble z = (gloCoords.num_elements() > 2 ? gloCoords[2] : 0.0);
1489  GetExp(0)->GetCoordim(), -1, x, y, z);
1490 
1491  // Get the list of elements whose bounding box contains the desired
1492  // point.
1493  std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
1494 
1495  NekDouble nearpt = 1e6;
1496  NekDouble nearpt_min = 1e6;
1497  int min_id = 0;
1498  Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
1499 
1500  // Check each element in turn to see if point lies within it.
1501  for (int i = 0; i < elmts.size(); ++i)
1502  {
1503  if ((*m_exp)[m_elmtToExpId[elmts[i]]]->
1504  GetGeom()->ContainsPoint(gloCoords,
1505  locCoords,
1506  tol, nearpt))
1507  {
1508  return m_elmtToExpId[elmts[i]];
1509  }
1510  else
1511  {
1512  // If it does not lie within, keep track of which element
1513  // is nearest.
1514  if(nearpt < nearpt_min)
1515  {
1516  min_id = m_elmtToExpId[elmts[i]];
1517  nearpt_min = nearpt;
1518  Vmath::Vcopy(locCoords.num_elements(),locCoords, 1,
1519  savLocCoords, 1);
1520  }
1521  }
1522  }
1523 
1524  // If the calling function is with just the nearest element, return
1525  // that. Otherwise return -1 to indicate no matching elemenet found.
1526  if(returnNearestElmt)
1527  {
1528 
1529  std::string msg = "Failed to find point within element to "
1530  "tolerance of "
1531  + boost::lexical_cast<std::string>(tol)
1532  + " using local point ("
1533  + boost::lexical_cast<std::string>(locCoords[0]) +","
1534  + boost::lexical_cast<std::string>(locCoords[1]) +","
1535  + boost::lexical_cast<std::string>(locCoords[1])
1536  + ") in element: "
1537  + boost::lexical_cast<std::string>(min_id);
1538  WARNINGL1(false,msg.c_str());
1539 
1540  Vmath::Vcopy(locCoords.num_elements(),savLocCoords, 1,
1541  locCoords, 1);
1542  return min_id;
1543  }
1544  else
1545  {
1546  return -1;
1547  }
1548  }
1549 
1550  /**
1551  * Given some coordinates, output the expansion field value at that
1552  * point
1553  */
1555  const Array<OneD, const NekDouble> &coords,
1556  const Array<OneD, const NekDouble> &phys)
1557  {
1558  int dim = GetCoordim(0);
1559  ASSERTL0(dim == coords.num_elements(),
1560  "Invalid coordinate dimension.");
1561 
1562  // Grab the element index corresponding to coords.
1563  Array<OneD, NekDouble> xi(dim);
1564  int elmtIdx = GetExpIndex(coords, xi);
1565  ASSERTL0(elmtIdx > 0, "Unable to find element containing point.");
1566 
1567  // Grab that element's physical storage.
1568  Array<OneD, NekDouble> elmtPhys = phys + m_phys_offset[elmtIdx];
1569 
1570  // Evaluate the element at the appropriate point.
1571  return (*m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
1572  }
1573 
1574  /**
1575  * Configures geometric info, such as tangent direction, on each
1576  * expansion.
1577  * @param graph2D Mesh
1578  */
1580  {
1581 
1582  }
1583 
1584  /**
1585  * @brief Reset geometry information, metrics, matrix managers and
1586  * geometry information.
1587  *
1588  * This routine clears all matrix managers and resets all geometry
1589  * information, which allows the geometry information to be dynamically
1590  * updated as the solver is run.
1591  */
1593  {
1594  // Reset matrix managers.
1596  DNekScalMat, LocalRegions::MatrixKey::opLess>::ClearManager();
1597  LibUtilities::NekManager<LocalRegions::MatrixKey,
1599 
1600  // Loop over all elements and reset geometry information.
1601  for (int i = 0; i < m_exp->size(); ++i)
1602  {
1603  (*m_exp)[i]->GetGeom()->Reset(m_graph->GetCurvedEdges(),
1604  m_graph->GetCurvedFaces());
1605  }
1606 
1607  // Loop over all elements and rebuild geometric factors.
1608  for (int i = 0; i < m_exp->size(); ++i)
1609  {
1610  (*m_exp)[i]->Reset();
1611  }
1612  }
1613 
1614  /**
1615  * Write Tecplot Files Header
1616  * @param outfile Output file name.
1617  * @param var variables names
1618  */
1619  void ExpList::v_WriteTecplotHeader(std::ostream &outfile,
1620  std::string var)
1621  {
1622  if (GetNumElmts() == 0)
1623  {
1624  return;
1625  }
1626 
1627  int coordim = GetExp(0)->GetCoordim();
1628  char vars[3] = { 'x', 'y', 'z' };
1629 
1630  if (m_expType == e3DH1D)
1631  {
1632  coordim += 1;
1633  }
1634  else if (m_expType == e3DH2D)
1635  {
1636  coordim += 2;
1637  }
1638 
1639  outfile << "Variables = x";
1640  for (int i = 1; i < coordim; ++i)
1641  {
1642  outfile << ", " << vars[i];
1643  }
1644 
1645  if (var.size() > 0)
1646  {
1647  outfile << ", " << var;
1648  }
1649 
1650  outfile << std::endl << std::endl;
1651  }
1652 
1653  /**
1654  * Write Tecplot Files Zone
1655  * @param outfile Output file name.
1656  * @param expansion Expansion that is considered
1657  */
1658  void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
1659  {
1660  int i, j;
1661  int coordim = GetCoordim(0);
1662  int nPoints = GetTotPoints();
1663  int nBases = (*m_exp)[0]->GetNumBases();
1664  int numBlocks = 0;
1665 
1667 
1668  if (expansion == -1)
1669  {
1670  nPoints = GetTotPoints();
1671 
1672  coords[0] = Array<OneD, NekDouble>(nPoints);
1673  coords[1] = Array<OneD, NekDouble>(nPoints);
1674  coords[2] = Array<OneD, NekDouble>(nPoints);
1675 
1676  GetCoords(coords[0], coords[1], coords[2]);
1677 
1678  for (i = 0; i < m_exp->size(); ++i)
1679  {
1680  int numInt = 1;
1681 
1682  for (j = 0; j < nBases; ++j)
1683  {
1684  numInt *= (*m_exp)[i]->GetNumPoints(j)-1;
1685  }
1686 
1687  numBlocks += numInt;
1688  }
1689  }
1690  else
1691  {
1692  nPoints = (*m_exp)[expansion]->GetTotPoints();
1693 
1694  coords[0] = Array<OneD, NekDouble>(nPoints);
1695  coords[1] = Array<OneD, NekDouble>(nPoints);
1696  coords[2] = Array<OneD, NekDouble>(nPoints);
1697 
1698  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
1699 
1700  numBlocks = 1;
1701  for (j = 0; j < nBases; ++j)
1702  {
1703  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j)-1;
1704  }
1705  }
1706 
1707  if (m_expType == e3DH1D)
1708  {
1709  nBases += 1;
1710  coordim += 1;
1711  int nPlanes = GetZIDs().num_elements();
1712  NekDouble tmp = numBlocks * (nPlanes-1.0) / nPlanes;
1713  numBlocks = (int)tmp;
1714  }
1715  else if (m_expType == e3DH2D)
1716  {
1717  nBases += 2;
1718  coordim += 1;
1719  }
1720 
1721  outfile << "Zone, N=" << nPoints << ", E="
1722  << numBlocks << ", F=FEBlock" ;
1723 
1724  switch(nBases)
1725  {
1726  case 2:
1727  outfile << ", ET=QUADRILATERAL" << std::endl;
1728  break;
1729  case 3:
1730  outfile << ", ET=BRICK" << std::endl;
1731  break;
1732  default:
1733  ASSERTL0(false,"Not set up for this type of output");
1734  break;
1735  }
1736 
1737  // Write out coordinates
1738  for (j = 0; j < coordim; ++j)
1739  {
1740  for (i = 0; i < nPoints; ++i)
1741  {
1742  outfile << coords[j][i] << " ";
1743  if (i % 1000 == 0 && i)
1744  {
1745  outfile << std::endl;
1746  }
1747  }
1748  outfile << std::endl;
1749  }
1750  }
1751 
1752  void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile,
1753  int expansion)
1754  {
1755  int i,j,k,l;
1756  int nbase = (*m_exp)[0]->GetNumBases();
1757  int cnt = 0;
1758 
1759  std::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
1760 
1761  if (expansion != -1)
1762  {
1763  exp = std::shared_ptr<LocalRegions::ExpansionVector>(
1765  (*exp)[0] = (*m_exp)[expansion];
1766  }
1767 
1768  if (nbase == 2)
1769  {
1770  for(i = 0; i < (*exp).size(); ++i)
1771  {
1772  const int np0 = (*exp)[i]->GetNumPoints(0);
1773  const int np1 = (*exp)[i]->GetNumPoints(1);
1774 
1775  for(j = 1; j < np1; ++j)
1776  {
1777  for(k = 1; k < np0; ++k)
1778  {
1779  outfile << cnt + (j-1)*np0 + k << " ";
1780  outfile << cnt + (j-1)*np0 + k+1 << " ";
1781  outfile << cnt + j *np0 + k+1 << " ";
1782  outfile << cnt + j *np0 + k << endl;
1783  }
1784  }
1785 
1786  cnt += np0*np1;
1787  }
1788  }
1789  else if (nbase == 3)
1790  {
1791  for(i = 0; i < (*exp).size(); ++i)
1792  {
1793  const int np0 = (*exp)[i]->GetNumPoints(0);
1794  const int np1 = (*exp)[i]->GetNumPoints(1);
1795  const int np2 = (*exp)[i]->GetNumPoints(2);
1796  const int np01 = np0*np1;
1797 
1798  for(j = 1; j < np2; ++j)
1799  {
1800  for(k = 1; k < np1; ++k)
1801  {
1802  for(l = 1; l < np0; ++l)
1803  {
1804  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l << " ";
1805  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l+1 << " ";
1806  outfile << cnt + (j-1)*np01 + k *np0 + l+1 << " ";
1807  outfile << cnt + (j-1)*np01 + k *np0 + l << " ";
1808  outfile << cnt + j *np01 + (k-1)*np0 + l << " ";
1809  outfile << cnt + j *np01 + (k-1)*np0 + l+1 << " ";
1810  outfile << cnt + j *np01 + k *np0 + l+1 << " ";
1811  outfile << cnt + j *np01 + k *np0 + l << endl;
1812  }
1813  }
1814  }
1815  cnt += np0*np1*np2;
1816  }
1817  }
1818  else
1819  {
1820  ASSERTL0(false,"Not set up for this dimension");
1821  }
1822  }
1823 
1824  /**
1825  * Write Tecplot Files Field
1826  * @param outfile Output file name.
1827  * @param expansion Expansion that is considered
1828  */
1829  void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
1830  {
1831  if (expansion == -1)
1832  {
1833  int totpoints = GetTotPoints();
1834  if(m_physState == false)
1835  {
1837  }
1838 
1839  for(int i = 0; i < totpoints; ++i)
1840  {
1841  outfile << m_phys[i] << " ";
1842  if(i % 1000 == 0 && i)
1843  {
1844  outfile << std::endl;
1845  }
1846  }
1847  outfile << std::endl;
1848 
1849  }
1850  else
1851  {
1852  int nPoints = (*m_exp)[expansion]->GetTotPoints();
1853 
1854  for (int i = 0; i < nPoints; ++i)
1855  {
1856  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
1857  }
1858 
1859  outfile << std::endl;
1860  }
1861  }
1862 
1863  void ExpList::WriteVtkHeader(std::ostream &outfile)
1864  {
1865  outfile << "<?xml version=\"1.0\"?>" << endl;
1866  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1867  << "byte_order=\"LittleEndian\">" << endl;
1868  outfile << " <UnstructuredGrid>" << endl;
1869  }
1870 
1871  void ExpList::WriteVtkFooter(std::ostream &outfile)
1872  {
1873  outfile << " </UnstructuredGrid>" << endl;
1874  outfile << "</VTKFile>" << endl;
1875  }
1876 
1877  void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
1878  {
1879  boost::ignore_unused(outfile, expansion, istrip);
1881  "Routine not implemented for this expansion.");
1882  }
1883 
1884  void ExpList::WriteVtkPieceFooter(std::ostream &outfile, int expansion)
1885  {
1886  boost::ignore_unused(expansion);
1887  outfile << " </PointData>" << endl;
1888  outfile << " </Piece>" << endl;
1889  }
1890 
1891  void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
1892  std::string var)
1893  {
1894  int i;
1895  int nq = (*m_exp)[expansion]->GetTotPoints();
1896 
1897  // printing the fields of that zone
1898  outfile << " <DataArray type=\"Float64\" Name=\""
1899  << var << "\">" << endl;
1900  outfile << " ";
1901  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
1902  for(i = 0; i < nq; ++i)
1903  {
1904  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
1905  }
1906  outfile << endl;
1907  outfile << " </DataArray>" << endl;
1908  }
1909 
1910  /**
1911  * Given a spectral/hp approximation
1912  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1913  * (which should be contained in #m_phys), this function calculates the
1914  * \f$L_\infty\f$ error of this approximation with respect to an exact
1915  * solution. The local distribution of the quadrature points allows an
1916  * elemental evaluation of this operation through the functions
1917  * StdRegions#StdExpansion#Linf.
1918  *
1919  * The exact solution, also evaluated at the quadrature
1920  * points, should be contained in the variable #m_phys of
1921  * the ExpList object \a Sol.
1922  *
1923  * @param soln A 1D array, containing the discrete
1924  * evaluation of the exact solution at the
1925  * quadrature points in its array #m_phys.
1926  * @return The \f$L_\infty\f$ error of the approximation.
1927  */
1929  const Array<OneD, const NekDouble> &inarray,
1930  const Array<OneD, const NekDouble> &soln)
1931  {
1932  NekDouble err = 0.0;
1933 
1934  if (soln == NullNekDouble1DArray)
1935  {
1936  err = Vmath::Vmax(m_npoints, inarray, 1);
1937  }
1938  else
1939  {
1940  for (int i = 0; i < m_npoints; ++i)
1941  {
1942  err = max(err, abs(inarray[i] - soln[i]));
1943  }
1944  }
1945 
1946  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
1947 
1948  return err;
1949  }
1950 
1951  /**
1952  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
1953  * evaluated at the quadrature points (which should be contained in
1954  * #m_phys), this function calculates the \f$L_2\f$ error of this
1955  * approximation with respect to an exact solution. The local
1956  * distribution of the quadrature points allows an elemental evaluation
1957  * of this operation through the functions StdRegions#StdExpansion#L2.
1958  *
1959  * The exact solution, also evaluated at the quadrature points, should
1960  * be contained in the variable #m_phys of the ExpList object \a Sol.
1961  *
1962  * @param Sol An ExpList, containing the discrete
1963  * evaluation of the exact solution at the
1964  * quadrature points in its array #m_phys.
1965  * @return The \f$L_2\f$ error of the approximation.
1966  */
1968  const Array<OneD, const NekDouble> &inarray,
1969  const Array<OneD, const NekDouble> &soln)
1970  {
1971  NekDouble err = 0.0, errl2;
1972  int i;
1973 
1974  if (soln == NullNekDouble1DArray)
1975  {
1976  for (i = 0; i < (*m_exp).size(); ++i)
1977  {
1978  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
1979  err += errl2*errl2;
1980  }
1981  }
1982  else
1983  {
1984  for (i = 0; i < (*m_exp).size(); ++i)
1985  {
1986  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
1987  soln + m_phys_offset[i]);
1988  err += errl2*errl2;
1989  }
1990  }
1991 
1992  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1993 
1994  return sqrt(err);
1995  }
1996 
1998  {
1999  NekDouble err = 0.0;
2000  int i = 0;
2001 
2002  for (i = 0; i < (*m_exp).size(); ++i)
2003  {
2004  err += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
2005  }
2006  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
2007 
2008  return err;
2009  }
2010 
2012  {
2013  NekDouble flux = 0.0;
2014  int i = 0;
2015  int j;
2016 
2017  for (i = 0; i < (*m_exp).size(); ++i)
2018  {
2019  Array<OneD, Array<OneD, NekDouble> > tmp(inarray.num_elements());
2020  for (j = 0; j < inarray.num_elements(); ++j)
2021  {
2022  tmp[j] = Array<OneD, NekDouble>(inarray[j] + m_phys_offset[i]);
2023  }
2024  flux += (*m_exp)[i]->VectorFlux(tmp);
2025  }
2026 
2027  return flux;
2028  }
2029 
2031  {
2033  "This method is not defined or valid for this class type");
2034  Array<OneD, NekDouble> NoEnergy(1,0.0);
2035  return NoEnergy;
2036  }
2037 
2039  {
2041  "This method is not defined or valid for this class type");
2043  return trans;
2044  }
2045 
2047  {
2049  "This method is not defined or valid for this class type");
2050  NekDouble len = 0.0;
2051  return len;
2052  }
2053 
2055  {
2056  boost::ignore_unused(lhom);
2058  "This method is not defined or valid for this class type");
2059  }
2060 
2062  {
2064  "This method is not defined or valid for this class type");
2065  Array<OneD, unsigned int> NoModes(1);
2066  return NoModes;
2067  }
2068 
2070  {
2072  "This method is not defined or valid for this class type");
2073  Array<OneD, unsigned int> NoModes(1);
2074  return NoModes;
2075  }
2076 
2077 
2079  const NekDouble scale,
2080  const Array<OneD, NekDouble> &inarray,
2081  Array<OneD, NekDouble> &outarray)
2082  {
2083  boost::ignore_unused(scale, inarray, outarray);
2085  "This method is not defined or valid for this class type");
2086  }
2087 
2089  const NekDouble scale,
2090  const Array<OneD, NekDouble> &inarray,
2091  Array<OneD, NekDouble> &outarray)
2092  {
2093  boost::ignore_unused(scale, inarray, outarray);
2095  "This method is not defined or valid for this class type");
2096  }
2097 
2099  {
2101  "This method is not defined or valid for this class type");
2102  }
2103 
2105  const std::string &fileName,
2107  const std::string &varName,
2108  const std::shared_ptr<ExpList> locExpList)
2109  {
2110  string varString = fileName.substr(0, fileName.find_last_of("."));
2111  int j, k, len = varString.length();
2112  varString = varString.substr(len-1, len);
2113 
2114  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2115  std::vector<std::vector<NekDouble> > FieldData;
2116 
2117  std::string ft = LibUtilities::FieldIO::GetFileType(fileName, comm);
2119  .CreateInstance(ft, comm, m_session->GetSharedFilesystem());
2120 
2121  f->Import(fileName, FieldDef, FieldData);
2122 
2123  bool found = false;
2124  for (j = 0; j < FieldDef.size(); ++j)
2125  {
2126  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
2127  {
2128  if (FieldDef[j]->m_fields[k] == varName)
2129  {
2130  // Copy FieldData into locExpList
2131  locExpList->ExtractDataToCoeffs(
2132  FieldDef[j], FieldData[j],
2133  FieldDef[j]->m_fields[k],
2134  locExpList->UpdateCoeffs());
2135  found = true;
2136  }
2137  }
2138  }
2139 
2140  ASSERTL0(found, "Could not find variable '"+varName+
2141  "' in file boundary condition "+fileName);
2142  locExpList->BwdTrans_IterPerExp(
2143  locExpList->GetCoeffs(),
2144  locExpList->UpdatePhys());
2145  }
2146 
2147  /**
2148  * Given a spectral/hp approximation
2149  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
2150  * (which should be contained in #m_phys), this function calculates the
2151  * \f$H^1_2\f$ error of this approximation with respect to an exact
2152  * solution. The local distribution of the quadrature points allows an
2153  * elemental evaluation of this operation through the functions
2154  * StdRegions#StdExpansion#H1.
2155  *
2156  * The exact solution, also evaluated at the quadrature points, should
2157  * be contained in the variable #m_phys of the ExpList object \a Sol.
2158  *
2159  * @param soln An 1D array, containing the discrete evaluation
2160  * of the exact solution at the quadrature points.
2161  *
2162  * @return The \f$H^1_2\f$ error of the approximation.
2163  */
2165  const Array<OneD, const NekDouble> &inarray,
2166  const Array<OneD, const NekDouble> &soln)
2167  {
2168  NekDouble err = 0.0, errh1;
2169  int i;
2170 
2171  for (i = 0; i < (*m_exp).size(); ++i)
2172  {
2173  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
2174  soln + m_phys_offset[i]);
2175  err += errh1*errh1;
2176  }
2177 
2178  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
2179 
2180  return sqrt(err);
2181  }
2182 
2183  void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2184  int NumHomoDir,
2186  std::vector<NekDouble> &HomoLen,
2187  bool homoStrips,
2188  std::vector<unsigned int> &HomoSIDs,
2189  std::vector<unsigned int> &HomoZIDs,
2190  std::vector<unsigned int> &HomoYIDs)
2191  {
2192  int startenum = (int) LibUtilities::eSegment;
2193  int endenum = (int) LibUtilities::eHexahedron;
2194  int s = 0;
2196 
2197  ASSERTL1(NumHomoDir == HomoBasis.num_elements(),"Homogeneous basis is not the same length as NumHomoDir");
2198  ASSERTL1(NumHomoDir == HomoLen.size(),"Homogeneous length vector is not the same length as NumHomDir");
2199 
2200  // count number of shapes
2201  switch((*m_exp)[0]->GetShapeDimension())
2202  {
2203  case 1:
2204  startenum = (int) LibUtilities::eSegment;
2205  endenum = (int) LibUtilities::eSegment;
2206  break;
2207  case 2:
2208  startenum = (int) LibUtilities::eTriangle;
2209  endenum = (int) LibUtilities::eQuadrilateral;
2210  break;
2211  case 3:
2212  startenum = (int) LibUtilities::eTetrahedron;
2213  endenum = (int) LibUtilities::eHexahedron;
2214  break;
2215  }
2216 
2217  for(s = startenum; s <= endenum; ++s)
2218  {
2219  std::vector<unsigned int> elementIDs;
2220  std::vector<LibUtilities::BasisType> basis;
2221  std::vector<unsigned int> numModes;
2222  std::vector<std::string> fields;
2223 
2224  bool first = true;
2225  bool UniOrder = true;
2226  int n;
2227 
2228  shape = (LibUtilities::ShapeType) s;
2229 
2230  for(int i = 0; i < (*m_exp).size(); ++i)
2231  {
2232  if((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
2233  {
2234  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
2235  if(first)
2236  {
2237  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
2238  {
2239  basis.push_back((*m_exp)[i]->GetBasis(j)->GetBasisType());
2240  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
2241  }
2242 
2243  // add homogeneous direction details if defined
2244  for(n = 0 ; n < NumHomoDir; ++n)
2245  {
2246  basis.push_back(HomoBasis[n]->GetBasisType());
2247  numModes.push_back(HomoBasis[n]->GetNumModes());
2248  }
2249 
2250  first = false;
2251  }
2252  else
2253  {
2254  ASSERTL0((*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],"Routine is not set up for multiple bases definitions");
2255 
2256  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
2257  {
2258  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
2259  if(numModes[j] != (*m_exp)[i]->GetBasis(j)->GetNumModes())
2260  {
2261  UniOrder = false;
2262  }
2263  }
2264  // add homogeneous direction details if defined
2265  for(n = 0 ; n < NumHomoDir; ++n)
2266  {
2267  numModes.push_back(HomoBasis[n]->GetNumModes());
2268  }
2269  }
2270  }
2271  }
2272 
2273 
2274  if(elementIDs.size() > 0)
2275  {
2278  AllocateSharedPtr(shape, elementIDs, basis,
2279  UniOrder, numModes,fields,
2280  NumHomoDir, HomoLen, homoStrips,
2281  HomoSIDs, HomoZIDs, HomoYIDs);
2282  fielddef.push_back(fdef);
2283  }
2284  }
2285  }
2286 
2287 
2288  //
2289  // Virtual functions
2290  //
2291  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::v_GetFieldDefinitions()
2292  {
2293  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
2294  v_GetFieldDefinitions(returnval);
2295  return returnval;
2296  }
2297 
2298  void ExpList::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
2299  {
2300  GeneralGetFieldDefinitions(fielddef);
2301  }
2302 
2303  //Append the element data listed in elements
2304  //fielddef->m_ElementIDs onto fielddata
2305  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
2306  {
2307  v_AppendFieldData(fielddef,fielddata,m_coeffs);
2308  }
2309 
2310  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
2311  {
2312  int i;
2313  // Determine mapping from element ids to location in
2314  // expansion list
2315  // Determine mapping from element ids to location in
2316  // expansion list
2317  map<int, int> ElmtID_to_ExpID;
2318  for(i = 0; i < (*m_exp).size(); ++i)
2319  {
2320  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2321  }
2322 
2323  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
2324  {
2325  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
2326  int datalen = (*m_exp)[eid]->GetNcoeffs();
2327  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]],&coeffs[m_coeff_offset[eid]]+datalen);
2328  }
2329 
2330  }
2331 
2332  /// Extract the data in fielddata into the coeffs
2335  std::vector<NekDouble> &fielddata,
2336  std::string &field,
2337  Array<OneD, NekDouble> &coeffs)
2338  {
2339  v_ExtractDataToCoeffs(fielddef,fielddata,field,coeffs);
2340  }
2341 
2342  void ExpList::ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2343  {
2344  v_ExtractCoeffsToCoeffs(fromExpList,fromCoeffs,toCoeffs);
2345  }
2346 
2347  /**
2348  * @brief Extract data from raw field data into expansion list.
2349  *
2350  * @param fielddef Field definitions.
2351  * @param fielddata Data for associated field.
2352  * @param field Field variable name.
2353  * @param coeffs Resulting coefficient array.
2354  */
2357  std::vector<NekDouble> &fielddata,
2358  std::string &field,
2359  Array<OneD, NekDouble> &coeffs)
2360  {
2361  int i, expId;
2362  int offset = 0;
2363  int modes_offset = 0;
2364  int datalen = fielddata.size()/fielddef->m_fields.size();
2365 
2366  // Find data location according to field definition
2367  for(i = 0; i < fielddef->m_fields.size(); ++i)
2368  {
2369  if(fielddef->m_fields[i] == field)
2370  {
2371  break;
2372  }
2373  offset += datalen;
2374  }
2375 
2376  ASSERTL0(i != fielddef->m_fields.size(),
2377  "Field (" + field + ") not found in file.");
2378 
2379  if (m_elmtToExpId.size() == 0)
2380  {
2381  // Loop in reverse order so that in case where using a
2382  // Homogeneous expansion it sets geometry ids to first part of
2383  // m_exp list. Otherwise will set to second (complex) expansion
2384  for(i = (*m_exp).size()-1; i >= 0; --i)
2385  {
2386  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2387  }
2388  }
2389 
2390  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
2391  {
2392  // Reset modes_offset in the case where all expansions of
2393  // the same order.
2394  if (fielddef->m_uniOrder == true)
2395  {
2396  modes_offset = 0;
2397  }
2398 
2400  fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
2401 
2402  const int elmtId = fielddef->m_elementIDs[i];
2403  auto eIt = m_elmtToExpId.find(elmtId);
2404 
2405  if (eIt == m_elmtToExpId.end())
2406  {
2407  offset += datalen;
2408  modes_offset += (*m_exp)[0]->GetNumBases();
2409  continue;
2410  }
2411 
2412  expId = eIt->second;
2413 
2414  bool sameBasis = true;
2415  for (int j = 0; j < fielddef->m_basis.size(); ++j)
2416  {
2417  if (fielddef->m_basis[j] != (*m_exp)[expId]->GetBasisType(j))
2418  {
2419  sameBasis = false;
2420  break;
2421  }
2422  }
2423 
2424  if (datalen == (*m_exp)[expId]->GetNcoeffs() && sameBasis)
2425  {
2426  Vmath::Vcopy(datalen, &fielddata[offset], 1,
2427  &coeffs[m_coeff_offset[expId]], 1);
2428  }
2429  else
2430  {
2431  (*m_exp)[expId]->ExtractDataToCoeffs(
2432  &fielddata[offset], fielddef->m_numModes,
2433  modes_offset, &coeffs[m_coeff_offset[expId]],
2434  fielddef->m_basis);
2435  }
2436 
2437  offset += datalen;
2438  modes_offset += (*m_exp)[0]->GetNumBases();
2439  }
2440 
2441  return;
2442  }
2443 
2444  void ExpList::v_ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2445  {
2446  int i;
2447  int offset = 0;
2448 
2449  for(i = 0; i < (*m_exp).size(); ++i)
2450  {
2451  std::vector<unsigned int> nummodes;
2452  vector<LibUtilities::BasisType> basisTypes;
2453  for(int j= 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
2454  {
2455  nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
2456  basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
2457  }
2458 
2459  (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
2460  &toCoeffs[m_coeff_offset[i]],
2461  basisTypes);
2462 
2463  offset += fromExpList->GetExp(i)->GetNcoeffs();
2464  }
2465  }
2466 
2467 
2469  const SpatialDomains::GeomMMF MMFdir,
2470  const Array<OneD, const NekDouble> &CircCentre,
2471  Array<OneD, Array<OneD, NekDouble> > &outarray)
2472  {
2473  int npts;
2474 
2475  int MFdim = 3;
2476  int nq = outarray[0].num_elements()/MFdim;
2477 
2478  // Assume whole array is of same coordinate dimension
2479  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
2480 
2481  Array<OneD, Array<OneD, NekDouble> > MFloc(MFdim*coordim);
2482  // Process each expansion.
2483  for(int i = 0; i < m_exp->size(); ++i)
2484  {
2485  npts = (*m_exp)[i]->GetTotPoints();
2486 
2487  for (int j=0; j< MFdim*coordim; ++j)
2488  {
2489  MFloc[j] = Array<OneD, NekDouble>(npts, 0.0);
2490  }
2491 
2492  // MF from LOCALREGIONS
2493  (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
2494  (*m_exp)[i]->GetPointsKeys(),
2495  MMFdir,
2496  CircCentre,
2497  MFloc );
2498 
2499  // Get the physical data offset for this expansion.
2500  for (int j = 0; j < MFdim; ++j)
2501  {
2502  for (int k = 0; k < coordim; ++k)
2503  {
2504  Vmath::Vcopy(npts,
2505  &MFloc[j*coordim+k][0], 1,
2506  &outarray[j][k*nq+m_phys_offset[i]], 1);
2507  }
2508  }
2509  }
2510 
2511  }
2512 
2513 
2514 
2515  /**
2516  * @brief Generate vector v such that v[i] = scalar1 if i is in the
2517  * element < ElementID. Otherwise, v[i] = scalar2.
2518  *
2519  */
2521  const int ElementID,
2522  const NekDouble scalar1,
2523  const NekDouble scalar2,
2524  Array<OneD, NekDouble> &outarray)
2525  {
2526  int npoints_e;
2527  NekDouble coeff;
2528 
2529  Array<OneD, NekDouble> outarray_e;
2530 
2531  for(int i = 0 ; i < (*m_exp).size(); ++i)
2532  {
2533  npoints_e = (*m_exp)[i]->GetTotPoints();
2534 
2535  if(i <= ElementID)
2536  {
2537  coeff = scalar1;
2538  }
2539  else
2540  {
2541  coeff = scalar2;
2542  }
2543 
2544  outarray_e = Array<OneD, NekDouble>(npoints_e, coeff);
2545  Vmath::Vcopy(npoints_e, &outarray_e[0], 1,
2546  &outarray[m_phys_offset[i]], 1);
2547  }
2548  }
2549 
2552  {
2553  ASSERTL0(false,
2554  "This method is not defined or valid for this class type");
2556  return result;
2557  }
2558 
2559  std::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
2560  {
2561  boost::ignore_unused(i);
2563  "This method is not defined or valid for this class type");
2564  static std::shared_ptr<ExpList> result;
2565  return result;
2566  }
2567 
2569  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
2570  const Array<OneD, const NekDouble> &Fwd,
2571  const Array<OneD, const NekDouble> &Bwd,
2573  {
2574  boost::ignore_unused(Vec, Fwd, Bwd, Upwind);
2576  "This method is not defined or valid for this class type");
2577  }
2578 
2580  const Array<OneD, const NekDouble> &Vn,
2581  const Array<OneD, const NekDouble> &Fwd,
2582  const Array<OneD, const NekDouble> &Bwd,
2584  {
2585  boost::ignore_unused(Vn, Fwd, Bwd, Upwind);
2587  "This method is not defined or valid for this class type");
2588  }
2589 
2590  std::shared_ptr<ExpList> &ExpList::v_GetTrace()
2591  {
2593  "This method is not defined or valid for this class type");
2594  static std::shared_ptr<ExpList> returnVal;
2595  return returnVal;
2596  }
2597 
2598  std::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
2599  {
2601  "This method is not defined or valid for this class type");
2602  static std::shared_ptr<AssemblyMapDG> result;
2603  return result;
2604  }
2605 
2607  {
2608  return GetTraceMap()->GetBndCondTraceToGlobalTraceMap();
2609  }
2610 
2612  Array<OneD, Array<OneD, NekDouble> > &normals)
2613  {
2614  boost::ignore_unused(normals);
2616  "This method is not defined or valid for this class type");
2617  }
2618 
2620  const Array<OneD, const NekDouble> &Fx,
2621  const Array<OneD, const NekDouble> &Fy,
2622  Array<OneD, NekDouble> &outarray)
2623  {
2624  boost::ignore_unused(Fx, Fy, outarray);
2626  "This method is not defined or valid for this class type");
2627  }
2628 
2630  const Array<OneD, const NekDouble> &Fn,
2631  Array<OneD, NekDouble> &outarray)
2632  {
2633  boost::ignore_unused(Fn, outarray);
2635  "This method is not defined or valid for this class type");
2636  }
2637 
2639  const Array<OneD, const NekDouble> &Fwd,
2640  const Array<OneD, const NekDouble> &Bwd,
2641  Array<OneD, NekDouble> &outarray)
2642  {
2643  boost::ignore_unused(Fwd, Bwd, outarray);
2645  "This method is not defined or valid for this class type");
2646  }
2647 
2649  Array<OneD,NekDouble> &Bwd)
2650  {
2651  boost::ignore_unused(Fwd, Bwd);
2653  "This method is not defined or valid for this class type");
2654  }
2655 
2657  const Array<OneD,const NekDouble> &field,
2658  Array<OneD,NekDouble> &Fwd,
2659  Array<OneD,NekDouble> &Bwd)
2660  {
2661  boost::ignore_unused(field, Fwd, Bwd);
2663  "This method is not defined or valid for this class type");
2664  }
2665 
2666  const vector<bool> &ExpList::v_GetLeftAdjacentFaces(void) const
2667  {
2669  "This method is not defined or valid for this class type");
2670  static vector<bool> tmp;
2671  return tmp;
2672  }
2673 
2674 
2676  {
2677  boost::ignore_unused(outarray);
2679  "This method is not defined or valid for this class type");
2680  }
2681 
2683  const Array<OneD, const NekDouble> &inarray,
2684  Array<OneD,NekDouble> &outarray)
2685  {
2686  boost::ignore_unused(inarray, outarray);
2688  "This method is not defined or valid for this class type");
2689  }
2690 
2692  const Array<OneD,const NekDouble> &inarray,
2693  Array<OneD, NekDouble> &outarray,
2694  CoeffState coeffstate)
2695  {
2696  boost::ignore_unused(inarray, outarray, coeffstate);
2698  "This method is not defined or valid for this class type");
2699  }
2700 
2702  const Array<OneD, const NekDouble> &inarray,
2703  Array<OneD, NekDouble> &outarray,
2704  const FlagList &flags,
2705  const StdRegions::ConstFactorMap &factors,
2706  const StdRegions::VarCoeffMap &varcoeff,
2707  const MultiRegions::VarFactorsMap &varfactors,
2708  const Array<OneD, const NekDouble> &dirForcing,
2709  const bool PhysSpaceForcing)
2710  {
2711  boost::ignore_unused(inarray, outarray, flags, factors, varcoeff,
2712  varfactors, dirForcing, PhysSpaceForcing);
2713  NEKERROR(ErrorUtil::efatal, "HelmSolve not implemented.");
2714  }
2715 
2717  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2718  const Array<OneD, const NekDouble> &inarray,
2719  Array<OneD, NekDouble> &outarray,
2720  const NekDouble lambda,
2721  CoeffState coeffstate,
2722  const Array<OneD, const NekDouble>& dirForcing)
2723  {
2724  boost::ignore_unused(velocity, inarray, outarray, lambda,
2725  coeffstate, dirForcing);
2727  "This method is not defined or valid for this class type");
2728  }
2729 
2731  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2732  const Array<OneD, const NekDouble> &inarray,
2733  Array<OneD, NekDouble> &outarray,
2734  const NekDouble lambda,
2735  CoeffState coeffstate,
2736  const Array<OneD, const NekDouble>& dirForcing)
2737  {
2738  boost::ignore_unused(velocity, inarray, outarray, lambda,
2739  coeffstate, dirForcing);
2741  "This method is not defined or valid for this class type");
2742  }
2743 
2745  Array<OneD, NekDouble> &outarray,
2746  CoeffState coeffstate,
2747  bool Shuff,
2748  bool UnShuff)
2749  {
2750  boost::ignore_unused(inarray, outarray, coeffstate, Shuff, UnShuff);
2752  "This method is not defined or valid for this class type");
2753  }
2754 
2756  Array<OneD, NekDouble> &outarray,
2757  CoeffState coeffstate,
2758  bool Shuff,
2759  bool UnShuff)
2760  {
2761  boost::ignore_unused(inarray, outarray, coeffstate, Shuff, UnShuff);
2763  "This method is not defined or valid for this class type");
2764  }
2765 
2767  const Array<OneD, NekDouble> &inarray2,
2768  Array<OneD, NekDouble> &outarray,
2769  CoeffState coeffstate)
2770  {
2771  boost::ignore_unused(inarray1, inarray2, outarray, coeffstate);
2773  "This method is not defined or valid for this class type");
2774  }
2775 
2777  const Array<OneD, Array<OneD, NekDouble> > &inarray1,
2778  const Array<OneD, Array<OneD, NekDouble> > &inarray2,
2779  Array<OneD, Array<OneD, NekDouble> > &outarray,
2780  CoeffState coeffstate)
2781  {
2782  boost::ignore_unused(inarray1, inarray2, outarray, coeffstate);
2784  "This method is not defined or valid for this class type");
2785  }
2786 
2788  const Array<OneD, NekDouble> &TotField,
2789  int BndID)
2790  {
2791  boost::ignore_unused(BndVals, TotField, BndID);
2793  "This method is not defined or valid for this class type");
2794  }
2795 
2798  Array<OneD, NekDouble> &outarray,
2799  int BndID)
2800  {
2801  boost::ignore_unused(V1, V2, outarray, BndID);
2803  "This method is not defined or valid for this class type");
2804  }
2805 
2808  Array<OneD, NekDouble> &outarray)
2809  {
2811  switch (GetCoordim(0))
2812  {
2813  case 1:
2814  {
2815  for(int i = 0; i < GetExpSize(); ++i)
2816  {
2817  (*m_exp)[i]->NormVectorIProductWRTBase(
2818  V[0] + GetPhys_Offset(i),
2819  tmp = outarray + GetCoeff_Offset(i));
2820  }
2821  }
2822  break;
2823  case 2:
2824  {
2825  for(int i = 0; i < GetExpSize(); ++i)
2826  {
2827  (*m_exp)[i]->NormVectorIProductWRTBase(
2828  V[0] + GetPhys_Offset(i),
2829  V[1] + GetPhys_Offset(i),
2830  tmp = outarray + GetCoeff_Offset(i));
2831  }
2832  }
2833  break;
2834  case 3:
2835  {
2836  for(int i = 0; i < GetExpSize(); ++i)
2837  {
2838  (*m_exp)[i]->NormVectorIProductWRTBase(
2839  V[0] + GetPhys_Offset(i),
2840  V[1] + GetPhys_Offset(i),
2841  V[2] + GetPhys_Offset(i),
2842  tmp = outarray + GetCoeff_Offset(i));
2843  }
2844  }
2845  break;
2846  default:
2847  ASSERTL0(false,"Dimension not supported");
2848  break;
2849  }
2850  }
2851 
2853  {
2854  boost::ignore_unused(outarray);
2856  "This method is not defined or valid for this class type");
2857  }
2858 
2859  /**
2860  */
2862  {
2864  "This method is not defined or valid for this class type");
2865  }
2866 
2867  /**
2868  */
2869  void ExpList::v_FillBndCondFromField(const int nreg)
2870  {
2871  boost::ignore_unused(nreg);
2873  "This method is not defined or valid for this class type");
2874  }
2875 
2876  void ExpList::v_LocalToGlobal(bool useComm)
2877  {
2878  boost::ignore_unused(useComm);
2880  "This method is not defined or valid for this class type");
2881  }
2882 
2883 
2885  Array<OneD,NekDouble> &outarray,
2886  bool useComm)
2887  {
2888  boost::ignore_unused(inarray, outarray, useComm);
2890  "This method is not defined or valid for this class type");
2891  }
2892 
2893 
2895  {
2897  "This method is not defined or valid for this class type");
2898  }
2899 
2900 
2902  Array<OneD,NekDouble> &outarray)
2903  {
2904  boost::ignore_unused(inarray, outarray);
2906  "This method is not defined or valid for this class type");
2907  }
2908 
2909 
2911  Array<OneD, NekDouble> &outarray,
2912  CoeffState coeffstate)
2913  {
2914  boost::ignore_unused(coeffstate);
2915  v_BwdTrans_IterPerExp(inarray,outarray);
2916  }
2917 
2919  Array<OneD, NekDouble> &outarray,
2920  CoeffState coeffstate)
2921  {
2922  boost::ignore_unused(coeffstate);
2923  v_FwdTrans_IterPerExp(inarray,outarray);
2924  }
2925 
2927  const Array<OneD, const NekDouble> &inarray,
2928  Array<OneD, NekDouble> &outarray,
2929  CoeffState coeffstate)
2930  {
2931  boost::ignore_unused(coeffstate);
2933  for (int i = 0; i < m_collections.size(); ++i)
2934  {
2935 
2937  inarray + m_coll_phys_offset[i],
2938  tmp = outarray + m_coll_coeff_offset[i]);
2939  }
2940  }
2941 
2943  const GlobalMatrixKey &gkey,
2944  const Array<OneD,const NekDouble> &inarray,
2945  Array<OneD, NekDouble> &outarray,
2946  CoeffState coeffstate)
2947  {
2948  boost::ignore_unused(coeffstate);
2949  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
2950  }
2951 
2952  /**
2953  * The operation is evaluated locally by the elemental
2954  * function StdRegions#StdExpansion#GetCoords.
2955  *
2956  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
2957  * will be stored in this array.
2958  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
2959  * will be stored in this array.
2960  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
2961  * will be stored in this array.
2962  */
2964  Array<OneD, NekDouble> &coord_1,
2965  Array<OneD, NekDouble> &coord_2)
2966  {
2967  if (GetNumElmts() == 0)
2968  {
2969  return;
2970  }
2971 
2972  int i;
2973  Array<OneD, NekDouble> e_coord_0;
2974  Array<OneD, NekDouble> e_coord_1;
2975  Array<OneD, NekDouble> e_coord_2;
2976 
2977  switch(GetExp(0)->GetCoordim())
2978  {
2979  case 1:
2980  for(i= 0; i < (*m_exp).size(); ++i)
2981  {
2982  e_coord_0 = coord_0 + m_phys_offset[i];
2983  (*m_exp)[i]->GetCoords(e_coord_0);
2984  }
2985  break;
2986  case 2:
2987  ASSERTL0(coord_1.num_elements() != 0,
2988  "output coord_1 is not defined");
2989 
2990  for(i= 0; i < (*m_exp).size(); ++i)
2991  {
2992  e_coord_0 = coord_0 + m_phys_offset[i];
2993  e_coord_1 = coord_1 + m_phys_offset[i];
2994  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1);
2995  }
2996  break;
2997  case 3:
2998  ASSERTL0(coord_1.num_elements() != 0,
2999  "output coord_1 is not defined");
3000  ASSERTL0(coord_2.num_elements() != 0,
3001  "output coord_2 is not defined");
3002 
3003  for(i= 0; i < (*m_exp).size(); ++i)
3004  {
3005  e_coord_0 = coord_0 + m_phys_offset[i];
3006  e_coord_1 = coord_1 + m_phys_offset[i];
3007  e_coord_2 = coord_2 + m_phys_offset[i];
3008  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1,e_coord_2);
3009  }
3010  break;
3011  }
3012  }
3013 
3014  /**
3015  */
3017  {
3019  "This method is not defined or valid for this class type");
3020  }
3021 
3022  /**
3023  */
3025  std::shared_ptr<ExpList> &result,
3026  const bool DeclareCoeffPhysArrays)
3027  {
3028  boost::ignore_unused(i, result, DeclareCoeffPhysArrays);
3030  "This method is not defined or valid for this class type");
3031  }
3032 
3033  /**
3034  */
3036  const Array<OneD, NekDouble> &element,
3037  Array<OneD, NekDouble> &boundary)
3038  {
3039  int n, cnt;
3040  Array<OneD, NekDouble> tmp1, tmp2;
3042 
3043  Array<OneD, int> ElmtID,EdgeID;
3044  GetBoundaryToElmtMap(ElmtID,EdgeID);
3045 
3046  // Initialise result
3047  boundary = Array<OneD, NekDouble>
3048  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3049 
3050  // Skip other boundary regions
3051  for (cnt = n = 0; n < i; ++n)
3052  {
3053  cnt += GetBndCondExpansions()[n]->GetExpSize();
3054  }
3055 
3056  int offsetBnd;
3057  int offsetElmt = 0;
3058  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3059  {
3060  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3061 
3062  elmt = GetExp(ElmtID[cnt+n]);
3063  elmt->GetTracePhysVals(EdgeID[cnt+n],
3064  GetBndCondExpansions()[i]->GetExp(n),
3065  tmp1 = element + offsetElmt,
3066  tmp2 = boundary + offsetBnd);
3067 
3068  offsetElmt += elmt->GetTotPoints();
3069  }
3070  }
3071 
3072  /**
3073  */
3075  const Array<OneD, const NekDouble> &phys,
3076  Array<OneD, NekDouble> &bndElmt)
3077  {
3078  int n, cnt, nq;
3079 
3080  Array<OneD, int> ElmtID,EdgeID;
3081  GetBoundaryToElmtMap(ElmtID,EdgeID);
3082 
3083  // Skip other boundary regions
3084  for (cnt = n = 0; n < i; ++n)
3085  {
3086  cnt += GetBndCondExpansions()[n]->GetExpSize();
3087  }
3088 
3089  // Count number of points
3090  int npoints = 0;
3091  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3092  {
3093  npoints += GetExp(ElmtID[cnt+n])->GetTotPoints();
3094  }
3095 
3096  // Initialise result
3097  bndElmt = Array<OneD, NekDouble> (npoints, 0.0);
3098 
3099  // Extract data
3100  int offsetPhys;
3101  int offsetElmt = 0;
3102  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3103  {
3104  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
3105  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
3106  Vmath::Vcopy(nq, &phys[offsetPhys], 1,
3107  &bndElmt[offsetElmt], 1);
3108  offsetElmt += nq;
3109  }
3110  }
3111 
3112  /**
3113  */
3115  const Array<OneD, const NekDouble> &phys,
3117  {
3118  int n, cnt;
3121 
3122  Array<OneD, int> ElmtID,EdgeID;
3123  GetBoundaryToElmtMap(ElmtID,EdgeID);
3124 
3125  // Initialise result
3127  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3128 
3129  // Skip other boundary regions
3130  for (cnt = n = 0; n < i; ++n)
3131  {
3132  cnt += GetBndCondExpansions()[n]->GetExpSize();
3133  }
3134 
3135  int offsetBnd;
3136  int offsetPhys;
3137  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3138  {
3139  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
3140  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3141 
3142  elmt = GetExp(ElmtID[cnt+n]);
3143  elmt->GetTracePhysVals(EdgeID[cnt+n],
3144  GetBndCondExpansions()[i]->GetExp(n),
3145  phys + offsetPhys,
3146  tmp1 = bnd + offsetBnd);
3147  }
3148  }
3149 
3150  /**
3151  */
3153  Array<OneD, Array<OneD, NekDouble> > &normals)
3154  {
3155  int j, n, cnt, nq;
3156  int coordim = GetCoordim(0);
3159 
3160  Array<OneD, int> ElmtID,EdgeID;
3161  GetBoundaryToElmtMap(ElmtID,EdgeID);
3162 
3163  // Initialise result
3164  normals = Array<OneD, Array<OneD, NekDouble> > (coordim);
3165  for (j = 0; j < coordim; ++j)
3166  {
3167  normals[j] = Array<OneD, NekDouble> (
3168  GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3169  }
3170 
3171  // Skip other boundary regions
3172  for (cnt = n = 0; n < i; ++n)
3173  {
3174  cnt += GetBndCondExpansions()[n]->GetExpSize();
3175  }
3176 
3177  int offset;
3178  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3179  {
3180  offset = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3181  nq = GetBndCondExpansions()[i]->GetExp(n)->GetTotPoints();
3182 
3183  elmt = GetExp(ElmtID[cnt+n]);
3184  const Array<OneD, const Array<OneD, NekDouble> > normalsElmt
3185  = elmt->GetSurfaceNormal(EdgeID[cnt+n]);
3186  // Copy to result
3187  for (j = 0; j < coordim; ++j)
3188  {
3189  Vmath::Vcopy(nq, normalsElmt[j], 1,
3190  tmp = normals[j] + offset, 1);
3191  }
3192  }
3193  }
3194 
3195  /**
3196  */
3198  Array<OneD,int> &EdgeID)
3199  {
3200  boost::ignore_unused(ElmtID, EdgeID);
3202  "This method is not defined or valid for this class type");
3203  }
3204 
3205  /**
3206  */
3208  {
3210  "This method is not defined or valid for this class type");
3211  }
3212 
3213  /**
3214  */
3217  {
3219  "This method is not defined or valid for this class type");
3221  result;
3222  return result;
3223  }
3224 
3225  /**
3226  */
3228  {
3230  "This method is not defined or valid for this class type");
3232  return result;
3233  }
3234 
3235  /**
3236  */
3238  const NekDouble time,
3239  const std::string varName,
3240  const NekDouble x2_in,
3241  const NekDouble x3_in)
3242  {
3243  boost::ignore_unused(time, varName, x2_in, x3_in);
3245  "This method is not defined or valid for this class type");
3246  }
3247 
3248  /**
3249  */
3250  map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
3251  {
3253  "This method is not defined or valid for this class type");
3254  static map<int,RobinBCInfoSharedPtr> result;
3255  return result;
3256  }
3257 
3258  /**
3259  */
3261  PeriodicMap &periodicVerts,
3262  PeriodicMap &periodicEdges,
3263  PeriodicMap &periodicFaces)
3264  {
3265  boost::ignore_unused(periodicVerts, periodicEdges, periodicFaces);
3267  "This method is not defined or valid for this class type");
3268  }
3269 
3272  unsigned int regionId,
3273  const std::string& variable)
3274  {
3275  auto collectionIter = collection.find(regionId);
3276  ASSERTL1(collectionIter != collection.end(),
3277  "Unable to locate collection " +
3278  boost::lexical_cast<string>(regionId));
3279 
3281  = (*collectionIter).second;
3282  auto conditionMapIter = bndCondMap->find(variable);
3283  ASSERTL1(conditionMapIter != bndCondMap->end(),
3284  "Unable to locate condition map.");
3285 
3286  const SpatialDomains::BoundaryConditionShPtr boundaryCondition
3287  = (*conditionMapIter).second;
3288 
3289  return boundaryCondition;
3290  }
3291 
3293  {
3294  boost::ignore_unused(n);
3296  "This method is not defined or valid for this class type");
3297  return NullExpListSharedPtr;
3298  }
3299 
3300  /**
3301  * @brief Construct collections of elements containing a single element
3302  * type and polynomial order from the list of expansions.
3303  */
3305  {
3307  vector<std::pair<LocalRegions::ExpansionSharedPtr,int> > > collections;
3308 
3309  // Figure out optimisation parameters if provided in
3310  // session file or default given
3312  ImpType = colOpt.GetDefaultImplementationType();
3313 
3314  bool autotuning = colOpt.IsUsingAutotuning();
3315  bool verbose = (m_session->DefinesCmdLineArgument("verbose")) &&
3316  (m_comm->GetRank() == 0);
3317  int collmax = (colOpt.GetMaxCollectionSize() > 0
3318  ? colOpt.GetMaxCollectionSize()
3319  : 2*m_exp->size());
3320 
3321  // clear vectors in case previously called
3322  m_collections.clear();
3323  m_coll_coeff_offset.clear();
3324  m_coll_phys_offset.clear();
3325 
3326  // Loop over expansions, and create collections for each element type
3327  for (int i = 0; i < m_exp->size(); ++i)
3328  {
3329  collections[(*m_exp)[i]->DetShapeType()].push_back(
3330  std::pair<LocalRegions::ExpansionSharedPtr,int> ((*m_exp)[i],i));
3331  }
3332 
3333  for (auto &it : collections)
3334  {
3335  LocalRegions::ExpansionSharedPtr exp = it.second[0].first;
3336 
3337  Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(exp);
3338  vector<StdRegions::StdExpansionSharedPtr> collExp;
3339 
3340  int prevCoeffOffset = m_coeff_offset[it.second[0].second];
3341  int prevPhysOffset = m_phys_offset [it.second[0].second];
3342  int collcnt;
3343 
3344  m_coll_coeff_offset.push_back(prevCoeffOffset);
3345  m_coll_phys_offset .push_back(prevPhysOffset);
3346 
3347  if(it.second.size() == 1) // single element case
3348  {
3349  collExp.push_back(it.second[0].first);
3350 
3351  // if no Imp Type provided and No settign in xml file.
3352  // reset impTypes using timings
3353  if(autotuning)
3354  {
3355  impTypes = colOpt.SetWithTimings(collExp,
3356  impTypes, verbose);
3357  }
3358 
3359  Collections::Collection tmp(collExp, impTypes);
3360  m_collections.push_back(tmp);
3361  }
3362  else
3363  {
3364  // set up first geometry
3365  collExp.push_back(it.second[0].first);
3366  int prevnCoeff = it.second[0].first->GetNcoeffs();
3367  int prevnPhys = it.second[0].first->GetTotPoints();
3368  collcnt = 1;
3369 
3370  for (int i = 1; i < it.second.size(); ++i)
3371  {
3372  int nCoeffs = it.second[i].first->GetNcoeffs();
3373  int nPhys = it.second[i].first->GetTotPoints();
3374  int coeffOffset = m_coeff_offset[it.second[i].second];
3375  int physOffset = m_phys_offset [it.second[i].second];
3376 
3377  // check to see if next elmt is different or
3378  // collmax reached and if so end collection
3379  // and start new one
3380  if(prevCoeffOffset + nCoeffs != coeffOffset ||
3381  prevnCoeff != nCoeffs ||
3382  prevPhysOffset + nPhys != physOffset ||
3383  prevnPhys != nPhys || collcnt >= collmax)
3384  {
3385 
3386  // if no Imp Type provided and No
3387  // settign in xml file. reset
3388  // impTypes using timings
3389  if(autotuning)
3390  {
3391  impTypes = colOpt.SetWithTimings(collExp,
3392  impTypes,
3393  verbose);
3394  }
3395 
3396  Collections::Collection tmp(collExp, impTypes);
3397  m_collections.push_back(tmp);
3398 
3399 
3400  // start new geom list
3401  collExp.clear();
3402 
3403  m_coll_coeff_offset.push_back(coeffOffset);
3404  m_coll_phys_offset .push_back(physOffset);
3405  collExp.push_back(it.second[i].first);
3406  collcnt = 1;
3407  }
3408  else // add to list of collections
3409  {
3410  collExp.push_back(it.second[i].first);
3411  collcnt++;
3412  }
3413 
3414  // if end of list finish up collection
3415  if (i == it.second.size() - 1)
3416  {
3417  // if no Imp Type provided and No
3418  // settign in xml file.
3419  if(autotuning)
3420  {
3421  impTypes = colOpt.SetWithTimings(collExp,
3422  impTypes,verbose);
3423  }
3424 
3425  Collections::Collection tmp(collExp, impTypes);
3426  m_collections.push_back(tmp);
3427  collExp.clear();
3428  collcnt = 0;
3429 
3430  }
3431 
3432  prevCoeffOffset = coeffOffset;
3433  prevPhysOffset = physOffset;
3434  prevnCoeff = nCoeffs;
3435  prevnPhys = nPhys;
3436  }
3437  }
3438  }
3439  }
3440 
3442  {
3444  }
3445 
3446  } //end of namespace
3447 } //end of namespace
3448 
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Definition: ExpList.cpp:1967
virtual void v_GlobalToLocal(void)
Definition: ExpList.cpp:2894
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2638
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
Definition: ExpList.cpp:826
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:553
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, v[i] = scalar2.
Definition: ExpList.cpp:2520
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2942
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
Definition: ExpList.cpp:1004
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:2355
COLLECTIONS_EXPORT OperatorImpMap SetWithTimings(std::vector< StdRegions::StdExpansionSharedPtr > pGeom, OperatorImpMap &impTypes, bool verbose=true)
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:429
int GetExpIndex(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false)
Definition: ExpList.cpp:1451
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1592
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:2342
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:163
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mod...
Definition: ErrorUtil.hpp:209
virtual void v_GetBoundaryNormals(int i, Array< OneD, Array< OneD, NekDouble > > &normals)
Definition: ExpList.cpp:3152
virtual std::shared_ptr< ExpList > & v_GetTrace()
Definition: ExpList.cpp:2590
void GeneralMatrixOp_IterPerExp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1019
ExpList()
The default constructor.
Definition: ExpList.cpp:96
static ExpListSharedPtr NullExpListSharedPtr
Definition: ExpList.h:1550
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
Definition: AssemblyMapCG.h:51
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1090
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:697
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
Definition: ExpList.h:1106
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::map< CoordType, NekDouble > COOMatType
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Definition: ExpList.h:968
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:782
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:73
virtual void v_FillBndCondFromField()
Definition: ExpList.cpp:2861
void ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: ExpList.cpp:729
NekDouble PhysIntegral(void)
This function integrates a function over the domain consisting of all the elements of the expansion...
Definition: ExpList.cpp:304
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:1658
ExpansionType GetExpType(void)
Returns the type of the expansion.
Definition: ExpList.cpp:271
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2926
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
virtual void v_BwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:1418
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:2183
Lagrange Polynomials using the Gauss points .
Definition: BasisType.h:55
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
Definition: GlobalLinSys.h:50
STL namespace.
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.h:1775
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Definition: ExpList.cpp:2611
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:589
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:294
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.h:2409
BlockMatrixMapShPtr m_blockMat
Definition: ExpList.h:1108
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1069
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:69
GeomMMF
Principle direction for MMF.
void ApplyGeomInfo()
Apply geometry information to each expansion.
Definition: ExpList.cpp:1579
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1052
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
Definition: ExpList.cpp:1877
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
Definition: ExpList.cpp:2796
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
Definition: ExpList.cpp:2038
int GetExpSize(void)
This function returns the number of elements in the expansion.
Definition: ExpList.h:2170
virtual void v_FwdTrans_BndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:804
void PhysDeriv(Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
Definition: ExpList.h:1897
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:71
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
Definition: ExpList.h:2191
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:1381
std::vector< ExpansionSharedPtr > ExpansionVector
Definition: Expansion.h:67
virtual void v_ExtractPhysToBndElmt(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
Definition: ExpList.cpp:3074
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:755
virtual void v_WriteTecplotHeader(std::ostream &outfile, std::string var="")
Definition: ExpList.cpp:1619
virtual NekDouble v_GetHomoLen(void)
Definition: ExpList.cpp:2046
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
Definition: ExpList.cpp:3114
void WriteVtkFooter(std::ostream &outfile)
Definition: ExpList.cpp:1871
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
Definition: ExpList.cpp:2963
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
Definition: AssemblyMap.h:52
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset)
Definition: ShapeType.hpp:313
void 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.h:2249
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
Definition: ExpList.cpp:2030
virtual const Array< OneD, const int > & v_GetTraceBndMap()
Definition: ExpList.cpp:2606
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
Definition: ExpList.cpp:1891
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Definition: ExpList.cpp:3024
const char *const GlobalSysSolnTypeMap[]
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
Definition: ExpList.cpp:2787
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:225
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Definition: ExpList.cpp:2011
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual void v_LocalToGlobal(bool UseComm)
Definition: ExpList.cpp:2876
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1101
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
int GetNumElmts(void)
This function returns the number of elements in the expansion which may be different for a homogeoeno...
Definition: ExpList.h:655
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:103
std::map< OperatorType, ImplementationType > OperatorImpMap
Definition: Operator.h:103
static const NekDouble kNekZeroTol
std::shared_ptr< AssemblyMapDG > & GetTraceMap(void)
Definition: ExpList.h:2272
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:1854
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2852
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:216
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
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:2568
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:1928
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
Definition: NekTypeDefs.hpp:65
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:1086
virtual void v_FwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:794
virtual ~ExpList()
The default destructor.
Definition: ExpList.cpp:284
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:409
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
Definition: StdRegions.hpp:264
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
Definition: GlobalMatrix.h:88
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1104
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
Definition: FieldIO.cpp:96
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:1725
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Definition: ExpList.cpp:1997
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
Definition: ExpList.cpp:2291
bool m_physState
The state of the array m_phys.
Definition: ExpList.h:1078
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:1752
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
Definition: ExpList.h:1030
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:2744
void SetCoeffPhysOffsets()
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
Definition: ExpList.cpp:249
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Definition: ExpList.h:1026
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:3237
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
Definition: ExpList.cpp:1554
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1023
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Definition: ExpList.cpp:2444
LibUtilities::ShapeType GetShapeType() const
Return the expansion type associated with key.
virtual void v_GetMovingFrames(const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: ExpList.cpp:2468
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:399
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap()
Definition: ExpList.cpp:2598
double NekDouble
virtual void v_SetHomoLen(const NekDouble lhom)
Definition: ExpList.cpp:2054
const DNekScalBlkMatSharedPtr GenBlockMatrix(const GlobalMatrixKey &gkey)
This function assembles the block diagonal matrix of local matrices of the type mtype.
Definition: ExpList.cpp:869
Defines a list of flags.
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
Definition: Conditions.h:226
std::shared_ptr< PointGeom > PointGeomSharedPtr
Definition: Geometry.h:59
virtual const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions(void)
Definition: ExpList.cpp:2551
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
Definition: ExpList.cpp:3227
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:2104
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_ClearGlobalLinSysManager(void)
Definition: ExpList.cpp:2098
Describe a linear system.
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
Definition: Conditions.h:219
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.cpp:2766
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
Definition: ExpList.cpp:3292
Describes a matrix with ordering defined by a local to global map.
std::map< int, RobinBCInfoSharedPtr > GetRobinBCInfo()
Definition: ExpList.h:884
const StdRegions::VarCoeffMap & GetVarCoeffs() const
virtual void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
Definition: ExpList.cpp:3260
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:346
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2078
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions()
Definition: ExpList.h:2239
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:65
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:88
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:125
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1020
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2619
virtual void v_CurlCurl(Array< OneD, Array< OneD, NekDouble > > &Vel, Array< OneD, Array< OneD, NekDouble > > &Q)
Definition: ExpList.cpp:625
virtual void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, CoeffState coeffstate=eLocal, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:2730
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
Definition: ExpList.cpp:2061
std::map< GlobalMatrixKey, DNekScalBlkMatSharedPtr > BlockMatrixMap
A map between global matrix keys and their associated block matrices.
Definition: ExpList.h:97
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo(void)
Definition: ExpList.cpp:3250
#define WARNINGL1(condition, msg)
Definition: ErrorUtil.hpp:251
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray, CoeffState coeffstate=eLocal)
Definition: ExpList.cpp:2776
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:2164
virtual void v_IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:382
void WriteVtkPieceFooter(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:1884
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:179
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
Definition: ExpList.cpp:2666
virtual void v_LinearAdvectionDiffusionReactionSolve(const Array< OneD, Array< OneD, NekDouble > > &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, CoeffState coeffstate=eLocal, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
Definition: ExpList.cpp:2716
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:1362
std::vector< int > m_coll_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1098
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
Definition: ExpList.cpp:2305
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:3270
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:3304
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Definition: ExpList.cpp:1222
void MultiplyByBlockMatrix(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:352
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2691
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1116
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
Definition: ExpList.cpp:3197
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
Definition: ExpList.cpp:2069
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:67
virtual void v_SetUpPhysNormals()
Definition: ExpList.cpp:3016
std::shared_ptr< Transposition > TranspositionSharedPtr
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Definition: ExpList.h:2013
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions()
Definition: ExpList.cpp:3216
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
Definition: ExpList.cpp:2559
Lagrange for SEM basis .
Definition: BasisType.h:54
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
Definition: ExpList.cpp:279
virtual void v_ExtractElmtToBndPhys(const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
Definition: ExpList.cpp:3035
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:2200
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
Definition: ExpList.h:1608
GlobalLinSysFactory & GetGlobalLinSysFactory()
std::vector< int > m_coll_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1095
void WriteVtkHeader(std::ostream &outfile)
Definition: ExpList.cpp:1863
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:306
std::pair< IndexType, IndexType > CoordType
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2675
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2088
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
std::shared_ptr< LibUtilities::Comm > GetComm()
Returns the comm object.
Definition: ExpList.h:980
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:2208
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2910
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
Definition: ExpList.cpp:2648
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:2333
std::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:1829
Collections::CollectionVector m_collections
Definition: ExpList.h:1092
virtual void v_ReadGlobalOptimizationParameters()
Definition: ExpList.cpp:3207
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Definition: ExpList.cpp:2755
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const FlagList &flags, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
Definition: ExpList.cpp:2701
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2918