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