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  return GetExp(GetExpIndex(gloCoord));
1434  }
1435 
1436 
1437  /**
1438  * @todo need a smarter search here that first just looks at bounding
1439  * vertices - suggest first seeing if point is within 10% of
1440  * region defined by vertices. The do point search.
1441  */
1443  const Array<OneD, const NekDouble> &gloCoord,
1444  NekDouble tol,
1445  bool returnNearestElmt)
1446  {
1447  Array<OneD, NekDouble> Lcoords(gloCoord.num_elements());
1448 
1449  return GetExpIndex(gloCoord,Lcoords,tol,returnNearestElmt);
1450  }
1451 
1452 
1454  const Array<OneD, const NekDouble> &gloCoords,
1455  Array<OneD, NekDouble> &locCoords,
1456  NekDouble tol,
1457  bool returnNearestElmt)
1458  {
1459  if (GetNumElmts() == 0)
1460  {
1461  return -1;
1462  }
1463 
1464  if (m_elmtToExpId.size() == 0)
1465  {
1466  // Loop in reverse order so that in case where using a
1467  // Homogeneous expansion it sets geometry ids to first part of
1468  // m_exp list. Otherwise will set to second (complex) expansion
1469  for(int i = (*m_exp).size()-1; i >= 0; --i)
1470  {
1471  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1472  }
1473  }
1474 
1475  NekDouble x = (gloCoords.num_elements() > 0 ? gloCoords[0] : 0.0);
1476  NekDouble y = (gloCoords.num_elements() > 1 ? gloCoords[1] : 0.0);
1477  NekDouble z = (gloCoords.num_elements() > 2 ? gloCoords[2] : 0.0);
1480  GetExp(0)->GetCoordim(), -1, x, y, z);
1481 
1482  // Get the list of elements whose bounding box contains the desired
1483  // point.
1484  std::vector<int> elmts = m_graph->GetElementsContainingPoint(p);
1485 
1486  NekDouble nearpt = 1e6;
1487  NekDouble nearpt_min = 1e6;
1488  int min_id = 0;
1489  Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
1490 
1491  // Check each element in turn to see if point lies within it.
1492  for (int i = 0; i < elmts.size(); ++i)
1493  {
1494  if ((*m_exp)[m_elmtToExpId[elmts[i]]]->
1495  GetGeom()->ContainsPoint(gloCoords,
1496  locCoords,
1497  tol, nearpt))
1498  {
1499  return m_elmtToExpId[elmts[i]];
1500  }
1501  else
1502  {
1503  // If it does not lie within, keep track of which element
1504  // is nearest.
1505  if(nearpt < nearpt_min)
1506  {
1507  min_id = m_elmtToExpId[elmts[i]];
1508  nearpt_min = nearpt;
1509  Vmath::Vcopy(locCoords.num_elements(),locCoords, 1,
1510  savLocCoords, 1);
1511  }
1512  }
1513  }
1514 
1515  // If the calling function is with just the nearest element, return
1516  // that. Otherwise return -1 to indicate no matching elemenet found.
1517  if(returnNearestElmt)
1518  {
1519 
1520  std::string msg = "Failed to find point within element to "
1521  "tolerance of "
1522  + boost::lexical_cast<std::string>(tol)
1523  + " using local point ("
1524  + boost::lexical_cast<std::string>(locCoords[0]) +","
1525  + boost::lexical_cast<std::string>(locCoords[1]) +","
1526  + boost::lexical_cast<std::string>(locCoords[1])
1527  + ") in element: "
1528  + boost::lexical_cast<std::string>(min_id);
1529  WARNINGL1(false,msg.c_str());
1530 
1531  Vmath::Vcopy(locCoords.num_elements(),savLocCoords, 1,
1532  locCoords, 1);
1533  return min_id;
1534  }
1535  else
1536  {
1537  return -1;
1538  }
1539  }
1540 
1541  /**
1542  * Given some coordinates, output the expansion field value at that
1543  * point
1544  */
1546  const Array<OneD, const NekDouble> &coords,
1547  const Array<OneD, const NekDouble> &phys)
1548  {
1549  int dim = GetCoordim(0);
1550  ASSERTL0(dim == coords.num_elements(),
1551  "Invalid coordinate dimension.");
1552 
1553  // Grab the element index corresponding to coords.
1554  Array<OneD, NekDouble> xi(dim);
1555  int elmtIdx = GetExpIndex(coords, xi);
1556  ASSERTL0(elmtIdx > 0, "Unable to find element containing point.");
1557 
1558  // Grab that element's physical storage.
1559  Array<OneD, NekDouble> elmtPhys = phys + m_phys_offset[elmtIdx];
1560 
1561  // Evaluate the element at the appropriate point.
1562  return (*m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
1563  }
1564 
1565  /**
1566  * Configures geometric info, such as tangent direction, on each
1567  * expansion.
1568  * @param graph2D Mesh
1569  */
1571  {
1572 
1573  }
1574 
1575  /**
1576  * @brief Reset geometry information, metrics, matrix managers and
1577  * geometry information.
1578  *
1579  * This routine clears all matrix managers and resets all geometry
1580  * information, which allows the geometry information to be dynamically
1581  * updated as the solver is run.
1582  */
1584  {
1585  // Reset matrix managers.
1587  DNekScalMat, LocalRegions::MatrixKey::opLess>::ClearManager();
1588  LibUtilities::NekManager<LocalRegions::MatrixKey,
1590 
1591  // Loop over all elements and reset geometry information.
1592  for (int i = 0; i < m_exp->size(); ++i)
1593  {
1594  (*m_exp)[i]->GetGeom()->Reset(m_graph->GetCurvedEdges(),
1595  m_graph->GetCurvedFaces());
1596  }
1597 
1598  // Loop over all elements and rebuild geometric factors.
1599  for (int i = 0; i < m_exp->size(); ++i)
1600  {
1601  (*m_exp)[i]->Reset();
1602  }
1603  }
1604 
1605  /**
1606  * Write Tecplot Files Header
1607  * @param outfile Output file name.
1608  * @param var variables names
1609  */
1610  void ExpList::v_WriteTecplotHeader(std::ostream &outfile,
1611  std::string var)
1612  {
1613  if (GetNumElmts() == 0)
1614  {
1615  return;
1616  }
1617 
1618  int coordim = GetExp(0)->GetCoordim();
1619  char vars[3] = { 'x', 'y', 'z' };
1620 
1621  if (m_expType == e3DH1D)
1622  {
1623  coordim += 1;
1624  }
1625  else if (m_expType == e3DH2D)
1626  {
1627  coordim += 2;
1628  }
1629 
1630  outfile << "Variables = x";
1631  for (int i = 1; i < coordim; ++i)
1632  {
1633  outfile << ", " << vars[i];
1634  }
1635 
1636  if (var.size() > 0)
1637  {
1638  outfile << ", " << var;
1639  }
1640 
1641  outfile << std::endl << std::endl;
1642  }
1643 
1644  /**
1645  * Write Tecplot Files Zone
1646  * @param outfile Output file name.
1647  * @param expansion Expansion that is considered
1648  */
1649  void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
1650  {
1651  int i, j;
1652  int coordim = GetCoordim(0);
1653  int nPoints = GetTotPoints();
1654  int nBases = (*m_exp)[0]->GetNumBases();
1655  int numBlocks = 0;
1656 
1658 
1659  if (expansion == -1)
1660  {
1661  nPoints = GetTotPoints();
1662 
1663  coords[0] = Array<OneD, NekDouble>(nPoints);
1664  coords[1] = Array<OneD, NekDouble>(nPoints);
1665  coords[2] = Array<OneD, NekDouble>(nPoints);
1666 
1667  GetCoords(coords[0], coords[1], coords[2]);
1668 
1669  for (i = 0; i < m_exp->size(); ++i)
1670  {
1671  int numInt = 1;
1672 
1673  for (j = 0; j < nBases; ++j)
1674  {
1675  numInt *= (*m_exp)[i]->GetNumPoints(j)-1;
1676  }
1677 
1678  numBlocks += numInt;
1679  }
1680  }
1681  else
1682  {
1683  nPoints = (*m_exp)[expansion]->GetTotPoints();
1684 
1685  coords[0] = Array<OneD, NekDouble>(nPoints);
1686  coords[1] = Array<OneD, NekDouble>(nPoints);
1687  coords[2] = Array<OneD, NekDouble>(nPoints);
1688 
1689  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
1690 
1691  numBlocks = 1;
1692  for (j = 0; j < nBases; ++j)
1693  {
1694  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j)-1;
1695  }
1696  }
1697 
1698  if (m_expType == e3DH1D)
1699  {
1700  nBases += 1;
1701  coordim += 1;
1702  int nPlanes = GetZIDs().num_elements();
1703  NekDouble tmp = numBlocks * (nPlanes-1.0) / nPlanes;
1704  numBlocks = (int)tmp;
1705  }
1706  else if (m_expType == e3DH2D)
1707  {
1708  nBases += 2;
1709  coordim += 1;
1710  }
1711 
1712  outfile << "Zone, N=" << nPoints << ", E="
1713  << numBlocks << ", F=FEBlock" ;
1714 
1715  switch(nBases)
1716  {
1717  case 2:
1718  outfile << ", ET=QUADRILATERAL" << std::endl;
1719  break;
1720  case 3:
1721  outfile << ", ET=BRICK" << std::endl;
1722  break;
1723  default:
1724  ASSERTL0(false,"Not set up for this type of output");
1725  break;
1726  }
1727 
1728  // Write out coordinates
1729  for (j = 0; j < coordim; ++j)
1730  {
1731  for (i = 0; i < nPoints; ++i)
1732  {
1733  outfile << coords[j][i] << " ";
1734  if (i % 1000 == 0 && i)
1735  {
1736  outfile << std::endl;
1737  }
1738  }
1739  outfile << std::endl;
1740  }
1741  }
1742 
1743  void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile,
1744  int expansion)
1745  {
1746  int i,j,k,l;
1747  int nbase = (*m_exp)[0]->GetNumBases();
1748  int cnt = 0;
1749 
1750  std::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
1751 
1752  if (expansion != -1)
1753  {
1754  exp = std::shared_ptr<LocalRegions::ExpansionVector>(
1756  (*exp)[0] = (*m_exp)[expansion];
1757  }
1758 
1759  if (nbase == 2)
1760  {
1761  for(i = 0; i < (*exp).size(); ++i)
1762  {
1763  const int np0 = (*exp)[i]->GetNumPoints(0);
1764  const int np1 = (*exp)[i]->GetNumPoints(1);
1765 
1766  for(j = 1; j < np1; ++j)
1767  {
1768  for(k = 1; k < np0; ++k)
1769  {
1770  outfile << cnt + (j-1)*np0 + k << " ";
1771  outfile << cnt + (j-1)*np0 + k+1 << " ";
1772  outfile << cnt + j *np0 + k+1 << " ";
1773  outfile << cnt + j *np0 + k << endl;
1774  }
1775  }
1776 
1777  cnt += np0*np1;
1778  }
1779  }
1780  else if (nbase == 3)
1781  {
1782  for(i = 0; i < (*exp).size(); ++i)
1783  {
1784  const int np0 = (*exp)[i]->GetNumPoints(0);
1785  const int np1 = (*exp)[i]->GetNumPoints(1);
1786  const int np2 = (*exp)[i]->GetNumPoints(2);
1787  const int np01 = np0*np1;
1788 
1789  for(j = 1; j < np2; ++j)
1790  {
1791  for(k = 1; k < np1; ++k)
1792  {
1793  for(l = 1; l < np0; ++l)
1794  {
1795  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l << " ";
1796  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l+1 << " ";
1797  outfile << cnt + (j-1)*np01 + k *np0 + l+1 << " ";
1798  outfile << cnt + (j-1)*np01 + k *np0 + l << " ";
1799  outfile << cnt + j *np01 + (k-1)*np0 + l << " ";
1800  outfile << cnt + j *np01 + (k-1)*np0 + l+1 << " ";
1801  outfile << cnt + j *np01 + k *np0 + l+1 << " ";
1802  outfile << cnt + j *np01 + k *np0 + l << endl;
1803  }
1804  }
1805  }
1806  cnt += np0*np1*np2;
1807  }
1808  }
1809  else
1810  {
1811  ASSERTL0(false,"Not set up for this dimension");
1812  }
1813  }
1814 
1815  /**
1816  * Write Tecplot Files Field
1817  * @param outfile Output file name.
1818  * @param expansion Expansion that is considered
1819  */
1820  void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
1821  {
1822  if (expansion == -1)
1823  {
1824  int totpoints = GetTotPoints();
1825  if(m_physState == false)
1826  {
1828  }
1829 
1830  for(int i = 0; i < totpoints; ++i)
1831  {
1832  outfile << m_phys[i] << " ";
1833  if(i % 1000 == 0 && i)
1834  {
1835  outfile << std::endl;
1836  }
1837  }
1838  outfile << std::endl;
1839 
1840  }
1841  else
1842  {
1843  int nPoints = (*m_exp)[expansion]->GetTotPoints();
1844 
1845  for (int i = 0; i < nPoints; ++i)
1846  {
1847  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
1848  }
1849 
1850  outfile << std::endl;
1851  }
1852  }
1853 
1854  void ExpList::WriteVtkHeader(std::ostream &outfile)
1855  {
1856  outfile << "<?xml version=\"1.0\"?>" << endl;
1857  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1858  << "byte_order=\"LittleEndian\">" << endl;
1859  outfile << " <UnstructuredGrid>" << endl;
1860  }
1861 
1862  void ExpList::WriteVtkFooter(std::ostream &outfile)
1863  {
1864  outfile << " </UnstructuredGrid>" << endl;
1865  outfile << "</VTKFile>" << endl;
1866  }
1867 
1868  void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
1869  {
1870  boost::ignore_unused(outfile, expansion, istrip);
1872  "Routine not implemented for this expansion.");
1873  }
1874 
1875  void ExpList::WriteVtkPieceFooter(std::ostream &outfile, int expansion)
1876  {
1877  boost::ignore_unused(expansion);
1878  outfile << " </PointData>" << endl;
1879  outfile << " </Piece>" << endl;
1880  }
1881 
1882  void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
1883  std::string var)
1884  {
1885  int i;
1886  int nq = (*m_exp)[expansion]->GetTotPoints();
1887 
1888  // printing the fields of that zone
1889  outfile << " <DataArray type=\"Float64\" Name=\""
1890  << var << "\">" << endl;
1891  outfile << " ";
1892  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
1893  for(i = 0; i < nq; ++i)
1894  {
1895  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
1896  }
1897  outfile << endl;
1898  outfile << " </DataArray>" << endl;
1899  }
1900 
1901  /**
1902  * Given a spectral/hp approximation
1903  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1904  * (which should be contained in #m_phys), this function calculates the
1905  * \f$L_\infty\f$ error of this approximation with respect to an exact
1906  * solution. The local distribution of the quadrature points allows an
1907  * elemental evaluation of this operation through the functions
1908  * StdRegions#StdExpansion#Linf.
1909  *
1910  * The exact solution, also evaluated at the quadrature
1911  * points, should be contained in the variable #m_phys of
1912  * the ExpList object \a Sol.
1913  *
1914  * @param soln A 1D array, containing the discrete
1915  * evaluation of the exact solution at the
1916  * quadrature points in its array #m_phys.
1917  * @return The \f$L_\infty\f$ error of the approximation.
1918  */
1920  const Array<OneD, const NekDouble> &inarray,
1921  const Array<OneD, const NekDouble> &soln)
1922  {
1923  NekDouble err = 0.0;
1924 
1925  if (soln == NullNekDouble1DArray)
1926  {
1927  err = Vmath::Vmax(m_npoints, inarray, 1);
1928  }
1929  else
1930  {
1931  for (int i = 0; i < m_npoints; ++i)
1932  {
1933  err = max(err, abs(inarray[i] - soln[i]));
1934  }
1935  }
1936 
1937  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
1938 
1939  return err;
1940  }
1941 
1942  /**
1943  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
1944  * evaluated at the quadrature points (which should be contained in
1945  * #m_phys), this function calculates the \f$L_2\f$ error of this
1946  * approximation with respect to an exact solution. The local
1947  * distribution of the quadrature points allows an elemental evaluation
1948  * of this operation through the functions StdRegions#StdExpansion#L2.
1949  *
1950  * The exact solution, also evaluated at the quadrature points, should
1951  * be contained in the variable #m_phys of the ExpList object \a Sol.
1952  *
1953  * @param Sol An ExpList, containing the discrete
1954  * evaluation of the exact solution at the
1955  * quadrature points in its array #m_phys.
1956  * @return The \f$L_2\f$ error of the approximation.
1957  */
1959  const Array<OneD, const NekDouble> &inarray,
1960  const Array<OneD, const NekDouble> &soln)
1961  {
1962  NekDouble err = 0.0, errl2;
1963  int i;
1964 
1965  if (soln == NullNekDouble1DArray)
1966  {
1967  for (i = 0; i < (*m_exp).size(); ++i)
1968  {
1969  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
1970  err += errl2*errl2;
1971  }
1972  }
1973  else
1974  {
1975  for (i = 0; i < (*m_exp).size(); ++i)
1976  {
1977  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
1978  soln + m_phys_offset[i]);
1979  err += errl2*errl2;
1980  }
1981  }
1982 
1983  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1984 
1985  return sqrt(err);
1986  }
1987 
1989  {
1990  NekDouble err = 0.0;
1991  int i = 0;
1992 
1993  for (i = 0; i < (*m_exp).size(); ++i)
1994  {
1995  err += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
1996  }
1997  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1998 
1999  return err;
2000  }
2001 
2003  {
2004  NekDouble flux = 0.0;
2005  int i = 0;
2006  int j;
2007 
2008  for (i = 0; i < (*m_exp).size(); ++i)
2009  {
2010  Array<OneD, Array<OneD, NekDouble> > tmp(inarray.num_elements());
2011  for (j = 0; j < inarray.num_elements(); ++j)
2012  {
2013  tmp[j] = Array<OneD, NekDouble>(inarray[j] + m_phys_offset[i]);
2014  }
2015  flux += (*m_exp)[i]->VectorFlux(tmp);
2016  }
2017 
2018  return flux;
2019  }
2020 
2022  {
2024  "This method is not defined or valid for this class type");
2025  Array<OneD, NekDouble> NoEnergy(1,0.0);
2026  return NoEnergy;
2027  }
2028 
2030  {
2032  "This method is not defined or valid for this class type");
2034  return trans;
2035  }
2036 
2038  {
2040  "This method is not defined or valid for this class type");
2041  NekDouble len = 0.0;
2042  return len;
2043  }
2044 
2046  {
2047  boost::ignore_unused(lhom);
2049  "This method is not defined or valid for this class type");
2050  }
2051 
2053  {
2055  "This method is not defined or valid for this class type");
2056  Array<OneD, unsigned int> NoModes(1);
2057  return NoModes;
2058  }
2059 
2061  {
2063  "This method is not defined or valid for this class type");
2064  Array<OneD, unsigned int> NoModes(1);
2065  return NoModes;
2066  }
2067 
2068 
2070  const NekDouble scale,
2071  const Array<OneD, NekDouble> &inarray,
2072  Array<OneD, NekDouble> &outarray)
2073  {
2074  boost::ignore_unused(scale, inarray, outarray);
2076  "This method is not defined or valid for this class type");
2077  }
2078 
2080  const NekDouble scale,
2081  const Array<OneD, NekDouble> &inarray,
2082  Array<OneD, NekDouble> &outarray)
2083  {
2084  boost::ignore_unused(scale, inarray, outarray);
2086  "This method is not defined or valid for this class type");
2087  }
2088 
2090  {
2092  "This method is not defined or valid for this class type");
2093  }
2094 
2096  const std::string &fileName,
2098  const std::string &varName,
2099  const std::shared_ptr<ExpList> locExpList)
2100  {
2101  string varString = fileName.substr(0, fileName.find_last_of("."));
2102  int j, k, len = varString.length();
2103  varString = varString.substr(len-1, len);
2104 
2105  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
2106  std::vector<std::vector<NekDouble> > FieldData;
2107 
2108  std::string ft = LibUtilities::FieldIO::GetFileType(fileName, comm);
2110  .CreateInstance(ft, comm, m_session->GetSharedFilesystem());
2111 
2112  f->Import(fileName, FieldDef, FieldData);
2113 
2114  bool found = false;
2115  for (j = 0; j < FieldDef.size(); ++j)
2116  {
2117  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
2118  {
2119  if (FieldDef[j]->m_fields[k] == varName)
2120  {
2121  // Copy FieldData into locExpList
2122  locExpList->ExtractDataToCoeffs(
2123  FieldDef[j], FieldData[j],
2124  FieldDef[j]->m_fields[k],
2125  locExpList->UpdateCoeffs());
2126  found = true;
2127  }
2128  }
2129  }
2130 
2131  ASSERTL0(found, "Could not find variable '"+varName+
2132  "' in file boundary condition "+fileName);
2133  locExpList->BwdTrans_IterPerExp(
2134  locExpList->GetCoeffs(),
2135  locExpList->UpdatePhys());
2136  }
2137 
2138  /**
2139  * Given a spectral/hp approximation
2140  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
2141  * (which should be contained in #m_phys), this function calculates the
2142  * \f$H^1_2\f$ error of this approximation with respect to an exact
2143  * solution. The local distribution of the quadrature points allows an
2144  * elemental evaluation of this operation through the functions
2145  * StdRegions#StdExpansion#H1.
2146  *
2147  * The exact solution, also evaluated at the quadrature points, should
2148  * be contained in the variable #m_phys of the ExpList object \a Sol.
2149  *
2150  * @param soln An 1D array, containing the discrete evaluation
2151  * of the exact solution at the quadrature points.
2152  *
2153  * @return The \f$H^1_2\f$ error of the approximation.
2154  */
2156  const Array<OneD, const NekDouble> &inarray,
2157  const Array<OneD, const NekDouble> &soln)
2158  {
2159  NekDouble err = 0.0, errh1;
2160  int i;
2161 
2162  for (i = 0; i < (*m_exp).size(); ++i)
2163  {
2164  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
2165  soln + m_phys_offset[i]);
2166  err += errh1*errh1;
2167  }
2168 
2169  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
2170 
2171  return sqrt(err);
2172  }
2173 
2174  void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
2175  int NumHomoDir,
2177  std::vector<NekDouble> &HomoLen,
2178  bool homoStrips,
2179  std::vector<unsigned int> &HomoSIDs,
2180  std::vector<unsigned int> &HomoZIDs,
2181  std::vector<unsigned int> &HomoYIDs)
2182  {
2183  int startenum = (int) LibUtilities::eSegment;
2184  int endenum = (int) LibUtilities::eHexahedron;
2185  int s = 0;
2187 
2188  ASSERTL1(NumHomoDir == HomoBasis.num_elements(),"Homogeneous basis is not the same length as NumHomoDir");
2189  ASSERTL1(NumHomoDir == HomoLen.size(),"Homogeneous length vector is not the same length as NumHomDir");
2190 
2191  // count number of shapes
2192  switch((*m_exp)[0]->GetShapeDimension())
2193  {
2194  case 1:
2195  startenum = (int) LibUtilities::eSegment;
2196  endenum = (int) LibUtilities::eSegment;
2197  break;
2198  case 2:
2199  startenum = (int) LibUtilities::eTriangle;
2200  endenum = (int) LibUtilities::eQuadrilateral;
2201  break;
2202  case 3:
2203  startenum = (int) LibUtilities::eTetrahedron;
2204  endenum = (int) LibUtilities::eHexahedron;
2205  break;
2206  }
2207 
2208  for(s = startenum; s <= endenum; ++s)
2209  {
2210  std::vector<unsigned int> elementIDs;
2211  std::vector<LibUtilities::BasisType> basis;
2212  std::vector<unsigned int> numModes;
2213  std::vector<std::string> fields;
2214 
2215  bool first = true;
2216  bool UniOrder = true;
2217  int n;
2218 
2219  shape = (LibUtilities::ShapeType) s;
2220 
2221  for(int i = 0; i < (*m_exp).size(); ++i)
2222  {
2223  if((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
2224  {
2225  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
2226  if(first)
2227  {
2228  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
2229  {
2230  basis.push_back((*m_exp)[i]->GetBasis(j)->GetBasisType());
2231  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
2232  }
2233 
2234  // add homogeneous direction details if defined
2235  for(n = 0 ; n < NumHomoDir; ++n)
2236  {
2237  basis.push_back(HomoBasis[n]->GetBasisType());
2238  numModes.push_back(HomoBasis[n]->GetNumModes());
2239  }
2240 
2241  first = false;
2242  }
2243  else
2244  {
2245  ASSERTL0((*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],"Routine is not set up for multiple bases definitions");
2246 
2247  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
2248  {
2249  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
2250  if(numModes[j] != (*m_exp)[i]->GetBasis(j)->GetNumModes())
2251  {
2252  UniOrder = false;
2253  }
2254  }
2255  // add homogeneous direction details if defined
2256  for(n = 0 ; n < NumHomoDir; ++n)
2257  {
2258  numModes.push_back(HomoBasis[n]->GetNumModes());
2259  }
2260  }
2261  }
2262  }
2263 
2264 
2265  if(elementIDs.size() > 0)
2266  {
2269  AllocateSharedPtr(shape, elementIDs, basis,
2270  UniOrder, numModes,fields,
2271  NumHomoDir, HomoLen, homoStrips,
2272  HomoSIDs, HomoZIDs, HomoYIDs);
2273  fielddef.push_back(fdef);
2274  }
2275  }
2276  }
2277 
2278 
2279  //
2280  // Virtual functions
2281  //
2282  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::v_GetFieldDefinitions()
2283  {
2284  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
2285  v_GetFieldDefinitions(returnval);
2286  return returnval;
2287  }
2288 
2289  void ExpList::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
2290  {
2291  GeneralGetFieldDefinitions(fielddef);
2292  }
2293 
2294  //Append the element data listed in elements
2295  //fielddef->m_ElementIDs onto fielddata
2296  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
2297  {
2298  v_AppendFieldData(fielddef,fielddata,m_coeffs);
2299  }
2300 
2301  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
2302  {
2303  int i;
2304  // Determine mapping from element ids to location in
2305  // expansion list
2306  // Determine mapping from element ids to location in
2307  // expansion list
2308  map<int, int> ElmtID_to_ExpID;
2309  for(i = 0; i < (*m_exp).size(); ++i)
2310  {
2311  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2312  }
2313 
2314  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
2315  {
2316  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
2317  int datalen = (*m_exp)[eid]->GetNcoeffs();
2318  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]],&coeffs[m_coeff_offset[eid]]+datalen);
2319  }
2320 
2321  }
2322 
2323  /// Extract the data in fielddata into the coeffs
2326  std::vector<NekDouble> &fielddata,
2327  std::string &field,
2328  Array<OneD, NekDouble> &coeffs)
2329  {
2330  v_ExtractDataToCoeffs(fielddef,fielddata,field,coeffs);
2331  }
2332 
2333  void ExpList::ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2334  {
2335  v_ExtractCoeffsToCoeffs(fromExpList,fromCoeffs,toCoeffs);
2336  }
2337 
2338  /**
2339  * @brief Extract data from raw field data into expansion list.
2340  *
2341  * @param fielddef Field definitions.
2342  * @param fielddata Data for associated field.
2343  * @param field Field variable name.
2344  * @param coeffs Resulting coefficient array.
2345  */
2348  std::vector<NekDouble> &fielddata,
2349  std::string &field,
2350  Array<OneD, NekDouble> &coeffs)
2351  {
2352  int i, expId;
2353  int offset = 0;
2354  int modes_offset = 0;
2355  int datalen = fielddata.size()/fielddef->m_fields.size();
2356 
2357  // Find data location according to field definition
2358  for(i = 0; i < fielddef->m_fields.size(); ++i)
2359  {
2360  if(fielddef->m_fields[i] == field)
2361  {
2362  break;
2363  }
2364  offset += datalen;
2365  }
2366 
2367  ASSERTL0(i != fielddef->m_fields.size(),
2368  "Field (" + field + ") not found in file.");
2369 
2370  if (m_elmtToExpId.size() == 0)
2371  {
2372  // Loop in reverse order so that in case where using a
2373  // Homogeneous expansion it sets geometry ids to first part of
2374  // m_exp list. Otherwise will set to second (complex) expansion
2375  for(i = (*m_exp).size()-1; i >= 0; --i)
2376  {
2377  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2378  }
2379  }
2380 
2381  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
2382  {
2383  // Reset modes_offset in the case where all expansions of
2384  // the same order.
2385  if (fielddef->m_uniOrder == true)
2386  {
2387  modes_offset = 0;
2388  }
2389 
2391  fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
2392 
2393  const int elmtId = fielddef->m_elementIDs[i];
2394  auto eIt = m_elmtToExpId.find(elmtId);
2395 
2396  if (eIt == m_elmtToExpId.end())
2397  {
2398  offset += datalen;
2399  modes_offset += (*m_exp)[0]->GetNumBases();
2400  continue;
2401  }
2402 
2403  expId = eIt->second;
2404 
2405  bool sameBasis = true;
2406  for (int j = 0; j < fielddef->m_basis.size(); ++j)
2407  {
2408  if (fielddef->m_basis[j] != (*m_exp)[expId]->GetBasisType(j))
2409  {
2410  sameBasis = false;
2411  break;
2412  }
2413  }
2414 
2415  if (datalen == (*m_exp)[expId]->GetNcoeffs() && sameBasis)
2416  {
2417  Vmath::Vcopy(datalen, &fielddata[offset], 1,
2418  &coeffs[m_coeff_offset[expId]], 1);
2419  }
2420  else
2421  {
2422  (*m_exp)[expId]->ExtractDataToCoeffs(
2423  &fielddata[offset], fielddef->m_numModes,
2424  modes_offset, &coeffs[m_coeff_offset[expId]],
2425  fielddef->m_basis);
2426  }
2427 
2428  offset += datalen;
2429  modes_offset += (*m_exp)[0]->GetNumBases();
2430  }
2431 
2432  return;
2433  }
2434 
2435  void ExpList::v_ExtractCoeffsToCoeffs(const std::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2436  {
2437  int i;
2438  int offset = 0;
2439 
2440  for(i = 0; i < (*m_exp).size(); ++i)
2441  {
2442  std::vector<unsigned int> nummodes;
2443  vector<LibUtilities::BasisType> basisTypes;
2444  for(int j= 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
2445  {
2446  nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
2447  basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
2448  }
2449 
2450  (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
2451  &toCoeffs[m_coeff_offset[i]],
2452  basisTypes);
2453 
2454  offset += fromExpList->GetExp(i)->GetNcoeffs();
2455  }
2456  }
2457 
2458 
2460  const SpatialDomains::GeomMMF MMFdir,
2461  const Array<OneD, const NekDouble> &CircCentre,
2462  Array<OneD, Array<OneD, NekDouble> > &outarray)
2463  {
2464  int npts;
2465 
2466  int MFdim = 3;
2467  int nq = outarray[0].num_elements()/MFdim;
2468 
2469  // Assume whole array is of same coordinate dimension
2470  int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
2471 
2472  Array<OneD, Array<OneD, NekDouble> > MFloc(MFdim*coordim);
2473  // Process each expansion.
2474  for(int i = 0; i < m_exp->size(); ++i)
2475  {
2476  npts = (*m_exp)[i]->GetTotPoints();
2477 
2478  for (int j=0; j< MFdim*coordim; ++j)
2479  {
2480  MFloc[j] = Array<OneD, NekDouble>(npts, 0.0);
2481  }
2482 
2483  // MF from LOCALREGIONS
2484  (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
2485  (*m_exp)[i]->GetPointsKeys(),
2486  MMFdir,
2487  CircCentre,
2488  MFloc );
2489 
2490  // Get the physical data offset for this expansion.
2491  for (int j = 0; j < MFdim; ++j)
2492  {
2493  for (int k = 0; k < coordim; ++k)
2494  {
2495  Vmath::Vcopy(npts,
2496  &MFloc[j*coordim+k][0], 1,
2497  &outarray[j][k*nq+m_phys_offset[i]], 1);
2498  }
2499  }
2500  }
2501 
2502  }
2503 
2504 
2505 
2506  /**
2507  * @brief Generate vector v such that v[i] = scalar1 if i is in the
2508  * element < ElementID. Otherwise, v[i] = scalar2.
2509  *
2510  */
2512  const int ElementID,
2513  const NekDouble scalar1,
2514  const NekDouble scalar2,
2515  Array<OneD, NekDouble> &outarray)
2516  {
2517  int npoints_e;
2518  NekDouble coeff;
2519 
2520  Array<OneD, NekDouble> outarray_e;
2521 
2522  for(int i = 0 ; i < (*m_exp).size(); ++i)
2523  {
2524  npoints_e = (*m_exp)[i]->GetTotPoints();
2525 
2526  if(i <= ElementID)
2527  {
2528  coeff = scalar1;
2529  }
2530  else
2531  {
2532  coeff = scalar2;
2533  }
2534 
2535  outarray_e = Array<OneD, NekDouble>(npoints_e, coeff);
2536  Vmath::Vcopy(npoints_e, &outarray_e[0], 1,
2537  &outarray[m_phys_offset[i]], 1);
2538  }
2539  }
2540 
2543  {
2544  ASSERTL0(false,
2545  "This method is not defined or valid for this class type");
2547  return result;
2548  }
2549 
2550  std::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
2551  {
2552  boost::ignore_unused(i);
2554  "This method is not defined or valid for this class type");
2555  static std::shared_ptr<ExpList> result;
2556  return result;
2557  }
2558 
2560  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
2561  const Array<OneD, const NekDouble> &Fwd,
2562  const Array<OneD, const NekDouble> &Bwd,
2564  {
2565  boost::ignore_unused(Vec, Fwd, Bwd, Upwind);
2567  "This method is not defined or valid for this class type");
2568  }
2569 
2571  const Array<OneD, const NekDouble> &Vn,
2572  const Array<OneD, const NekDouble> &Fwd,
2573  const Array<OneD, const NekDouble> &Bwd,
2575  {
2576  boost::ignore_unused(Vn, Fwd, Bwd, Upwind);
2578  "This method is not defined or valid for this class type");
2579  }
2580 
2581  std::shared_ptr<ExpList> &ExpList::v_GetTrace()
2582  {
2584  "This method is not defined or valid for this class type");
2585  static std::shared_ptr<ExpList> returnVal;
2586  return returnVal;
2587  }
2588 
2589  std::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
2590  {
2592  "This method is not defined or valid for this class type");
2593  static std::shared_ptr<AssemblyMapDG> result;
2594  return result;
2595  }
2596 
2598  {
2599  return GetTraceMap()->GetBndCondTraceToGlobalTraceMap();
2600  }
2601 
2603  Array<OneD, Array<OneD, NekDouble> > &normals)
2604  {
2605  boost::ignore_unused(normals);
2607  "This method is not defined or valid for this class type");
2608  }
2609 
2611  const Array<OneD, const NekDouble> &Fx,
2612  const Array<OneD, const NekDouble> &Fy,
2613  Array<OneD, NekDouble> &outarray)
2614  {
2615  boost::ignore_unused(Fx, Fy, outarray);
2617  "This method is not defined or valid for this class type");
2618  }
2619 
2621  const Array<OneD, const NekDouble> &Fn,
2622  Array<OneD, NekDouble> &outarray)
2623  {
2624  boost::ignore_unused(Fn, outarray);
2626  "This method is not defined or valid for this class type");
2627  }
2628 
2630  const Array<OneD, const NekDouble> &Fwd,
2631  const Array<OneD, const NekDouble> &Bwd,
2632  Array<OneD, NekDouble> &outarray)
2633  {
2634  boost::ignore_unused(Fwd, Bwd, outarray);
2636  "This method is not defined or valid for this class type");
2637  }
2638 
2640  Array<OneD,NekDouble> &Bwd)
2641  {
2642  boost::ignore_unused(Fwd, Bwd);
2644  "This method is not defined or valid for this class type");
2645  }
2646 
2648  const Array<OneD,const NekDouble> &field,
2649  Array<OneD,NekDouble> &Fwd,
2650  Array<OneD,NekDouble> &Bwd)
2651  {
2652  boost::ignore_unused(field, Fwd, Bwd);
2654  "This method is not defined or valid for this class type");
2655  }
2656 
2657  const vector<bool> &ExpList::v_GetLeftAdjacentFaces(void) const
2658  {
2660  "This method is not defined or valid for this class type");
2661  static vector<bool> tmp;
2662  return tmp;
2663  }
2664 
2665 
2667  {
2668  boost::ignore_unused(outarray);
2670  "This method is not defined or valid for this class type");
2671  }
2672 
2674  const Array<OneD, const NekDouble> &inarray,
2675  Array<OneD,NekDouble> &outarray)
2676  {
2677  boost::ignore_unused(inarray, 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  CoeffState coeffstate)
2686  {
2687  boost::ignore_unused(inarray, outarray, coeffstate);
2689  "This method is not defined or valid for this class type");
2690  }
2691 
2693  const Array<OneD, const NekDouble> &inarray,
2694  Array<OneD, NekDouble> &outarray,
2695  const FlagList &flags,
2696  const StdRegions::ConstFactorMap &factors,
2697  const StdRegions::VarCoeffMap &varcoeff,
2698  const MultiRegions::VarFactorsMap &varfactors,
2699  const Array<OneD, const NekDouble> &dirForcing,
2700  const bool PhysSpaceForcing)
2701  {
2702  boost::ignore_unused(inarray, outarray, flags, factors, varcoeff,
2703  varfactors, dirForcing, PhysSpaceForcing);
2704  NEKERROR(ErrorUtil::efatal, "HelmSolve not implemented.");
2705  }
2706 
2708  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2709  const Array<OneD, const NekDouble> &inarray,
2710  Array<OneD, NekDouble> &outarray,
2711  const NekDouble lambda,
2712  CoeffState coeffstate,
2713  const Array<OneD, const NekDouble>& dirForcing)
2714  {
2715  boost::ignore_unused(velocity, inarray, outarray, lambda,
2716  coeffstate, dirForcing);
2718  "This method is not defined or valid for this class type");
2719  }
2720 
2722  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2723  const Array<OneD, const NekDouble> &inarray,
2724  Array<OneD, NekDouble> &outarray,
2725  const NekDouble lambda,
2726  CoeffState coeffstate,
2727  const Array<OneD, const NekDouble>& dirForcing)
2728  {
2729  boost::ignore_unused(velocity, inarray, outarray, lambda,
2730  coeffstate, dirForcing);
2732  "This method is not defined or valid for this class type");
2733  }
2734 
2736  Array<OneD, NekDouble> &outarray,
2737  CoeffState coeffstate,
2738  bool Shuff,
2739  bool UnShuff)
2740  {
2741  boost::ignore_unused(inarray, outarray, coeffstate, Shuff, UnShuff);
2743  "This method is not defined or valid for this class type");
2744  }
2745 
2747  Array<OneD, NekDouble> &outarray,
2748  CoeffState coeffstate,
2749  bool Shuff,
2750  bool UnShuff)
2751  {
2752  boost::ignore_unused(inarray, outarray, coeffstate, Shuff, UnShuff);
2754  "This method is not defined or valid for this class type");
2755  }
2756 
2758  const Array<OneD, NekDouble> &inarray2,
2759  Array<OneD, NekDouble> &outarray,
2760  CoeffState coeffstate)
2761  {
2762  boost::ignore_unused(inarray1, inarray2, outarray, coeffstate);
2764  "This method is not defined or valid for this class type");
2765  }
2766 
2768  const Array<OneD, Array<OneD, NekDouble> > &inarray1,
2769  const Array<OneD, Array<OneD, NekDouble> > &inarray2,
2770  Array<OneD, Array<OneD, NekDouble> > &outarray,
2771  CoeffState coeffstate)
2772  {
2773  boost::ignore_unused(inarray1, inarray2, outarray, coeffstate);
2775  "This method is not defined or valid for this class type");
2776  }
2777 
2779  const Array<OneD, NekDouble> &TotField,
2780  int BndID)
2781  {
2782  boost::ignore_unused(BndVals, TotField, BndID);
2784  "This method is not defined or valid for this class type");
2785  }
2786 
2789  Array<OneD, NekDouble> &outarray,
2790  int BndID)
2791  {
2792  boost::ignore_unused(V1, V2, outarray, BndID);
2794  "This method is not defined or valid for this class type");
2795  }
2796 
2799  Array<OneD, NekDouble> &outarray)
2800  {
2802  switch (GetCoordim(0))
2803  {
2804  case 1:
2805  {
2806  for(int i = 0; i < GetExpSize(); ++i)
2807  {
2808  (*m_exp)[i]->NormVectorIProductWRTBase(
2809  V[0] + GetPhys_Offset(i),
2810  tmp = outarray + GetCoeff_Offset(i));
2811  }
2812  }
2813  break;
2814  case 2:
2815  {
2816  for(int i = 0; i < GetExpSize(); ++i)
2817  {
2818  (*m_exp)[i]->NormVectorIProductWRTBase(
2819  V[0] + GetPhys_Offset(i),
2820  V[1] + GetPhys_Offset(i),
2821  tmp = outarray + GetCoeff_Offset(i));
2822  }
2823  }
2824  break;
2825  case 3:
2826  {
2827  for(int i = 0; i < GetExpSize(); ++i)
2828  {
2829  (*m_exp)[i]->NormVectorIProductWRTBase(
2830  V[0] + GetPhys_Offset(i),
2831  V[1] + GetPhys_Offset(i),
2832  V[2] + GetPhys_Offset(i),
2833  tmp = outarray + GetCoeff_Offset(i));
2834  }
2835  }
2836  break;
2837  default:
2838  ASSERTL0(false,"Dimension not supported");
2839  break;
2840  }
2841  }
2842 
2844  {
2845  boost::ignore_unused(outarray);
2847  "This method is not defined or valid for this class type");
2848  }
2849 
2850  /**
2851  */
2853  {
2855  "This method is not defined or valid for this class type");
2856  }
2857 
2858  /**
2859  */
2860  void ExpList::v_FillBndCondFromField(const int nreg)
2861  {
2862  boost::ignore_unused(nreg);
2864  "This method is not defined or valid for this class type");
2865  }
2866 
2867  void ExpList::v_LocalToGlobal(bool useComm)
2868  {
2869  boost::ignore_unused(useComm);
2871  "This method is not defined or valid for this class type");
2872  }
2873 
2874 
2876  Array<OneD,NekDouble> &outarray,
2877  bool useComm)
2878  {
2879  boost::ignore_unused(inarray, outarray, useComm);
2881  "This method is not defined or valid for this class type");
2882  }
2883 
2884 
2886  {
2888  "This method is not defined or valid for this class type");
2889  }
2890 
2891 
2893  Array<OneD,NekDouble> &outarray)
2894  {
2895  boost::ignore_unused(inarray, outarray);
2897  "This method is not defined or valid for this class type");
2898  }
2899 
2900 
2902  Array<OneD, NekDouble> &outarray,
2903  CoeffState coeffstate)
2904  {
2905  boost::ignore_unused(coeffstate);
2906  v_BwdTrans_IterPerExp(inarray,outarray);
2907  }
2908 
2910  Array<OneD, NekDouble> &outarray,
2911  CoeffState coeffstate)
2912  {
2913  boost::ignore_unused(coeffstate);
2914  v_FwdTrans_IterPerExp(inarray,outarray);
2915  }
2916 
2918  const Array<OneD, const NekDouble> &inarray,
2919  Array<OneD, NekDouble> &outarray,
2920  CoeffState coeffstate)
2921  {
2922  boost::ignore_unused(coeffstate);
2924  for (int i = 0; i < m_collections.size(); ++i)
2925  {
2926 
2928  inarray + m_coll_phys_offset[i],
2929  tmp = outarray + m_coll_coeff_offset[i]);
2930  }
2931  }
2932 
2934  const GlobalMatrixKey &gkey,
2935  const Array<OneD,const NekDouble> &inarray,
2936  Array<OneD, NekDouble> &outarray,
2937  CoeffState coeffstate)
2938  {
2939  boost::ignore_unused(coeffstate);
2940  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
2941  }
2942 
2943  /**
2944  * The operation is evaluated locally by the elemental
2945  * function StdRegions#StdExpansion#GetCoords.
2946  *
2947  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
2948  * will be stored in this array.
2949  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
2950  * will be stored in this array.
2951  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
2952  * will be stored in this array.
2953  */
2955  Array<OneD, NekDouble> &coord_1,
2956  Array<OneD, NekDouble> &coord_2)
2957  {
2958  if (GetNumElmts() == 0)
2959  {
2960  return;
2961  }
2962 
2963  int i;
2964  Array<OneD, NekDouble> e_coord_0;
2965  Array<OneD, NekDouble> e_coord_1;
2966  Array<OneD, NekDouble> e_coord_2;
2967 
2968  switch(GetExp(0)->GetCoordim())
2969  {
2970  case 1:
2971  for(i= 0; i < (*m_exp).size(); ++i)
2972  {
2973  e_coord_0 = coord_0 + m_phys_offset[i];
2974  (*m_exp)[i]->GetCoords(e_coord_0);
2975  }
2976  break;
2977  case 2:
2978  ASSERTL0(coord_1.num_elements() != 0,
2979  "output coord_1 is not defined");
2980 
2981  for(i= 0; i < (*m_exp).size(); ++i)
2982  {
2983  e_coord_0 = coord_0 + m_phys_offset[i];
2984  e_coord_1 = coord_1 + m_phys_offset[i];
2985  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1);
2986  }
2987  break;
2988  case 3:
2989  ASSERTL0(coord_1.num_elements() != 0,
2990  "output coord_1 is not defined");
2991  ASSERTL0(coord_2.num_elements() != 0,
2992  "output coord_2 is not defined");
2993 
2994  for(i= 0; i < (*m_exp).size(); ++i)
2995  {
2996  e_coord_0 = coord_0 + m_phys_offset[i];
2997  e_coord_1 = coord_1 + m_phys_offset[i];
2998  e_coord_2 = coord_2 + m_phys_offset[i];
2999  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1,e_coord_2);
3000  }
3001  break;
3002  }
3003  }
3004 
3005  /**
3006  */
3008  {
3010  "This method is not defined or valid for this class type");
3011  }
3012 
3013  /**
3014  */
3016  std::shared_ptr<ExpList> &result,
3017  const bool DeclareCoeffPhysArrays)
3018  {
3019  boost::ignore_unused(i, result, DeclareCoeffPhysArrays);
3021  "This method is not defined or valid for this class type");
3022  }
3023 
3024  /**
3025  */
3027  const Array<OneD, NekDouble> &element,
3028  Array<OneD, NekDouble> &boundary)
3029  {
3030  int n, cnt;
3031  Array<OneD, NekDouble> tmp1, tmp2;
3033 
3034  Array<OneD, int> ElmtID,EdgeID;
3035  GetBoundaryToElmtMap(ElmtID,EdgeID);
3036 
3037  // Initialise result
3038  boundary = Array<OneD, NekDouble>
3039  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3040 
3041  // Skip other boundary regions
3042  for (cnt = n = 0; n < i; ++n)
3043  {
3044  cnt += GetBndCondExpansions()[n]->GetExpSize();
3045  }
3046 
3047  int offsetBnd;
3048  int offsetElmt = 0;
3049  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3050  {
3051  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3052 
3053  elmt = GetExp(ElmtID[cnt+n]);
3054  elmt->GetTracePhysVals(EdgeID[cnt+n],
3055  GetBndCondExpansions()[i]->GetExp(n),
3056  tmp1 = element + offsetElmt,
3057  tmp2 = boundary + offsetBnd);
3058 
3059  offsetElmt += elmt->GetTotPoints();
3060  }
3061  }
3062 
3063  /**
3064  */
3066  const Array<OneD, const NekDouble> &phys,
3067  Array<OneD, NekDouble> &bndElmt)
3068  {
3069  int n, cnt, nq;
3070 
3071  Array<OneD, int> ElmtID,EdgeID;
3072  GetBoundaryToElmtMap(ElmtID,EdgeID);
3073 
3074  // Skip other boundary regions
3075  for (cnt = n = 0; n < i; ++n)
3076  {
3077  cnt += GetBndCondExpansions()[n]->GetExpSize();
3078  }
3079 
3080  // Count number of points
3081  int npoints = 0;
3082  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3083  {
3084  npoints += GetExp(ElmtID[cnt+n])->GetTotPoints();
3085  }
3086 
3087  // Initialise result
3088  bndElmt = Array<OneD, NekDouble> (npoints, 0.0);
3089 
3090  // Extract data
3091  int offsetPhys;
3092  int offsetElmt = 0;
3093  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3094  {
3095  nq = GetExp(ElmtID[cnt+n])->GetTotPoints();
3096  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
3097  Vmath::Vcopy(nq, &phys[offsetPhys], 1,
3098  &bndElmt[offsetElmt], 1);
3099  offsetElmt += nq;
3100  }
3101  }
3102 
3103  /**
3104  */
3106  const Array<OneD, const NekDouble> &phys,
3108  {
3109  int n, cnt;
3112 
3113  Array<OneD, int> ElmtID,EdgeID;
3114  GetBoundaryToElmtMap(ElmtID,EdgeID);
3115 
3116  // Initialise result
3118  (GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3119 
3120  // Skip other boundary regions
3121  for (cnt = n = 0; n < i; ++n)
3122  {
3123  cnt += GetBndCondExpansions()[n]->GetExpSize();
3124  }
3125 
3126  int offsetBnd;
3127  int offsetPhys;
3128  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3129  {
3130  offsetPhys = GetPhys_Offset(ElmtID[cnt+n]);
3131  offsetBnd = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3132 
3133  elmt = GetExp(ElmtID[cnt+n]);
3134  elmt->GetTracePhysVals(EdgeID[cnt+n],
3135  GetBndCondExpansions()[i]->GetExp(n),
3136  phys + offsetPhys,
3137  tmp1 = bnd + offsetBnd);
3138  }
3139  }
3140 
3141  /**
3142  */
3144  Array<OneD, Array<OneD, NekDouble> > &normals)
3145  {
3146  int j, n, cnt, nq;
3147  int coordim = GetCoordim(0);
3150 
3151  Array<OneD, int> ElmtID,EdgeID;
3152  GetBoundaryToElmtMap(ElmtID,EdgeID);
3153 
3154  // Initialise result
3155  normals = Array<OneD, Array<OneD, NekDouble> > (coordim);
3156  for (j = 0; j < coordim; ++j)
3157  {
3158  normals[j] = Array<OneD, NekDouble> (
3159  GetBndCondExpansions()[i]->GetTotPoints(), 0.0);
3160  }
3161 
3162  // Skip other boundary regions
3163  for (cnt = n = 0; n < i; ++n)
3164  {
3165  cnt += GetBndCondExpansions()[n]->GetExpSize();
3166  }
3167 
3168  int offset;
3169  for (n = 0; n < GetBndCondExpansions()[i]->GetExpSize(); ++n)
3170  {
3171  offset = GetBndCondExpansions()[i]->GetPhys_Offset(n);
3172  nq = GetBndCondExpansions()[i]->GetExp(n)->GetTotPoints();
3173 
3174  elmt = GetExp(ElmtID[cnt+n]);
3175  const Array<OneD, const Array<OneD, NekDouble> > normalsElmt
3176  = elmt->GetSurfaceNormal(EdgeID[cnt+n]);
3177  // Copy to result
3178  for (j = 0; j < coordim; ++j)
3179  {
3180  Vmath::Vcopy(nq, normalsElmt[j], 1,
3181  tmp = normals[j] + offset, 1);
3182  }
3183  }
3184  }
3185 
3186  /**
3187  */
3189  Array<OneD,int> &EdgeID)
3190  {
3191  boost::ignore_unused(ElmtID, EdgeID);
3193  "This method is not defined or valid for this class type");
3194  }
3195 
3196  /**
3197  */
3199  {
3201  "This method is not defined or valid for this class type");
3202  }
3203 
3204  /**
3205  */
3208  {
3210  "This method is not defined or valid for this class type");
3212  result;
3213  return result;
3214  }
3215 
3216  /**
3217  */
3219  {
3221  "This method is not defined or valid for this class type");
3223  return result;
3224  }
3225 
3226  /**
3227  */
3229  const NekDouble time,
3230  const std::string varName,
3231  const NekDouble x2_in,
3232  const NekDouble x3_in)
3233  {
3234  boost::ignore_unused(time, varName, x2_in, x3_in);
3236  "This method is not defined or valid for this class type");
3237  }
3238 
3239  /**
3240  */
3241  map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
3242  {
3244  "This method is not defined or valid for this class type");
3245  static map<int,RobinBCInfoSharedPtr> result;
3246  return result;
3247  }
3248 
3249  /**
3250  */
3252  PeriodicMap &periodicVerts,
3253  PeriodicMap &periodicEdges,
3254  PeriodicMap &periodicFaces)
3255  {
3256  boost::ignore_unused(periodicVerts, periodicEdges, periodicFaces);
3258  "This method is not defined or valid for this class type");
3259  }
3260 
3263  unsigned int regionId,
3264  const std::string& variable)
3265  {
3266  auto collectionIter = collection.find(regionId);
3267  ASSERTL1(collectionIter != collection.end(),
3268  "Unable to locate collection " +
3269  boost::lexical_cast<string>(regionId));
3270 
3272  = (*collectionIter).second;
3273  auto conditionMapIter = bndCondMap->find(variable);
3274  ASSERTL1(conditionMapIter != bndCondMap->end(),
3275  "Unable to locate condition map.");
3276 
3277  const SpatialDomains::BoundaryConditionShPtr boundaryCondition
3278  = (*conditionMapIter).second;
3279 
3280  return boundaryCondition;
3281  }
3282 
3284  {
3285  boost::ignore_unused(n);
3287  "This method is not defined or valid for this class type");
3288  return NullExpListSharedPtr;
3289  }
3290 
3291  /**
3292  * @brief Construct collections of elements containing a single element
3293  * type and polynomial order from the list of expansions.
3294  */
3296  {
3298  vector<std::pair<LocalRegions::ExpansionSharedPtr,int> > > collections;
3299 
3300  // Figure out optimisation parameters if provided in
3301  // session file or default given
3303  ImpType = colOpt.GetDefaultImplementationType();
3304 
3305  bool autotuning = colOpt.IsUsingAutotuning();
3306  bool verbose = (m_session->DefinesCmdLineArgument("verbose")) &&
3307  (m_comm->GetRank() == 0);
3308  int collmax = (colOpt.GetMaxCollectionSize() > 0
3309  ? colOpt.GetMaxCollectionSize()
3310  : 2*m_exp->size());
3311 
3312  // clear vectors in case previously called
3313  m_collections.clear();
3314  m_coll_coeff_offset.clear();
3315  m_coll_phys_offset.clear();
3316 
3317  // Loop over expansions, and create collections for each element type
3318  for (int i = 0; i < m_exp->size(); ++i)
3319  {
3320  collections[(*m_exp)[i]->DetShapeType()].push_back(
3321  std::pair<LocalRegions::ExpansionSharedPtr,int> ((*m_exp)[i],i));
3322  }
3323 
3324  for (auto &it : collections)
3325  {
3326  LocalRegions::ExpansionSharedPtr exp = it.second[0].first;
3327 
3328  Collections::OperatorImpMap impTypes = colOpt.GetOperatorImpMap(exp);
3329  vector<StdRegions::StdExpansionSharedPtr> collExp;
3330 
3331  int prevCoeffOffset = m_coeff_offset[it.second[0].second];
3332  int prevPhysOffset = m_phys_offset [it.second[0].second];
3333  int collcnt;
3334 
3335  m_coll_coeff_offset.push_back(prevCoeffOffset);
3336  m_coll_phys_offset .push_back(prevPhysOffset);
3337 
3338  if(it.second.size() == 1) // single element case
3339  {
3340  collExp.push_back(it.second[0].first);
3341 
3342  // if no Imp Type provided and No settign in xml file.
3343  // reset impTypes using timings
3344  if(autotuning)
3345  {
3346  impTypes = colOpt.SetWithTimings(collExp,
3347  impTypes, verbose);
3348  }
3349 
3350  Collections::Collection tmp(collExp, impTypes);
3351  m_collections.push_back(tmp);
3352  }
3353  else
3354  {
3355  // set up first geometry
3356  collExp.push_back(it.second[0].first);
3357  int prevnCoeff = it.second[0].first->GetNcoeffs();
3358  int prevnPhys = it.second[0].first->GetTotPoints();
3359  collcnt = 1;
3360 
3361  for (int i = 1; i < it.second.size(); ++i)
3362  {
3363  int nCoeffs = it.second[i].first->GetNcoeffs();
3364  int nPhys = it.second[i].first->GetTotPoints();
3365  int coeffOffset = m_coeff_offset[it.second[i].second];
3366  int physOffset = m_phys_offset [it.second[i].second];
3367 
3368  // check to see if next elmt is different or
3369  // collmax reached and if so end collection
3370  // and start new one
3371  if(prevCoeffOffset + nCoeffs != coeffOffset ||
3372  prevnCoeff != nCoeffs ||
3373  prevPhysOffset + nPhys != physOffset ||
3374  prevnPhys != nPhys || collcnt >= collmax)
3375  {
3376 
3377  // if no Imp Type provided and No
3378  // settign in xml file. reset
3379  // impTypes using timings
3380  if(autotuning)
3381  {
3382  impTypes = colOpt.SetWithTimings(collExp,
3383  impTypes,
3384  verbose);
3385  }
3386 
3387  Collections::Collection tmp(collExp, impTypes);
3388  m_collections.push_back(tmp);
3389 
3390 
3391  // start new geom list
3392  collExp.clear();
3393 
3394  m_coll_coeff_offset.push_back(coeffOffset);
3395  m_coll_phys_offset .push_back(physOffset);
3396  collExp.push_back(it.second[i].first);
3397  collcnt = 1;
3398  }
3399  else // add to list of collections
3400  {
3401  collExp.push_back(it.second[i].first);
3402  collcnt++;
3403  }
3404 
3405  // if end of list finish up collection
3406  if (i == it.second.size() - 1)
3407  {
3408  // if no Imp Type provided and No
3409  // settign in xml file.
3410  if(autotuning)
3411  {
3412  impTypes = colOpt.SetWithTimings(collExp,
3413  impTypes,verbose);
3414  }
3415 
3416  Collections::Collection tmp(collExp, impTypes);
3417  m_collections.push_back(tmp);
3418  collExp.clear();
3419  collcnt = 0;
3420 
3421  }
3422 
3423  prevCoeffOffset = coeffOffset;
3424  prevPhysOffset = physOffset;
3425  prevnCoeff = nCoeffs;
3426  prevnPhys = nPhys;
3427  }
3428  }
3429  }
3430  }
3431 
3433  {
3435  }
3436 
3437  } //end of namespace
3438 } //end of namespace
3439 
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Definition: ExpList.cpp:1958
virtual void v_GlobalToLocal(void)
Definition: ExpList.cpp:2885
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2629
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:2511
virtual void v_GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2933
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:2346
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:1442
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Definition: ExpList.cpp:1583
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:2333
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:3143
virtual std::shared_ptr< ExpList > & v_GetTrace()
Definition: ExpList.cpp:2581
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:2852
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:1649
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:2917
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:2174
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:2602
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:1570
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:1868
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
Definition: ExpList.cpp:2787
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
Definition: ExpList.cpp:2029
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:72
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:3065
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:1610
virtual NekDouble v_GetHomoLen(void)
Definition: ExpList.cpp:2037
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
Definition: ExpList.cpp:3105
void WriteVtkFooter(std::ostream &outfile)
Definition: ExpList.cpp:1862
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
Definition: ExpList.cpp:2954
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:2021
virtual const Array< OneD, const int > & v_GetTraceBndMap()
Definition: ExpList.cpp:2597
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
Definition: ExpList.cpp:1882
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
Definition: ExpList.cpp:3015
const char *const GlobalSysSolnTypeMap[]
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
Definition: ExpList.cpp:2778
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
Definition: Conditions.h:225
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble > > &inarray)
Definition: ExpList.cpp:2002
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
virtual void v_LocalToGlobal(bool UseComm)
Definition: ExpList.cpp:2867
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:2843
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:2559
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:1919
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:97
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:1988
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
Definition: ExpList.cpp:2282
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:1743
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:2735
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:3228
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
Definition: ExpList.cpp:1545
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:2435
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:2459
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:2589
double NekDouble
virtual void v_SetHomoLen(const NekDouble lhom)
Definition: ExpList.cpp:2045
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:2542
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
Definition: ExpList.cpp:3218
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
Definition: ExpList.cpp:2095
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
virtual void v_ClearGlobalLinSysManager(void)
Definition: ExpList.cpp:2089
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:2757
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
Definition: ExpList.cpp:3283
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:3251
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:2069
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:2610
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:2721
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
Definition: ExpList.cpp:2052
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:3241
#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:2767
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:2155
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:1875
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:179
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
Definition: ExpList.cpp:2657
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:2707
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:2296
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
Definition: ExpList.cpp:3261
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:3295
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:2682
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:3188
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
Definition: ExpList.cpp:2060
Used to lookup the create function in NekManager.
Definition: MatrixKey.h:67
virtual void v_SetUpPhysNormals()
Definition: ExpList.cpp:3007
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:3207
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
Definition: ExpList.cpp:2550
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:3026
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:1854
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:2666
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: ExpList.cpp:2079
#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:2901
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:2639
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:2324
std::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
Definition: ExpList.cpp:1820
Collections::CollectionVector m_collections
Definition: ExpList.h:1092
virtual void v_ReadGlobalOptimizationParameters()
Definition: ExpList.cpp:3198
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:2746
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:2692
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Definition: ExpList.cpp:2909