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