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