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