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>
37 #include <LibUtilities/Communication/Comm.h>
38 #include <MultiRegions/GlobalLinSys.h>
39 
40 #include <LocalRegions/MatrixKey.h> // for MatrixKey
41 #include <LocalRegions/Expansion.h> // for Expansion
42 
43 #include <MultiRegions/AssemblyMap/AssemblyMapCG.h> // for AssemblyMapCG, etc
44 #include <MultiRegions/AssemblyMap/AssemblyMapDG.h> // for AssemblyMapCG, etc
45 #include <MultiRegions/GlobalLinSysKey.h> // for GlobalLinSysKey
46 #include <MultiRegions/GlobalMatrix.h> // for GlobalMatrix, etc
47 #include <MultiRegions/GlobalMatrixKey.h> // for GlobalMatrixKey
48 
49 #include <LibUtilities/LinearAlgebra/SparseMatrixFwd.hpp>
50 #include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>
51 #include <LibUtilities/LinearAlgebra/NekMatrix.hpp>
52 
53 
54 namespace Nektar
55 {
56  namespace MultiRegions
57  {
58  /**
59  * @class ExpList
60  * All multi-elemental expansions \f$u^{\delta}(\boldsymbol{x})\f$ can
61  * be considered as the assembly of the various elemental contributions.
62  * On a discrete level, this yields,
63  * \f[u^{\delta}(\boldsymbol{x}_i)=\sum_{e=1}^{{N_{\mathrm{el}}}}
64  * \sum_{n=0}^{N^{e}_m-1}\hat{u}_n^e\phi_n^e(\boldsymbol{x}_i).\f]
65  * where \f${N_{\mathrm{el}}}\f$ is the number of elements and
66  * \f$N^{e}_m\f$ is the local elemental number of expansion modes.
67  * As it is the lowest level class, it contains the definition of the
68  * common data and common routines to all multi-elemental expansions.
69  *
70  * The class stores a vector of expansions, \a m_exp, (each derived from
71  * StdRegions#StdExpansion) which define the constituent components of
72  * the domain. The coefficients from these expansions are concatenated
73  * in \a m_coeffs, while the expansion evaluated at the quadrature
74  * points is stored in \a m_phys.
75  */
76 
77  /**
78  * Creates an empty expansion list. The expansion list will typically be
79  * populated by a derived class (namely one of MultiRegions#ExpList1D,
80  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
81  */
82  ExpList::ExpList():
83  m_comm(),
84  m_session(),
85  m_graph(),
86  m_ncoeffs(0),
87  m_npoints(0),
88  m_coeffs(),
89  m_phys(),
90  m_physState(false),
91  m_exp(MemoryManager<LocalRegions::ExpansionVector>
92  ::AllocateSharedPtr()),
93  m_coeff_offset(),
94  m_phys_offset(),
95  m_offset_elmt_id(),
96  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
97  m_WaveSpace(false)
98  {
99  SetExpType(eNoType);
100  }
101 
102 
103  /**
104  * Creates an empty expansion list. The expansion list will typically be
105  * populated by a derived class (namely one of MultiRegions#ExpList1D,
106  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
107  */
108  ExpList::ExpList(
109  const LibUtilities::SessionReaderSharedPtr &pSession):
110  m_comm(pSession->GetComm()),
111  m_session(pSession),
112  m_graph(),
113  m_ncoeffs(0),
114  m_npoints(0),
115  m_coeffs(),
116  m_phys(),
117  m_physState(false),
118  m_exp(MemoryManager<LocalRegions::ExpansionVector>
119  ::AllocateSharedPtr()),
120  m_coeff_offset(),
121  m_phys_offset(),
122  m_offset_elmt_id(),
123  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
124  m_WaveSpace(false)
125  {
126  SetExpType(eNoType);
127  }
128 
129 
130  /**
131  * Creates an empty expansion list. The expansion list will typically be
132  * populated by a derived class (namely one of MultiRegions#ExpList1D,
133  * MultiRegions#ExpList2D or MultiRegions#ExpList3D).
134  */
135  ExpList::ExpList(
136  const LibUtilities::SessionReaderSharedPtr &pSession,
137  const SpatialDomains::MeshGraphSharedPtr &pGraph):
138  m_comm(pSession->GetComm()),
139  m_session(pSession),
140  m_graph(pGraph),
141  m_ncoeffs(0),
142  m_npoints(0),
143  m_coeffs(),
144  m_phys(),
145  m_physState(false),
146  m_exp(MemoryManager<LocalRegions::ExpansionVector>
147  ::AllocateSharedPtr()),
148  m_coeff_offset(),
149  m_phys_offset(),
150  m_offset_elmt_id(),
151  m_blockMat(MemoryManager<BlockMatrixMap>::AllocateSharedPtr()),
152  m_WaveSpace(false)
153  {
154  SetExpType(eNoType);
155  }
156 
157 
158  /**
159  * Copies an existing expansion list.
160  * @param in Source expansion list.
161  */
162  ExpList::ExpList(const ExpList &in, const bool DeclareCoeffPhysArrays):
163  m_comm(in.m_comm),
164  m_session(in.m_session),
165  m_graph(in.m_graph),
166  m_ncoeffs(in.m_ncoeffs),
167  m_npoints(in.m_npoints),
168  m_physState(false),
169  m_exp(in.m_exp),
170  m_coeff_offset(in.m_coeff_offset),
171  m_phys_offset(in.m_phys_offset),
172  m_offset_elmt_id(in.m_offset_elmt_id),
173  m_globalOptParam(in.m_globalOptParam),
174  m_blockMat(in.m_blockMat),
175  m_WaveSpace(false)
176  {
177  SetExpType(eNoType);
178 
179  if(DeclareCoeffPhysArrays)
180  {
181  m_coeffs = Array<OneD, NekDouble>(m_ncoeffs);
182  m_phys = Array<OneD, NekDouble>(m_npoints);
183  }
184  }
185 
186  /**
187  *
188  */
189  ExpansionType ExpList::GetExpType(void)
190  {
191  return m_expType;
192  }
193 
194  /**
195  *
196  */
197  void ExpList::SetExpType(ExpansionType Type)
198  {
199  m_expType = Type;
200  }
201 
202  ExpList::~ExpList()
203  {
204  }
205 
206  /**
207  * The integration is evaluated locally, that is
208  * \f[\int
209  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
210  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
211  * where the integration over the separate elements is done by the
212  * function StdRegions#StdExpansion#Integral, which discretely
213  * evaluates the integral using Gaussian quadrature.
214  *
215  * Note that the array #m_phys should be filled with the values of the
216  * function \f$f(\boldsymbol{x})\f$ at the quadrature points
217  * \f$\boldsymbol{x}_i\f$.
218  *
219  * @return The value of the discretely evaluated integral
220  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
221  */
222  NekDouble ExpList::PhysIntegral()
223  {
224  ASSERTL1(m_physState == true,
225  "local physical space is not true ");
226 
227  return PhysIntegral(m_phys);
228  }
229 
230 
231  /**
232  * The integration is evaluated locally, that is
233  * \f[\int
234  * f(\boldsymbol{x})d\boldsymbol{x}=\sum_{e=1}^{{N_{\mathrm{el}}}}
235  * \left\{\int_{\Omega_e}f(\boldsymbol{x})d\boldsymbol{x}\right\}, \f]
236  * where the integration over the separate elements is done by the
237  * function StdRegions#StdExpansion#Integral, which discretely
238  * evaluates the integral using Gaussian quadrature.
239  *
240  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
241  * containing the values of the function
242  * \f$f(\boldsymbol{x})\f$ at the quadrature
243  * points \f$\boldsymbol{x}_i\f$.
244  * @return The value of the discretely evaluated integral
245  * \f$\int f(\boldsymbol{x})d\boldsymbol{x}\f$.
246  */
247  NekDouble ExpList::PhysIntegral(
248  const Array<OneD, const NekDouble> &inarray)
249  {
250  int i;
251  NekDouble sum = 0.0;
252 
253  for(i = 0; i < (*m_exp).size(); ++i)
254  {
255  sum += (*m_exp)[i]->Integral(inarray + m_phys_offset[i]);
256  }
257 
258  return sum;
259  }
260 
261 
262  /**
263  * Retrieves the block matrix specified by \a bkey, and computes
264  * \f$ y=Mx \f$.
265  * @param gkey GlobalMatrixKey specifying the block matrix to
266  * use in the matrix-vector multiply.
267  * @param inarray Input vector \f$ x \f$.
268  * @param outarray Output vector \f$ y \f$.
269  */
270  void ExpList::MultiplyByBlockMatrix(
271  const GlobalMatrixKey &gkey,
272  const Array<OneD,const NekDouble> &inarray,
273  Array<OneD, NekDouble> &outarray)
274  {
275  // Retrieve the block matrix using the given key.
276  const DNekScalBlkMatSharedPtr& blockmat = GetBlockMatrix(gkey);
277  int nrows = blockmat->GetRows();
278  int ncols = blockmat->GetColumns();
279 
280  // Create NekVectors from the given data arrays
281  NekVector<NekDouble> in (ncols,inarray, eWrapper);
282  NekVector< NekDouble> out(nrows,outarray,eWrapper);
283 
284  // Perform matrix-vector multiply.
285  out = (*blockmat)*in;
286  }
287 
288 
289  /**
290  * The operation is evaluated locally for every element by the function
291  * StdRegions#StdExpansion#IProductWRTBase.
292  *
293  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
294  * containing the values of the function
295  * \f$f(\boldsymbol{x})\f$ at the quadrature
296  * points \f$\boldsymbol{x}_i\f$.
297  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
298  * used to store the result.
299  */
300  void ExpList::v_IProductWRTBase_IterPerExp(
301  const Array<OneD, const NekDouble> &inarray,
302  Array<OneD, NekDouble> &outarray)
303  {
304  // get optimisation information about performing block
305  // matrix multiplies
306  const Array<OneD, const bool> doBlockMatOp
307  = m_globalOptParam->DoBlockMatOp(StdRegions::eIProductWRTBase);
308  const Array<OneD, LibUtilities::ShapeType> shape = m_globalOptParam->GetShapeList();
309  const Array<OneD, const int> num_elmts = m_globalOptParam->GetShapeNumElements();
310 
311  Array<OneD,NekDouble> tmp_outarray;
312  int cnt = 0,eid;
313 
314  for(int n = 0; n < shape.num_elements(); ++n)
315  {
316  if(doBlockMatOp[n])
317  {
318  if(num_elmts[n])
319  {
320  GlobalMatrixKey mkey(StdRegions::eIProductWRTBase,
321  shape[n]);
322  eid = m_offset_elmt_id[cnt];
323  MultiplyByBlockMatrix(mkey,inarray + m_phys_offset[eid],
324  tmp_outarray = outarray + m_coeff_offset[eid]);
325  cnt += num_elmts[n];
326  }
327  }
328  else
329  {
330  int i;
331  for(i = 0; i < num_elmts[n]; ++i)
332  {
333  eid = m_offset_elmt_id[cnt++];
334  (*m_exp)[eid]->IProductWRTBase(inarray+m_phys_offset[eid],
335  tmp_outarray = outarray+m_coeff_offset[eid]);
336  }
337  }
338  }
339  }
340 
341  /**
342  * The operation is evaluated locally for every element by the function
343  * StdRegions#StdExpansion#IProductWRTDerivBase.
344  *
345  * @param dir {0,1} is the direction in which the
346  * derivative of the basis should be taken
347  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
348  * containing the values of the function
349  * \f$f(\boldsymbol{x})\f$ at the quadrature
350  * points \f$\boldsymbol{x}_i\f$.
351  * @param outarray An array of size \f$N_{\mathrm{eof}}\f$
352  * used to store the result.
353  */
354  void ExpList::IProductWRTDerivBase(const int dir,
355  const Array<OneD, const NekDouble> &inarray,
356  Array<OneD, NekDouble> &outarray)
357  {
358  int i;
359 
360  Array<OneD,NekDouble> e_outarray;
361 
362  for(i = 0; i < (*m_exp).size(); ++i)
363  {
364  (*m_exp)[i]->IProductWRTDerivBase(dir,inarray+m_phys_offset[i],
365  e_outarray = outarray+m_coeff_offset[i]);
366  }
367  }
368 
369  /**
370  * Given a function \f$f(\boldsymbol{x})\f$ evaluated at
371  * the quadrature points, this function calculates the
372  * derivatives \f$\frac{d}{dx_1}\f$, \f$\frac{d}{dx_2}\f$
373  * and \f$\frac{d}{dx_3}\f$ of the function
374  * \f$f(\boldsymbol{x})\f$ at the same quadrature
375  * points. The local distribution of the quadrature points
376  * allows an elemental evaluation of the derivative. This
377  * is done by a call to the function
378  * StdRegions#StdExpansion#PhysDeriv.
379  *
380  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
381  * containing the values of the function
382  * \f$f(\boldsymbol{x})\f$ at the quadrature
383  * points \f$\boldsymbol{x}_i\f$.
384  * @param out_d0 The discrete evaluation of the
385  * derivative\f$\frac{d}{dx_1}\f$ will
386  * be stored in this array of size
387  * \f$Q_{\mathrm{tot}}\f$.
388  * @param out_d1 The discrete evaluation of the
389  * derivative\f$\frac{d}{dx_2}\f$ will be
390  * stored in this array of size
391  * \f$Q_{\mathrm{tot}}\f$. Note that if no
392  * memory is allocated for \a out_d1, the
393  * derivative \f$\frac{d}{dx_2}\f$ will not be
394  * calculated.
395  * @param out_d2 The discrete evaluation of the
396  * derivative\f$\frac{d}{dx_3}\f$ will be
397  * stored in this array of size
398  * \f$Q_{\mathrm{tot}}\f$. Note that if no
399  * memory is allocated for \a out_d2, the
400  * derivative \f$\frac{d}{dx_3}\f$ will not be
401  * calculated.
402  */
403  void ExpList::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
404  Array<OneD, NekDouble> &out_d0,
405  Array<OneD, NekDouble> &out_d1,
406  Array<OneD, NekDouble> &out_d2)
407  {
408  int i;
409  Array<OneD, NekDouble> e_out_d0;
410  Array<OneD, NekDouble> e_out_d1;
411  Array<OneD, NekDouble> e_out_d2;
412 
413  for(i= 0; i < (*m_exp).size(); ++i)
414  {
415  e_out_d0 = out_d0 + m_phys_offset[i];
416  if(out_d1.num_elements())
417  {
418  e_out_d1 = out_d1 + m_phys_offset[i];
419  }
420 
421  if(out_d2.num_elements())
422  {
423  e_out_d2 = out_d2 + m_phys_offset[i];
424  }
425  (*m_exp)[i]->PhysDeriv(inarray+m_phys_offset[i],e_out_d0,e_out_d1,e_out_d2);
426  }
427  }
428 
429  void ExpList::v_PhysDeriv(const int dir,
430  const Array<OneD, const NekDouble> &inarray,
431  Array<OneD, NekDouble> &out_d)
432  {
433  Direction edir = DirCartesianMap[dir];
434  v_PhysDeriv(edir, inarray,out_d);
435  }
436 
437  void ExpList::v_PhysDeriv(Direction edir, const Array<OneD, const NekDouble> &inarray,
438  Array<OneD, NekDouble> &out_d)
439  {
440  int i;
441  if(edir==MultiRegions::eS)
442  {
443  Array<OneD, NekDouble> e_out_ds;
444  for(i=0; i<(*m_exp).size(); ++i)
445  {
446  e_out_ds = out_d + m_phys_offset[i];
447  (*m_exp)[i]->PhysDeriv_s(inarray+m_phys_offset[i],e_out_ds);
448  }
449  }
450  else if(edir==MultiRegions::eN)
451  {
452  Array<OneD, NekDouble > e_out_dn;
453  for(i=0; i<(*m_exp).size(); i++)
454  {
455  e_out_dn = out_d +m_phys_offset[i];
456  (*m_exp)[i]->PhysDeriv_n(inarray+m_phys_offset[i],e_out_dn);
457  }
458  }
459  else
460  {
461  // convert enum into int
462  int intdir= (int)edir;
463  Array<OneD, NekDouble> e_out_d;
464  for(i= 0; i < (*m_exp).size(); ++i)
465  {
466  e_out_d = out_d + m_phys_offset[i];
467  (*m_exp)[i]->PhysDeriv(intdir, inarray+m_phys_offset[i], e_out_d);
468  }
469 
470  }
471  }
472 
473 
474  /**
475  * The coefficients of the function to be acted upon
476  * should be contained in the \param inarray. The
477  * resulting coefficients are stored in \param outarray
478  *
479  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
480  * containing the inner product.
481  */
482  void ExpList::MultiplyByElmtInvMass(
483  const Array<OneD, const NekDouble> &inarray,
484  Array<OneD, NekDouble> &outarray)
485  {
486  GlobalMatrixKey mkey(StdRegions::eInvMass);
487  const DNekScalBlkMatSharedPtr& InvMass = GetBlockMatrix(mkey);
488 
489  // Inverse mass matrix
490  NekVector<NekDouble> out(m_ncoeffs,outarray,eWrapper);
491  if(inarray.get() == outarray.get())
492  {
493  NekVector<NekDouble> in(m_ncoeffs,inarray); // copy data
494  out = (*InvMass)*in;
495  }
496  else
497  {
498  NekVector<NekDouble> in(m_ncoeffs,inarray,eWrapper);
499  out = (*InvMass)*in;
500  }
501  }
502 
503  /**
504  * Given a function \f$u(\boldsymbol{x})\f$ defined at the
505  * quadrature points, this function determines the
506  * transformed elemental coefficients \f$\hat{u}_n^e\f$
507  * employing a discrete elemental Galerkin projection from
508  * physical space to coefficient space. For each element,
509  * the operation is evaluated locally by the function
510  * StdRegions#StdExpansion#IproductWRTBase followed by a
511  * call to #MultiRegions#MultiplyByElmtInvMass.
512  *
513  * @param inarray An array of size \f$Q_{\mathrm{tot}}\f$
514  * containing the values of the function
515  * \f$f(\boldsymbol{x})\f$ at the quadrature
516  * points \f$\boldsymbol{x}_i\f$.
517  * @param outarray The resulting coefficients
518  * \f$\hat{u}_n^e\f$ will be stored in this
519  * array of size \f$N_{\mathrm{eof}}\f$.
520  */
521  void ExpList::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
522  Array<OneD, NekDouble> &outarray)
523  {
524  Array<OneD,NekDouble> f(m_ncoeffs);
525 
526  IProductWRTBase_IterPerExp(inarray,f);
527  MultiplyByElmtInvMass(f,outarray);
528 
529  }
530 
531  void ExpList::FwdTrans_BndConstrained(
532  const Array<OneD, const NekDouble>& inarray,
533  Array<OneD, NekDouble> &outarray)
534  {
535  int i;
536 
537  Array<OneD,NekDouble> e_outarray;
538 
539  for(i= 0; i < (*m_exp).size(); ++i)
540  {
541  (*m_exp)[i]->FwdTrans_BndConstrained(inarray+m_phys_offset[i],
542  e_outarray = outarray+m_coeff_offset[i]);
543  }
544  }
545 
546  /**
547  * This function smooth a field after some calculaitons which have
548  * been done elementally.
549  *
550  * @param field An array containing the field in physical space
551  *
552  */
553  void ExpList::v_SmoothField(Array<OneD, NekDouble> &field)
554  {
555  // Do nothing unless the method is implemented in the appropriate
556  // class, i.e. ContField1D,ContField2D, etc.
557 
558  // So far it has been implemented just for ContField2D and
559  // ContField3DHomogeneous1D
560 
561  // Block in case users try the smoothing with a modal expansion.
562  // Maybe a different techique for the smoothing require
563  // implementation for modal basis.
564 
565  ASSERTL0((*m_exp)[0]->GetBasisType(0)
566  == LibUtilities::eGLL_Lagrange ||
567  (*m_exp)[0]->GetBasisType(0)
568  == LibUtilities::eGauss_Lagrange,
569  "Smoothing is currently not allowed unless you are using "
570  "a nodal base for efficiency reasons. The implemented "
571  "smoothing technique requires the mass matrix inversion "
572  "which is trivial just for GLL_LAGRANGE_SEM and "
573  "GAUSS_LAGRANGE_SEMexpansions.");
574  }
575 
576 
577  /**
578  * This function assembles the block diagonal matrix
579  * \f$\underline{\boldsymbol{M}}^e\f$, which is the
580  * concatenation of the local matrices
581  * \f$\boldsymbol{M}^e\f$ of the type \a mtype, that is
582  *
583  * \f[
584  * \underline{\boldsymbol{M}}^e = \left[
585  * \begin{array}{cccc}
586  * \boldsymbol{M}^1 & 0 & \hspace{3mm}0 \hspace{3mm}& 0 \\
587  * 0 & \boldsymbol{M}^2 & 0 & 0 \\
588  * 0 & 0 & \ddots & 0 \\
589  * 0 & 0 & 0 & \boldsymbol{M}^{N_{\mathrm{el}}} \end{array}\right].\f]
590  *
591  * @param mtype the type of matrix to be assembled
592  * @param scalar an optional parameter
593  * @param constant an optional parameter
594  */
595  const DNekScalBlkMatSharedPtr ExpList::GenBlockMatrix(
596  const GlobalMatrixKey &gkey)
597  {
598  int i,cnt1;
599  int n_exp = 0;
600  DNekScalMatSharedPtr loc_mat;
601  DNekScalBlkMatSharedPtr BlkMatrix;
602  map<int,int> elmt_id;
603  LibUtilities::ShapeType ShapeType = gkey.GetShapeType();
604 
605  if(ShapeType != LibUtilities::eNoShapeType)
606  {
607  for(i = 0 ; i < (*m_exp).size(); ++i)
608  {
609  if((*m_exp)[m_offset_elmt_id[i]]->DetShapeType()
610  == ShapeType)
611  {
612  elmt_id[n_exp++] = m_offset_elmt_id[i];
613  }
614  }
615  }
616  else
617  {
618  n_exp = (*m_exp).size();
619  for(i = 0; i < n_exp; ++i)
620  {
621  elmt_id[i] = m_offset_elmt_id[i];
622  }
623  }
624 
625  Array<OneD,unsigned int> nrows(n_exp);
626  Array<OneD,unsigned int> ncols(n_exp);
627 
628  switch(gkey.GetMatrixType())
629  {
630  case StdRegions::eBwdTrans:
631  {
632  // set up an array of integers for block matrix construction
633  for(i = 0; i < n_exp; ++i)
634  {
635  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
636  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
637  }
638  }
639  break;
640  case StdRegions::eIProductWRTBase:
641  {
642  // set up an array of integers for block matrix construction
643  for(i = 0; i < n_exp; ++i)
644  {
645  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
646  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
647  }
648  }
649  break;
650  case StdRegions::eMass:
651  case StdRegions::eInvMass:
652  case StdRegions::eHelmholtz:
653  case StdRegions::eLaplacian:
654  case StdRegions::eInvHybridDGHelmholtz:
655  {
656  // set up an array of integers for block matrix construction
657  for(i = 0; i < n_exp; ++i)
658  {
659  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
660  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
661  }
662  }
663  break;
664 
665  case StdRegions::eHybridDGLamToU:
666  {
667  // set up an array of integers for block matrix construction
668  for(i = 0; i < n_exp; ++i)
669  {
670  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
671  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
672  }
673  }
674  break;
675 
676  case StdRegions::eHybridDGHelmBndLam:
677  {
678  // set up an array of integers for block matrix construction
679  for(i = 0; i < n_exp; ++i)
680  {
681  nrows[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
682  ncols[i] = (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
683  }
684  }
685  break;
686 
687  default:
688  {
689  NEKERROR(ErrorUtil::efatal,
690  "Global Matrix creation not defined for this type "
691  "of matrix");
692  }
693  }
694 
695  MatrixStorage blkmatStorage = eDIAGONAL;
696  BlkMatrix = MemoryManager<DNekScalBlkMat>
697  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
698 
699  int nvarcoeffs = gkey.GetNVarCoeffs();
700  int eid;
701  Array<OneD, NekDouble> varcoeffs_wk;
702 
703  for(i = cnt1 = 0; i < n_exp; ++i)
704  {
705  // need to be initialised with zero size for non variable coefficient case
706  StdRegions::VarCoeffMap varcoeffs;
707 
708  eid = elmt_id[i];
709  if(nvarcoeffs>0)
710  {
711  StdRegions::VarCoeffMap::const_iterator x;
712  for (x = gkey.GetVarCoeffs().begin(); x != gkey.GetVarCoeffs().end(); ++x)
713  {
714  varcoeffs[x->first] = x->second + m_phys_offset[eid];
715  }
716  }
717 
718  LocalRegions::MatrixKey matkey(gkey.GetMatrixType(),
719  (*m_exp)[eid]->DetShapeType(),
720  *(*m_exp)[eid],
721  gkey.GetConstFactors(),
722  varcoeffs );
723 
724  loc_mat = boost::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[elmt_id.find(i)->second])->GetLocMatrix(matkey);
725  BlkMatrix->SetBlock(i,i,loc_mat);
726  }
727 
728  return BlkMatrix;
729  }
730 
731  const DNekScalBlkMatSharedPtr& ExpList::GetBlockMatrix(
732  const GlobalMatrixKey &gkey)
733  {
734  BlockMatrixMap::iterator matrixIter = m_blockMat->find(gkey);
735 
736  if(matrixIter == m_blockMat->end())
737  {
738  return ((*m_blockMat)[gkey] = GenBlockMatrix(gkey));
739  }
740  else
741  {
742  return matrixIter->second;
743  }
744  }
745 
746  void ExpList::GeneralMatrixOp_IterPerExp(
747  const GlobalMatrixKey &gkey,
748  const Array<OneD,const NekDouble> &inarray,
749  Array<OneD, NekDouble> &outarray)
750  {
751  const Array<OneD, const bool> doBlockMatOp
752  = m_globalOptParam->DoBlockMatOp(gkey.GetMatrixType());
753  const Array<OneD, const int> num_elmts
754  = m_globalOptParam->GetShapeNumElements();
755 
756  Array<OneD,NekDouble> tmp_outarray;
757  int cnt = 0;
758  int eid;
759  for(int n = 0; n < num_elmts.num_elements(); ++n)
760  {
761  if(doBlockMatOp[n])
762  {
763  const LibUtilities::ShapeType vType
764  = m_globalOptParam->GetShapeList()[n];
765  const MultiRegions::GlobalMatrixKey vKey(gkey, vType);
766  if (cnt < m_offset_elmt_id.num_elements())
767  {
768  eid = m_offset_elmt_id[cnt];
769  MultiplyByBlockMatrix(vKey,inarray + m_coeff_offset[eid],
770  tmp_outarray = outarray + m_coeff_offset[eid]);
771  cnt += num_elmts[n];
772  }
773  }
774  else
775  {
776  int i;
777  int nvarcoeffs = gkey.GetNVarCoeffs();
778 
779  for(i= 0; i < num_elmts[n]; ++i)
780  {
781  // need to be initialised with zero size for non variable coefficient case
782  StdRegions::VarCoeffMap varcoeffs;
783 
784  eid = m_offset_elmt_id[cnt++];
785  if(nvarcoeffs>0)
786  {
787  StdRegions::VarCoeffMap::const_iterator x;
788  for (x = gkey.GetVarCoeffs().begin(); x != gkey.GetVarCoeffs().end(); ++x)
789  {
790  varcoeffs[x->first] = x->second + m_phys_offset[eid];
791  }
792  }
793 
794  StdRegions::StdMatrixKey mkey(gkey.GetMatrixType(),
795  (*m_exp)[eid]->DetShapeType(),
796  *((*m_exp)[eid]),
797  gkey.GetConstFactors(),varcoeffs);
798 
799  (*m_exp)[eid]->GeneralMatrixOp(inarray + m_coeff_offset[eid],
800  tmp_outarray = outarray+m_coeff_offset[eid],
801  mkey);
802  }
803  }
804  }
805  }
806 
807  /**
808  * Retrieves local matrices from each expansion in the expansion list
809  * and combines them together to generate a global matrix system.
810  * @param mkey Matrix key for the matrix to be generated.
811  * @param locToGloMap Local to global mapping.
812  * @returns Shared pointer to the generated global matrix.
813  */
814  GlobalMatrixSharedPtr ExpList::GenGlobalMatrix(
815  const GlobalMatrixKey &mkey,
816  const AssemblyMapCGSharedPtr &locToGloMap)
817  {
818  int i,j,n,gid1,gid2,cntdim1,cntdim2;
819  NekDouble sign1,sign2;
820  DNekScalMatSharedPtr loc_mat;
821 
822  unsigned int glob_rows;
823  unsigned int glob_cols;
824  unsigned int loc_rows;
825  unsigned int loc_cols;
826 
827  bool assembleFirstDim;
828  bool assembleSecondDim;
829 
830  switch(mkey.GetMatrixType())
831  {
832  case StdRegions::eBwdTrans:
833  {
834  glob_rows = m_npoints;
835  glob_cols = locToGloMap->GetNumGlobalCoeffs();
836 
837  assembleFirstDim = false;
838  assembleSecondDim = true;
839  }
840  break;
841  case StdRegions::eIProductWRTBase:
842  {
843  glob_rows = locToGloMap->GetNumGlobalCoeffs();
844  glob_cols = m_npoints;
845 
846  assembleFirstDim = true;
847  assembleSecondDim = false;
848  }
849  break;
850  case StdRegions::eMass:
851  case StdRegions::eHelmholtz:
852  case StdRegions::eLaplacian:
853  case StdRegions::eHybridDGHelmBndLam:
854  {
855  glob_rows = locToGloMap->GetNumGlobalCoeffs();
856  glob_cols = locToGloMap->GetNumGlobalCoeffs();
857 
858  assembleFirstDim = true;
859  assembleSecondDim = true;
860  }
861  break;
862  default:
863  {
864  NEKERROR(ErrorUtil::efatal,
865  "Global Matrix creation not defined for this type "
866  "of matrix");
867  }
868  }
869 
870  COOMatType spcoomat;
871  CoordType coord;
872 
873  int nvarcoeffs = mkey.GetNVarCoeffs();
874  int eid;
875 
876  // fill global matrix
877  for(n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
878  {
879  // need to be initialised with zero size for non variable coefficient case
880  StdRegions::VarCoeffMap varcoeffs;
881 
882  eid = m_offset_elmt_id[n];
883  if(nvarcoeffs>0)
884  {
885  StdRegions::VarCoeffMap::const_iterator x;
886  for (x = mkey.GetVarCoeffs().begin(); x != mkey.GetVarCoeffs().end(); ++x)
887  {
888  varcoeffs[x->first] = x->second + m_phys_offset[eid];
889  }
890  }
891 
892  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
893  (*m_exp)[eid]->DetShapeType(),
894  *((*m_exp)[eid]),
895  mkey.GetConstFactors(),varcoeffs);
896 
897  loc_mat = boost::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[m_offset_elmt_id[n]])->GetLocMatrix(matkey);
898 
899  loc_rows = loc_mat->GetRows();
900  loc_cols = loc_mat->GetColumns();
901 
902  for(i = 0; i < loc_rows; ++i)
903  {
904  if(assembleFirstDim)
905  {
906  gid1 = locToGloMap->GetLocalToGlobalMap (cntdim1 + i);
907  sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
908  }
909  else
910  {
911  gid1 = cntdim1 + i;
912  sign1 = 1.0;
913  }
914 
915  for(j = 0; j < loc_cols; ++j)
916  {
917  if(assembleSecondDim)
918  {
919  gid2 = locToGloMap
920  ->GetLocalToGlobalMap(cntdim2 + j);
921  sign2 = locToGloMap
922  ->GetLocalToGlobalSign(cntdim2 + j);
923  }
924  else
925  {
926  gid2 = cntdim2 + j;
927  sign2 = 1.0;
928  }
929 
930  // sparse matrix fill
931  coord = make_pair(gid1,gid2);
932  if( spcoomat.count(coord) == 0 )
933  {
934  spcoomat[coord] = sign1*sign2*(*loc_mat)(i,j);
935  }
936  else
937  {
938  spcoomat[coord] += sign1*sign2*(*loc_mat)(i,j);
939  }
940  }
941  }
942  cntdim1 += loc_rows;
943  cntdim2 += loc_cols;
944  }
945 
946  return MemoryManager<GlobalMatrix>
947  ::AllocateSharedPtr(m_session,glob_rows,glob_cols,spcoomat);
948  }
949 
950 
951  DNekMatSharedPtr ExpList::GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const AssemblyMapCGSharedPtr &locToGloMap)
952  {
953  int i,j,n,gid1,gid2,loc_lda,eid;
954  NekDouble sign1,sign2,value;
955  DNekScalMatSharedPtr loc_mat;
956 
957  int totDofs = locToGloMap->GetNumGlobalCoeffs();
958  int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
959 
960  unsigned int rows = totDofs - NumDirBCs;
961  unsigned int cols = totDofs - NumDirBCs;
962  NekDouble zero = 0.0;
963 
964  DNekMatSharedPtr Gmat;
965  int bwidth = locToGloMap->GetFullSystemBandWidth();
966 
967  int nvarcoeffs = mkey.GetNVarCoeffs();
968  MatrixStorage matStorage;
969 
970  map<int, RobinBCInfoSharedPtr> RobinBCInfo = GetRobinBCInfo();
971 
972  switch(mkey.GetMatrixType())
973  {
974  // case for all symmetric matices
975  case StdRegions::eHelmholtz:
976  case StdRegions::eLaplacian:
977  if( (2*(bwidth+1)) < rows)
978  {
979  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC_BANDED;
980  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage,bwidth,bwidth);
981  }
982  else
983  {
984  matStorage = ePOSITIVE_DEFINITE_SYMMETRIC;
985  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
986  }
987 
988  break;
989  default: // Assume general matrix - currently only set up for full invert
990  {
991  matStorage = eFULL;
992  Gmat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,zero,matStorage);
993  }
994  }
995 
996  // fill global symmetric matrix
997  for(n = 0; n < (*m_exp).size(); ++n)
998  {
999  // need to be initialised with zero size for non variable coefficient case
1000  StdRegions::VarCoeffMap varcoeffs;
1001 
1002  eid = m_offset_elmt_id[n];
1003  if(nvarcoeffs>0)
1004  {
1005  StdRegions::VarCoeffMap::const_iterator x;
1006  for (x = mkey.GetVarCoeffs().begin(); x != mkey.GetVarCoeffs().end(); ++x)
1007  {
1008  varcoeffs[x->first] = x->second + m_phys_offset[eid];
1009  }
1010  }
1011 
1012  LocalRegions::MatrixKey matkey(mkey.GetMatrixType(),
1013  (*m_exp)[eid]->DetShapeType(),
1014  *((*m_exp)[eid]),
1015  mkey.GetConstFactors(),varcoeffs);
1016 
1017  loc_mat = boost::dynamic_pointer_cast<LocalRegions::Expansion>((*m_exp)[n])->GetLocMatrix(matkey);
1018 
1019 
1020  if(RobinBCInfo.count(n) != 0) // add robin mass matrix
1021  {
1022  RobinBCInfoSharedPtr rBC;
1023 
1024  // declare local matrix from scaled matrix.
1025  int rows = loc_mat->GetRows();
1026  int cols = loc_mat->GetColumns();
1027  const NekDouble *dat = loc_mat->GetRawPtr();
1028  DNekMatSharedPtr new_mat = MemoryManager<DNekMat>::AllocateSharedPtr(rows,cols,dat);
1029  Blas::Dscal(rows*cols,loc_mat->Scale(),new_mat->GetRawPtr(),1);
1030 
1031  // add local matrix contribution
1032  for(rBC = RobinBCInfo.find(n)->second;rBC; rBC = rBC->next)
1033  {
1034  (*m_exp)[n]->AddRobinMassMatrix(rBC->m_robinID,rBC->m_robinPrimitiveCoeffs,new_mat);
1035  }
1036 
1037  NekDouble one = 1.0;
1038  // redeclare loc_mat to point to new_mat plus the scalar.
1039  loc_mat = MemoryManager<DNekScalMat>::AllocateSharedPtr(one,new_mat);
1040  }
1041 
1042  loc_lda = loc_mat->GetColumns();
1043 
1044  for(i = 0; i < loc_lda; ++i)
1045  {
1046  gid1 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + i) - NumDirBCs;
1047  sign1 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + i);
1048  if(gid1 >= 0)
1049  {
1050  for(j = 0; j < loc_lda; ++j)
1051  {
1052  gid2 = locToGloMap->GetLocalToGlobalMap(m_coeff_offset[n] + j) - NumDirBCs;
1053  sign2 = locToGloMap->GetLocalToGlobalSign(m_coeff_offset[n] + j);
1054  if(gid2 >= 0)
1055  {
1056  // When global matrix is symmetric,
1057  // only add the value for the upper
1058  // triangular part in order to avoid
1059  // entries to be entered twice
1060  if((matStorage == eFULL)||(gid2 >= gid1))
1061  {
1062  value = Gmat->GetValue(gid1,gid2) + sign1*sign2*(*loc_mat)(i,j);
1063  Gmat->SetValue(gid1,gid2,value);
1064  }
1065  }
1066  }
1067  }
1068  }
1069  }
1070 
1071  return Gmat;
1072  }
1073 
1074 
1075  /**
1076  * Consider a linear system
1077  * \f$\boldsymbol{M\hat{u}}_g=\boldsymbol{f}\f$ to be solved. Dependent
1078  * on the solution method, this function constructs
1079  * - <b>The full linear system</b><BR>
1080  * A call to the function #GenGlobalLinSysFullDirect
1081  * - <b>The statically condensed linear system</b><BR>
1082  * A call to the function #GenGlobalLinSysStaticCond
1083  *
1084  * @param mkey A key which uniquely defines the global
1085  * matrix to be constructed.
1086  * @param locToGloMap Contains the mapping array and required
1087  * information for the transformation from
1088  * local to global degrees of freedom.
1089  * @return (A shared pointer to) the global linear system in
1090  * required format.
1091  */
1092  GlobalLinSysSharedPtr ExpList::GenGlobalLinSys(
1093  const GlobalLinSysKey &mkey,
1094  const AssemblyMapCGSharedPtr &locToGloMap)
1095  {
1096  GlobalLinSysSharedPtr returnlinsys;
1097  boost::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
1098 
1099  MultiRegions::GlobalSysSolnType vType = mkey.GetGlobalSysSolnType();
1100 
1101  if (vType >= eSIZE_GlobalSysSolnType)
1102  {
1103  ASSERTL0(false,"Matrix solution type not defined");
1104  }
1105  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
1106 
1107  return GetGlobalLinSysFactory().CreateInstance( vSolnType, mkey,
1108  vExpList, locToGloMap);
1109  }
1110 
1111  GlobalLinSysSharedPtr ExpList::GenGlobalBndLinSys(
1112  const GlobalLinSysKey &mkey,
1113  const AssemblyMapSharedPtr &locToGloMap)
1114  {
1115  boost::shared_ptr<ExpList> vExpList = GetSharedThisPtr();
1116  const map<int,RobinBCInfoSharedPtr> vRobinBCInfo = GetRobinBCInfo();
1117 
1118  MultiRegions::GlobalSysSolnType vType = mkey.GetGlobalSysSolnType();
1119 
1120  if (vType >= eSIZE_GlobalSysSolnType)
1121  {
1122  ASSERTL0(false,"Matrix solution type not defined");
1123  }
1124  std::string vSolnType = MultiRegions::GlobalSysSolnTypeMap[vType];
1125 
1126  return GetGlobalLinSysFactory().CreateInstance(vSolnType,mkey,
1127  vExpList,locToGloMap);
1128  }
1129 
1130 
1131  /**
1132  * Given the elemental coefficients \f$\hat{u}_n^e\f$ of
1133  * an expansion, this function evaluates the spectral/hp
1134  * expansion \f$u^{\delta}(\boldsymbol{x})\f$ at the
1135  * quadrature points \f$\boldsymbol{x}_i\f$. The operation
1136  * is evaluated locally by the elemental function
1137  * StdRegions#StdExpansion#BwdTrans.
1138  *
1139  * @param inarray An array of size \f$N_{\mathrm{eof}}\f$
1140  * containing the local coefficients
1141  * \f$\hat{u}_n^e\f$.
1142  * @param outarray The resulting physical values at the
1143  * quadrature points
1144  * \f$u^{\delta}(\boldsymbol{x}_i)\f$
1145  * will be stored in this array of size
1146  * \f$Q_{\mathrm{tot}}\f$.
1147  */
1148  void ExpList::v_BwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray,
1149  Array<OneD, NekDouble> &outarray)
1150  {
1151  // get optimisation information about performing block
1152  // matrix multiplies
1153  const Array<OneD, const bool> doBlockMatOp
1154  = m_globalOptParam->DoBlockMatOp(StdRegions::eBwdTrans);
1155  const Array<OneD, LibUtilities::ShapeType> shape = m_globalOptParam->GetShapeList();
1156  const Array<OneD, const int> num_elmts = m_globalOptParam->GetShapeNumElements();
1157 
1158  Array<OneD,NekDouble> tmp_outarray;
1159  int cnt = 0,eid;
1160 
1161  for(int n = 0; n < num_elmts.num_elements(); ++n)
1162  {
1163  if(doBlockMatOp[n])
1164  {
1165  if(num_elmts[n])
1166  {
1167  GlobalMatrixKey mkey(StdRegions::eBwdTrans, shape[n]);
1168  eid = m_offset_elmt_id[cnt];
1169  MultiplyByBlockMatrix(mkey,inarray + m_coeff_offset[eid],
1170  tmp_outarray = outarray + m_phys_offset[eid]);
1171  cnt += num_elmts[n];
1172  }
1173  }
1174  else
1175  {
1176  int i;
1177 
1178  for(i= 0; i < num_elmts[n]; ++i)
1179  {
1180  eid = m_offset_elmt_id[cnt++];
1181  (*m_exp)[eid]->BwdTrans(inarray + m_coeff_offset[eid],
1182  tmp_outarray = outarray+m_phys_offset[eid]);
1183  }
1184  }
1185  }
1186  }
1187 
1188  LocalRegions::ExpansionSharedPtr& ExpList::GetExp(
1189  const Array<OneD, const NekDouble> &gloCoord)
1190  {
1191  Array<OneD, NekDouble> stdCoord(GetCoordim(0),0.0);
1192  for (int i = 0; i < (*m_exp).size(); ++i)
1193  {
1194  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoord))
1195  {
1196  return (*m_exp)[i];
1197  }
1198  }
1199  ASSERTL0(false, "Cannot find element for this point.");
1200  return (*m_exp)[0]; // avoid warnings
1201  }
1202 
1203 
1204  /**
1205  * @todo need a smarter search here that first just looks at bounding
1206  * vertices - suggest first seeing if point is within 10% of
1207  * region defined by vertices. The do point search.
1208  */
1209  int ExpList::GetExpIndex(
1210  const Array<OneD, const NekDouble> &gloCoord,
1211  NekDouble tol,
1212  bool returnNearestElmt)
1213  {
1214  Array<OneD, NekDouble> Lcoords(gloCoord.num_elements());
1215 
1216  return GetExpIndex(gloCoord,Lcoords,tol,returnNearestElmt);
1217  }
1218 
1219 
1220  int ExpList::GetExpIndex(const Array<OneD, const NekDouble> &gloCoords,
1221  Array<OneD, NekDouble> &locCoords,
1222  NekDouble tol,
1223  bool returnNearestElmt)
1224  {
1225  NekDouble resid;
1226 
1227  if (GetNumElmts() == 0)
1228  {
1229  return -1;
1230  }
1231  std::vector<std::pair<int,NekDouble> > elmtIdDist;
1232 
1233  // Manifold case (point may match multiple elements)
1234  if (GetExp(0)->GetCoordim() > GetExp(0)->GetShapeDimension())
1235  {
1236  SpatialDomains::PointGeomSharedPtr v;
1237  SpatialDomains::PointGeom w;
1238  NekDouble x, y, z;
1239 
1240  // Scan all elements and store those which may contain the point
1241  for (int i = 0; i < (*m_exp).size(); ++i)
1242  {
1243  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
1244  tol, resid))
1245  {
1246  v = m_graph->GetVertex((*m_exp)[i]->GetGeom()->GetVid(0));
1247 
1248  w.SetX(gloCoords[0]);
1249  w.SetY(gloCoords[1]);
1250  w.SetZ(gloCoords[2]);
1251  v->GetCoords(x,y,z);
1252 
1253  elmtIdDist.push_back(std::pair<int, NekDouble>(i, v->dist(w)));
1254  }
1255  }
1256 
1257  // Find nearest element
1258  if (!elmtIdDist.empty())
1259  {
1260  NekDouble min_d = elmtIdDist[0].first;
1261  int min_id = elmtIdDist[0].second;
1262 
1263  for (int i = 1; i < elmtIdDist.size(); ++i)
1264  {
1265  if (elmtIdDist[i].second < min_d) {
1266  min_d = elmtIdDist[i].second;
1267  min_id = elmtIdDist[i].first;
1268  }
1269  }
1270 
1271  return min_id;
1272  }
1273  else
1274  {
1275  return -1;
1276  }
1277  }
1278  // non-embedded mesh (point can only match one element)
1279  else
1280  {
1281  static int start = 0;
1282  int min_id = 0;
1283  NekDouble resid_min = 1e6;
1284  Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
1285 
1286  // restart search from last found value
1287  for (int i = start; i < (*m_exp).size(); ++i)
1288  {
1289  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
1290  tol, resid))
1291  {
1292  start = i;
1293  return i;
1294  }
1295  else
1296  {
1297  if(resid < resid_min)
1298  {
1299  min_id = i;
1300  resid_min = resid;
1301  Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
1302  }
1303  }
1304  }
1305 
1306  for (int i = 0; i < start; ++i)
1307  {
1308  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
1309  tol, resid))
1310  {
1311  start = i;
1312  return i;
1313  }
1314  else
1315  {
1316  if(resid < resid_min)
1317  {
1318  min_id = i;
1319  resid_min = resid;
1320  Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
1321  }
1322  }
1323  }
1324 
1325  std::string msg = "Failed to find point in element to tolerance of "
1326  + boost::lexical_cast<std::string>(resid)
1327  + " using nearest point found";
1328  WARNINGL0(true,msg.c_str());
1329 
1330  if(returnNearestElmt)
1331  {
1332  Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
1333  return min_id;
1334  }
1335  else
1336  {
1337  return -1;
1338  }
1339 
1340  }
1341  }
1342 
1343 
1344  /**
1345  * Configures geometric info, such as tangent direction, on each
1346  * expansion.
1347  * @param graph2D Mesh
1348  */
1349  void ExpList::ApplyGeomInfo()
1350  {
1351 
1352  }
1353 
1354  /**
1355  * Write Tecplot Files Header
1356  * @param outfile Output file name.
1357  * @param var variables names
1358  */
1359  void ExpList::v_WriteTecplotHeader(std::ofstream &outfile,
1360  std::string var)
1361  {
1362  int coordim = GetExp(0)->GetCoordim();
1363  char vars[3] = { 'x', 'y', 'z' };
1364 
1365  if (m_expType == e3DH1D)
1366  {
1367  coordim += 1;
1368  }
1369  else if (m_expType == e3DH2D)
1370  {
1371  coordim += 2;
1372  }
1373 
1374  outfile << "Variables = x";
1375  for (int i = 1; i < coordim; ++i)
1376  {
1377  outfile << ", " << vars[i];
1378  }
1379 
1380  if (var.size() > 0)
1381  {
1382  outfile << ", " << var;
1383  }
1384 
1385  outfile << std::endl << std::endl;
1386  }
1387 
1388  /**
1389  * Write Tecplot Files Zone
1390  * @param outfile Output file name.
1391  * @param expansion Expansion that is considered
1392  */
1393  void ExpList::v_WriteTecplotZone(std::ofstream &outfile, int expansion)
1394  {
1395  int i, j;
1396  int coordim = GetCoordim(0);
1397  int nPoints = GetTotPoints();
1398  int nBases = (*m_exp)[0]->GetNumBases();
1399  int numBlocks = 0;
1400 
1401  Array<OneD, Array<OneD, NekDouble> > coords(3);
1402 
1403  if (expansion == -1)
1404  {
1405  nPoints = GetTotPoints();
1406 
1407  coords[0] = Array<OneD, NekDouble>(nPoints);
1408  coords[1] = Array<OneD, NekDouble>(nPoints);
1409  coords[2] = Array<OneD, NekDouble>(nPoints);
1410 
1411  GetCoords(coords[0], coords[1], coords[2]);
1412 
1413  for (i = 0; i < m_exp->size(); ++i)
1414  {
1415  int numInt = 1;
1416 
1417  for (j = 0; j < nBases; ++j)
1418  {
1419  numInt *= (*m_exp)[i]->GetNumPoints(j)-1;
1420  }
1421 
1422  numBlocks += numInt;
1423  }
1424  }
1425  else
1426  {
1427  nPoints = (*m_exp)[expansion]->GetTotPoints();
1428 
1429  coords[0] = Array<OneD, NekDouble>(nPoints);
1430  coords[1] = Array<OneD, NekDouble>(nPoints);
1431  coords[2] = Array<OneD, NekDouble>(nPoints);
1432 
1433  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
1434 
1435  numBlocks = 1;
1436  for (j = 0; j < nBases; ++j)
1437  {
1438  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j)-1;
1439  }
1440  }
1441 
1442  if (m_expType == e3DH1D)
1443  {
1444  nBases += 1;
1445  coordim += 1;
1446  int nPlanes = GetZIDs().num_elements();
1447  NekDouble tmp = numBlocks * (nPlanes-1.0) / nPlanes;
1448  numBlocks = (int)tmp;
1449  }
1450  else if (m_expType == e3DH2D)
1451  {
1452  nBases += 2;
1453  coordim += 1;
1454  }
1455 
1456  outfile << "Zone, N=" << nPoints << ", E="
1457  << numBlocks << ", F=FEBlock" ;
1458 
1459  switch(nBases)
1460  {
1461  case 2:
1462  outfile << ", ET=QUADRILATERAL" << std::endl;
1463  break;
1464  case 3:
1465  outfile << ", ET=BRICK" << std::endl;
1466  break;
1467  default:
1468  ASSERTL0(false,"Not set up for this type of output");
1469  break;
1470  }
1471 
1472  // Write out coordinates
1473  for (j = 0; j < coordim; ++j)
1474  {
1475  for (i = 0; i < nPoints; ++i)
1476  {
1477  outfile << coords[j][i] << " ";
1478  if (i % 1000 == 0 && i)
1479  {
1480  outfile << std::endl;
1481  }
1482  }
1483  outfile << std::endl;
1484  }
1485  }
1486 
1487  void ExpList::v_WriteTecplotConnectivity(std::ofstream &outfile,
1488  int expansion)
1489  {
1490  int i,j,k,l;
1491  int nbase = (*m_exp)[0]->GetNumBases();
1492  int cnt = 0;
1493 
1494  boost::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
1495 
1496  if (expansion != -1)
1497  {
1498  exp = boost::shared_ptr<LocalRegions::ExpansionVector>(
1499  new LocalRegions::ExpansionVector(1));
1500  (*exp)[0] = (*m_exp)[expansion];
1501  }
1502 
1503  if (nbase == 2)
1504  {
1505  for(i = 0; i < (*exp).size(); ++i)
1506  {
1507  const int np0 = (*exp)[i]->GetNumPoints(0);
1508  const int np1 = (*exp)[i]->GetNumPoints(1);
1509 
1510  for(j = 1; j < np1; ++j)
1511  {
1512  for(k = 1; k < np0; ++k)
1513  {
1514  outfile << cnt + (j-1)*np0 + k << " ";
1515  outfile << cnt + (j-1)*np0 + k+1 << " ";
1516  outfile << cnt + j *np0 + k+1 << " ";
1517  outfile << cnt + j *np0 + k << endl;
1518  }
1519  }
1520 
1521  cnt += np0*np1;
1522  }
1523  }
1524  else if (nbase == 3)
1525  {
1526  for(i = 0; i < (*exp).size(); ++i)
1527  {
1528  const int np0 = (*exp)[i]->GetNumPoints(0);
1529  const int np1 = (*exp)[i]->GetNumPoints(1);
1530  const int np2 = (*exp)[i]->GetNumPoints(2);
1531  const int np01 = np0*np1;
1532 
1533  for(j = 1; j < np2; ++j)
1534  {
1535  for(k = 1; k < np1; ++k)
1536  {
1537  for(l = 1; l < np0; ++l)
1538  {
1539  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l << " ";
1540  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l+1 << " ";
1541  outfile << cnt + (j-1)*np01 + k *np0 + l+1 << " ";
1542  outfile << cnt + (j-1)*np01 + k *np0 + l << " ";
1543  outfile << cnt + j *np01 + (k-1)*np0 + l << " ";
1544  outfile << cnt + j *np01 + (k-1)*np0 + l+1 << " ";
1545  outfile << cnt + j *np01 + k *np0 + l+1 << " ";
1546  outfile << cnt + j *np01 + k *np0 + l << endl;
1547  }
1548  }
1549  }
1550  cnt += np0*np1*np2;
1551  }
1552  }
1553  else
1554  {
1555  ASSERTL0(false,"Not set up for this dimension");
1556  }
1557  }
1558 
1559  /**
1560  * Write Tecplot Files Field
1561  * @param outfile Output file name.
1562  * @param expansion Expansion that is considered
1563  */
1564  void ExpList::v_WriteTecplotField(std::ofstream &outfile, int expansion)
1565  {
1566  if (expansion == -1)
1567  {
1568  int totpoints = GetTotPoints();
1569  if(m_physState == false)
1570  {
1571  BwdTrans(m_coeffs,m_phys);
1572  }
1573 
1574  for(int i = 0; i < totpoints; ++i)
1575  {
1576  outfile << m_phys[i] << " ";
1577  if(i % 1000 == 0 && i)
1578  {
1579  outfile << std::endl;
1580  }
1581  }
1582  outfile << std::endl;
1583 
1584  }
1585  else
1586  {
1587  int nPoints = (*m_exp)[expansion]->GetTotPoints();
1588 
1589  for (int i = 0; i < nPoints; ++i)
1590  {
1591  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
1592  }
1593 
1594  outfile << std::endl;
1595  }
1596  }
1597 
1598  void ExpList::WriteVtkHeader(std::ofstream &outfile)
1599  {
1600  outfile << "<?xml version=\"1.0\"?>" << endl;
1601  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1602  << "byte_order=\"LittleEndian\">" << endl;
1603  outfile << " <UnstructuredGrid>" << endl;
1604  }
1605 
1606  void ExpList::WriteVtkFooter(std::ofstream &outfile)
1607  {
1608  outfile << " </UnstructuredGrid>" << endl;
1609  outfile << "</VTKFile>" << endl;
1610  }
1611 
1612  void ExpList::v_WriteVtkPieceHeader(std::ofstream &outfile, int expansion)
1613  {
1614  ASSERTL0(false, "Routine not implemented for this expansion.");
1615  }
1616 
1617  void ExpList::WriteVtkPieceFooter(std::ofstream &outfile, int expansion)
1618  {
1619  outfile << " </PointData>" << endl;
1620  outfile << " </Piece>" << endl;
1621  }
1622 
1623  void ExpList::v_WriteVtkPieceData(std::ofstream &outfile, int expansion,
1624  std::string var)
1625  {
1626  int i;
1627  int nq = (*m_exp)[expansion]->GetTotPoints();
1628 
1629  // printing the fields of that zone
1630  outfile << " <DataArray type=\"Float64\" Name=\""
1631  << var << "\">" << endl;
1632  outfile << " ";
1633  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
1634  for(i = 0; i < nq; ++i)
1635  {
1636  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
1637  }
1638  outfile << endl;
1639  outfile << " </DataArray>" << endl;
1640  }
1641 
1642  /**
1643  * Given a spectral/hp approximation
1644  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1645  * (which should be contained in #m_phys), this function calculates the
1646  * \f$L_\infty\f$ error of this approximation with respect to an exact
1647  * solution. The local distribution of the quadrature points allows an
1648  * elemental evaluation of this operation through the functions
1649  * StdRegions#StdExpansion#Linf.
1650  *
1651  * The exact solution, also evaluated at the quadrature
1652  * points, should be contained in the variable #m_phys of
1653  * the ExpList object \a Sol.
1654  *
1655  * @param soln A 1D array, containing the discrete
1656  * evaluation of the exact solution at the
1657  * quadrature points in its array #m_phys.
1658  * @return The \f$L_\infty\f$ error of the approximation.
1659  */
1660  NekDouble ExpList::Linf(
1661  const Array<OneD, const NekDouble> &inarray,
1662  const Array<OneD, const NekDouble> &soln)
1663  {
1664  NekDouble err = 0.0;
1665 
1666  if (soln == NullNekDouble1DArray)
1667  {
1668  err = Vmath::Vmax(m_npoints, inarray, 1);
1669  }
1670  else
1671  {
1672  for (int i = 0; i < m_npoints; ++i)
1673  {
1674  err = max(err, abs(inarray[i] - soln[i]));
1675  }
1676  }
1677 
1678  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
1679 
1680  return err;
1681  }
1682 
1683  /**
1684  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
1685  * evaluated at the quadrature points (which should be contained in
1686  * #m_phys), this function calculates the \f$L_2\f$ error of this
1687  * approximation with respect to an exact solution. The local
1688  * distribution of the quadrature points allows an elemental evaluation
1689  * of this operation through the functions StdRegions#StdExpansion#L2.
1690  *
1691  * The exact solution, also evaluated at the quadrature points, should
1692  * be contained in the variable #m_phys of the ExpList object \a Sol.
1693  *
1694  * @param Sol An ExpList, containing the discrete
1695  * evaluation of the exact solution at the
1696  * quadrature points in its array #m_phys.
1697  * @return The \f$L_2\f$ error of the approximation.
1698  */
1699  NekDouble ExpList::v_L2(
1700  const Array<OneD, const NekDouble> &inarray,
1701  const Array<OneD, const NekDouble> &soln)
1702  {
1703  NekDouble err = 0.0, errl2;
1704  int i;
1705 
1706  if (soln == NullNekDouble1DArray)
1707  {
1708  for (i = 0; i < (*m_exp).size(); ++i)
1709  {
1710  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
1711  err += errl2*errl2;
1712  }
1713  }
1714  else
1715  {
1716  for (i = 0; i < (*m_exp).size(); ++i)
1717  {
1718  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
1719  soln + m_phys_offset[i]);
1720  err += errl2*errl2;
1721  }
1722  }
1723 
1724  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1725 
1726  return sqrt(err);
1727  }
1728 
1729  NekDouble ExpList::v_Integral(const Array<OneD, const NekDouble> &inarray)
1730  {
1731  NekDouble err = 0.0;
1732  int i = 0;
1733 
1734  for (i = 0; i < (*m_exp).size(); ++i)
1735  {
1736  err += (*m_exp)[m_offset_elmt_id[i]]->Integral(inarray + m_phys_offset[i]);
1737  }
1738  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1739 
1740  return err;
1741  }
1742 
1743  Array<OneD, const NekDouble> ExpList::v_HomogeneousEnergy (void)
1744  {
1745  ASSERTL0(false,
1746  "This method is not defined or valid for this class type");
1747  Array<OneD, NekDouble> NoEnergy(1,0.0);
1748  return NoEnergy;
1749  }
1750 
1751  LibUtilities::TranspositionSharedPtr ExpList::v_GetTransposition(void)
1752  {
1753  ASSERTL0(false,
1754  "This method is not defined or valid for this class type");
1755  LibUtilities::TranspositionSharedPtr trans;
1756 
1757  return trans;
1758  }
1759 
1760  NekDouble ExpList::v_GetHomoLen(void)
1761  {
1762  ASSERTL0(false,
1763  "This method is not defined or valid for this class type");
1764  NekDouble len = 0.0;
1765  return len;
1766  }
1767 
1768  Array<OneD, const unsigned int> ExpList::v_GetZIDs(void)
1769  {
1770  ASSERTL0(false,
1771  "This method is not defined or valid for this class type");
1772  Array<OneD, unsigned int> NoModes(1);
1773 
1774  return NoModes;
1775  }
1776 
1777  Array<OneD, const unsigned int> ExpList::v_GetYIDs(void)
1778  {
1779  ASSERTL0(false,
1780  "This method is not defined or valid for this class type");
1781  Array<OneD, unsigned int> NoModes(1);
1782 
1783  return NoModes;
1784  }
1785 
1786 
1787  void ExpList::v_PhysInterp1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
1788  {
1789  ASSERTL0(false,
1790  "This method is not defined or valid for this class type");
1791  }
1792 
1793  void ExpList::v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray) {
1794  ASSERTL0(false,
1795  "This method is not defined or valid for this class type");
1796  }
1797 
1798  void ExpList::ExtractFileBCs(
1799  const std::string &fileName,
1800  const std::string &varName,
1801  const boost::shared_ptr<ExpList> locExpList)
1802  {
1803  string varString = fileName.substr(0, fileName.find_last_of("."));
1804  int j, k, len = varString.length();
1805  varString = varString.substr(len-1, len);
1806 
1807  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1808  std::vector<std::vector<NekDouble> > FieldData;
1809 
1810  LibUtilities::FieldIO f(m_session->GetComm());
1811  f.Import(fileName, FieldDef, FieldData);
1812 
1813  bool found = false;
1814  for (j = 0; j < FieldDef.size(); ++j)
1815  {
1816  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
1817  {
1818  if (FieldDef[j]->m_fields[k] == varName)
1819  {
1820  // Copy FieldData into locExpList
1821  locExpList->ExtractDataToCoeffs(
1822  FieldDef[j], FieldData[j],
1823  FieldDef[j]->m_fields[k],
1824  locExpList->UpdateCoeffs());
1825  found = true;
1826  }
1827  }
1828  }
1829 
1830  ASSERTL0(found, "Could not find variable '"+varName+
1831  "' in file boundary condition "+fileName);
1832  locExpList->BwdTrans_IterPerExp(
1833  locExpList->GetCoeffs(),
1834  locExpList->UpdatePhys());
1835  }
1836 
1837  /**
1838  * Given a spectral/hp approximation
1839  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1840  * (which should be contained in #m_phys), this function calculates the
1841  * \f$H^1_2\f$ error of this approximation with respect to an exact
1842  * solution. The local distribution of the quadrature points allows an
1843  * elemental evaluation of this operation through the functions
1844  * StdRegions#StdExpansion#H1.
1845  *
1846  * The exact solution, also evaluated at the quadrature points, should
1847  * be contained in the variable #m_phys of the ExpList object \a Sol.
1848  *
1849  * @param soln An 1D array, containing the discrete evaluation
1850  * of the exact solution at the quadrature points.
1851  *
1852  * @return The \f$H^1_2\f$ error of the approximation.
1853  */
1854  NekDouble ExpList::H1(
1855  const Array<OneD, const NekDouble> &inarray,
1856  const Array<OneD, const NekDouble> &soln)
1857  {
1858  NekDouble err = 0.0, errh1;
1859  int i;
1860 
1861  for (i = 0; i < (*m_exp).size(); ++i)
1862  {
1863  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
1864  soln + m_phys_offset[i]);
1865  err += errh1*errh1;
1866  }
1867 
1868  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1869 
1870  return sqrt(err);
1871  }
1872 
1873  void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
1874  int NumHomoDir,
1875  Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
1876  std::vector<NekDouble> &HomoLen,
1877  std::vector<unsigned int> &HomoZIDs,
1878  std::vector<unsigned int> &HomoYIDs)
1879  {
1880  int startenum, endenum, s;
1881 
1882  ASSERTL1(NumHomoDir == HomoBasis.num_elements(),"Homogeneous basis is not the same length as NumHomoDir");
1883  ASSERTL1(NumHomoDir == HomoLen.size(),"Homogeneous length vector is not the same length as NumHomDir");
1884 
1885  // count number of shapes
1886  switch((*m_exp)[0]->GetShapeDimension())
1887  {
1888  case 1:
1889  startenum = (int) LibUtilities::eSegment;
1890  endenum = (int) LibUtilities::eSegment;
1891  break;
1892  case 2:
1893  startenum = (int) LibUtilities::eTriangle;
1894  endenum = (int) LibUtilities::eQuadrilateral;
1895  break;
1896  case 3:
1897  startenum = (int) LibUtilities::eTetrahedron;
1898  endenum = (int) LibUtilities::eHexahedron;
1899  break;
1900  }
1901 
1902  for(s = startenum; s <= endenum; ++s)
1903  {
1904  LibUtilities::ShapeType shape;
1905  std::vector<unsigned int> elementIDs;
1906  std::vector<LibUtilities::BasisType> basis;
1907  std::vector<unsigned int> numModes;
1908  std::vector<std::string> fields;
1909 
1910  bool first = true;
1911  bool UniOrder = true;
1912  int n;
1913 
1914  for(int i = 0; i < (*m_exp).size(); ++i)
1915  {
1916  if((*m_exp)[i]->GetGeom()->GetShapeType() == (LibUtilities::ShapeType) s)
1917  {
1918  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
1919  if(first)
1920  {
1921  shape = (LibUtilities::ShapeType) s;
1922  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
1923  {
1924  basis.push_back((*m_exp)[i]->GetBasis(j)->GetBasisType());
1925  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
1926  }
1927 
1928  // add homogeneous direction details if defined
1929  for(n = 0 ; n < NumHomoDir; ++n)
1930  {
1931  basis.push_back(HomoBasis[n]->GetBasisType());
1932  numModes.push_back(HomoBasis[n]->GetNumModes());
1933  }
1934 
1935  first = false;
1936  }
1937  else
1938  {
1939  ASSERTL0((*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],"Routine is not set up for multiple bases definitions");
1940 
1941  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
1942  {
1943  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
1944  if(numModes[j] != (*m_exp)[i]->GetBasis(j)->GetNumModes())
1945  {
1946  UniOrder = false;
1947  }
1948  }
1949  // add homogeneous direction details if defined
1950  for(n = 0 ; n < NumHomoDir; ++n)
1951  {
1952  numModes.push_back(HomoBasis[n]->GetNumModes());
1953  }
1954  }
1955  }
1956  }
1957 
1958 
1959  if(elementIDs.size() > 0)
1960  {
1961  LibUtilities::FieldDefinitionsSharedPtr fdef = MemoryManager<LibUtilities::FieldDefinitions>::AllocateSharedPtr(shape, elementIDs, basis, UniOrder, numModes,fields, NumHomoDir, HomoLen, HomoZIDs, HomoYIDs);
1962  fielddef.push_back(fdef);
1963  }
1964  }
1965  }
1966 
1967 
1968  //
1969  // Virtual functions
1970  //
1971  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::v_GetFieldDefinitions()
1972  {
1973  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
1974  v_GetFieldDefinitions(returnval);
1975  return returnval;
1976  }
1977 
1978  void ExpList::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
1979  {
1980  GeneralGetFieldDefinitions(fielddef);
1981  }
1982 
1983  //Append the element data listed in elements
1984  //fielddef->m_ElementIDs onto fielddata
1985  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
1986  {
1987  v_AppendFieldData(fielddef,fielddata,m_coeffs);
1988  }
1989 
1990  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
1991  {
1992  int i;
1993  // Determine mapping from element ids to location in
1994  // expansion list
1995  // Determine mapping from element ids to location in
1996  // expansion list
1997  map<int, int> ElmtID_to_ExpID;
1998  for(i = 0; i < (*m_exp).size(); ++i)
1999  {
2000  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2001  }
2002 
2003  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
2004  {
2005  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
2006  int datalen = (*m_exp)[eid]->GetNcoeffs();
2007  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]],&coeffs[m_coeff_offset[eid]]+datalen);
2008  }
2009 
2010  }
2011 
2012  /// Extract the data in fielddata into the coeffs
2013  void ExpList::ExtractDataToCoeffs(
2014  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
2015  std::vector<NekDouble> &fielddata,
2016  std::string &field,
2017  Array<OneD, NekDouble> &coeffs)
2018  {
2019  v_ExtractDataToCoeffs(fielddef,fielddata,field,coeffs);
2020  }
2021 
2022  void ExpList::ExtractCoeffsToCoeffs(const boost::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2023  {
2024  v_ExtractCoeffsToCoeffs(fromExpList,fromCoeffs,toCoeffs);
2025  }
2026 
2027  /**
2028  * @brief Extract data from raw field data into expansion list.
2029  *
2030  * @param fielddef Field definitions.
2031  * @param fielddata Data for associated field.
2032  * @param field Field variable name.
2033  * @param coeffs Resulting coefficient array.
2034  */
2035  void ExpList::v_ExtractDataToCoeffs(
2036  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
2037  std::vector<NekDouble> &fielddata,
2038  std::string &field,
2039  Array<OneD, NekDouble> &coeffs)
2040  {
2041  int i, expId;
2042  int offset = 0;
2043  int modes_offset = 0;
2044  int datalen = fielddata.size()/fielddef->m_fields.size();
2045 
2046  // Find data location according to field definition
2047  for(i = 0; i < fielddef->m_fields.size(); ++i)
2048  {
2049  if(fielddef->m_fields[i] == field)
2050  {
2051  break;
2052  }
2053  offset += datalen;
2054  }
2055 
2056  ASSERTL0(i != fielddef->m_fields.size(),
2057  "Field (" + field + ") not found in file.");
2058 
2059  // Determine mapping from element ids to location in expansion list
2060  map<int, int> elmtToExpId;
2061 
2062  // Loop in reverse order so that in case where using a Homogeneous
2063  // expansion it sets geometry ids to first part of m_exp
2064  // list. Otherwise will set to second (complex) expansion
2065  for(i = (*m_exp).size()-1; i >= 0; --i)
2066  {
2067  elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2068  }
2069 
2070  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
2071  {
2072  // Reset modes_offset in the case where all expansions of
2073  // the same order.
2074  if (fielddef->m_uniOrder == true)
2075  {
2076  modes_offset = 0;
2077  }
2078 
2079  datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType,
2080  fielddef->m_numModes, modes_offset);
2081 
2082  const int elmtId = fielddef->m_elementIDs[i];
2083  if (elmtToExpId.count(elmtId) == 0)
2084  {
2085  offset += datalen;
2086  continue;
2087  }
2088 
2089  expId = elmtToExpId[elmtId];
2090 
2091  if (datalen == (*m_exp)[expId]->GetNcoeffs())
2092  {
2093  Vmath::Vcopy(datalen, &fielddata[offset], 1,
2094  &coeffs[m_coeff_offset[expId]], 1);
2095  }
2096  else
2097  {
2098  (*m_exp)[expId]->ExtractDataToCoeffs(
2099  &fielddata[offset], fielddef->m_numModes,
2100  modes_offset, &coeffs[m_coeff_offset[expId]]);
2101  }
2102 
2103  offset += datalen;
2104  modes_offset += (*m_exp)[0]->GetNumBases();
2105  }
2106 
2107  return;
2108  }
2109 
2110  void ExpList::v_ExtractCoeffsToCoeffs(const boost::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2111  {
2112  int i;
2113  int offset = 0;
2114 
2115  // check if the same and if so just copy over coeffs
2116  if(fromExpList->GetNcoeffs() == m_ncoeffs)
2117  {
2118  Vmath::Vcopy(m_ncoeffs,fromCoeffs,1,toCoeffs,1);
2119  }
2120  else
2121  {
2122  std::vector<unsigned int> nummodes;
2123  for(i = 0; i < (*m_exp).size(); ++i)
2124  {
2125  int eid = m_offset_elmt_id[i];
2126  for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
2127  {
2128  nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j));
2129  }
2130 
2131  (*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
2132  &toCoeffs[m_coeff_offset[eid]]);
2133 
2134  offset += fromExpList->GetExp(eid)->GetNcoeffs();
2135  }
2136  }
2137  }
2138 
2139 
2140  const Array<OneD,const boost::shared_ptr<ExpList> >
2141  &ExpList::v_GetBndCondExpansions(void)
2142  {
2143  ASSERTL0(false,
2144  "This method is not defined or valid for this class type");
2145  static Array<OneD,const boost::shared_ptr<ExpList> > result;
2146  return result;
2147  }
2148 
2149  boost::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
2150  {
2151  ASSERTL0(false,
2152  "This method is not defined or valid for this class type");
2153  static boost::shared_ptr<ExpList> result;
2154  return result;
2155  }
2156 
2157  void ExpList::v_Upwind(
2158  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
2159  const Array<OneD, const NekDouble> &Fwd,
2160  const Array<OneD, const NekDouble> &Bwd,
2161  Array<OneD, NekDouble> &Upwind)
2162  {
2163  ASSERTL0(false,
2164  "This method is not defined or valid for this class type");
2165  }
2166 
2167  void ExpList::v_Upwind(
2168  const Array<OneD, const NekDouble> &Vn,
2169  const Array<OneD, const NekDouble> &Fwd,
2170  const Array<OneD, const NekDouble> &Bwd,
2171  Array<OneD, NekDouble> &Upwind)
2172  {
2173  ASSERTL0(false,
2174  "This method is not defined or valid for this class type");
2175  }
2176 
2177  boost::shared_ptr<ExpList> &ExpList::v_GetTrace()
2178  {
2179  ASSERTL0(false,
2180  "This method is not defined or valid for this class type");
2181  static boost::shared_ptr<ExpList> returnVal;
2182  return returnVal;
2183  }
2184 
2185  boost::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
2186  {
2187  ASSERTL0(false,
2188  "This method is not defined or valid for this class type");
2189  static boost::shared_ptr<AssemblyMapDG> result;
2190  return result;
2191  }
2192 
2193  const Array<OneD, const int> &ExpList::v_GetTraceBndMap()
2194  {
2195  return GetTraceMap()->GetBndCondTraceToGlobalTraceMap();
2196  }
2197 
2198  void ExpList::v_GetNormals(
2199  Array<OneD, Array<OneD, NekDouble> > &normals)
2200  {
2201  ASSERTL0(false,
2202  "This method is not defined or valid for this class type");
2203  }
2204 
2205  void ExpList::v_AddTraceIntegral(
2206  const Array<OneD, const NekDouble> &Fx,
2207  const Array<OneD, const NekDouble> &Fy,
2208  Array<OneD, NekDouble> &outarray)
2209  {
2210  ASSERTL0(false,
2211  "This method is not defined or valid for this class type");
2212  }
2213 
2214  void ExpList::v_AddTraceIntegral(
2215  const Array<OneD, const NekDouble> &Fn,
2216  Array<OneD, NekDouble> &outarray)
2217  {
2218  ASSERTL0(false,
2219  "This method is not defined or valid for this class type");
2220  }
2221 
2222  void ExpList::v_AddFwdBwdTraceIntegral(
2223  const Array<OneD, const NekDouble> &Fwd,
2224  const Array<OneD, const NekDouble> &Bwd,
2225  Array<OneD, NekDouble> &outarray)
2226  {
2227  ASSERTL0(false,
2228  "This method is not defined or valid for this class type");
2229  }
2230 
2231  void ExpList::v_GetFwdBwdTracePhys(Array<OneD,NekDouble> &Fwd,
2232  Array<OneD,NekDouble> &Bwd)
2233  {
2234  ASSERTL0(false,
2235  "This method is not defined or valid for this class type");
2236  }
2237 
2238  void ExpList::v_GetFwdBwdTracePhys(
2239  const Array<OneD,const NekDouble> &field,
2240  Array<OneD,NekDouble> &Fwd,
2241  Array<OneD,NekDouble> &Bwd)
2242  {
2243  ASSERTL0(false,
2244  "This method is not defined or valid for this class type");
2245  }
2246 
2247  void ExpList::v_ExtractTracePhys(Array<OneD,NekDouble> &outarray)
2248  {
2249  ASSERTL0(false,
2250  "This method is not defined or valid for this class type");
2251  }
2252 
2253  void ExpList::v_ExtractTracePhys(
2254  const Array<OneD, const NekDouble> &inarray,
2255  Array<OneD,NekDouble> &outarray)
2256  {
2257  ASSERTL0(false,
2258  "This method is not defined or valid for this class type");
2259  }
2260 
2261  void ExpList::v_MultiplyByInvMassMatrix(
2262  const Array<OneD,const NekDouble> &inarray,
2263  Array<OneD, NekDouble> &outarray,
2264  CoeffState coeffstate)
2265  {
2266  ASSERTL0(false,
2267  "This method is not defined or valid for this class type");
2268  }
2269 
2270  void ExpList::v_HelmSolve(
2271  const Array<OneD, const NekDouble> &inarray,
2272  Array<OneD, NekDouble> &outarray,
2273  const FlagList &flags,
2274  const StdRegions::ConstFactorMap &factors,
2275  const StdRegions::VarCoeffMap &varcoeff,
2276  const Array<OneD, const NekDouble> &dirForcing)
2277  {
2278  ASSERTL0(false, "HelmSolve not implemented.");
2279  }
2280 
2281  void ExpList::v_LinearAdvectionDiffusionReactionSolve(
2282  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2283  const Array<OneD, const NekDouble> &inarray,
2284  Array<OneD, NekDouble> &outarray,
2285  const NekDouble lambda,
2286  CoeffState coeffstate,
2287  const Array<OneD, const NekDouble>& dirForcing)
2288  {
2289  ASSERTL0(false,
2290  "This method is not defined or valid for this class type");
2291  }
2292 
2293  void ExpList::v_LinearAdvectionReactionSolve(
2294  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2295  const Array<OneD, const NekDouble> &inarray,
2296  Array<OneD, NekDouble> &outarray,
2297  const NekDouble lambda,
2298  CoeffState coeffstate,
2299  const Array<OneD, const NekDouble>& dirForcing)
2300  {
2301  ASSERTL0(false,
2302  "This method is not defined or valid for this class type");
2303  }
2304 
2305  void ExpList::v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
2306  Array<OneD, NekDouble> &outarray,
2307  CoeffState coeffstate,
2308  bool Shuff,
2309  bool UnShuff)
2310  {
2311  ASSERTL0(false,
2312  "This method is not defined or valid for this class type");
2313  }
2314 
2315  void ExpList::v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
2316  Array<OneD, NekDouble> &outarray,
2317  CoeffState coeffstate,
2318  bool Shuff,
2319  bool UnShuff)
2320  {
2321  ASSERTL0(false,
2322  "This method is not defined or valid for this class type");
2323  }
2324 
2325  void ExpList::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,const Array<OneD, NekDouble> &inarray2,Array<OneD, NekDouble> &outarray,CoeffState coeffstate)
2326  {
2327  ASSERTL0(false,
2328  "This method is not defined or valid for this class type");
2329  }
2330 
2331 
2332  void ExpList::v_GetBCValues(Array<OneD, NekDouble> &BndVals,
2333  const Array<OneD, NekDouble> &TotField,
2334  int BndID)
2335  {
2336  ASSERTL0(false,
2337  "This method is not defined or valid for this class type");
2338  }
2339 
2340  void ExpList::v_NormVectorIProductWRTBase(Array<OneD, const NekDouble> &V1,
2341  Array<OneD, const NekDouble> &V2,
2342  Array<OneD, NekDouble> &outarray,
2343  int BndID)
2344  {
2345  ASSERTL0(false,
2346  "This method is not defined or valid for this class type");
2347  }
2348 
2349  void ExpList::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
2350  {
2351  ASSERTL0(false,
2352  "This method is not defined or valid for this class type");
2353  }
2354 
2355  /**
2356  */
2357  void ExpList::v_FillBndCondFromField()
2358  {
2359  ASSERTL0(false,
2360  "This method is not defined or valid for this class type");
2361  }
2362 
2363  void ExpList::v_LocalToGlobal(void)
2364  {
2365  ASSERTL0(false,
2366  "This method is not defined or valid for this class type");
2367  }
2368 
2369  void ExpList::v_GlobalToLocal(void)
2370  {
2371  ASSERTL0(false,
2372  "This method is not defined or valid for this class type");
2373  }
2374 
2375 
2376  void ExpList::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
2377  Array<OneD, NekDouble> &outarray,
2378  CoeffState coeffstate)
2379  {
2380  v_BwdTrans_IterPerExp(inarray,outarray);
2381  }
2382 
2383  void ExpList::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
2384  Array<OneD, NekDouble> &outarray,
2385  CoeffState coeffstate)
2386  {
2387  v_FwdTrans_IterPerExp(inarray,outarray);
2388  }
2389 
2390  void ExpList::v_IProductWRTBase(
2391  const Array<OneD, const NekDouble> &inarray,
2392  Array<OneD, NekDouble> &outarray,
2393  CoeffState coeffstate)
2394  {
2395  v_IProductWRTBase_IterPerExp(inarray,outarray);
2396  }
2397 
2398  void ExpList::v_GeneralMatrixOp(
2399  const GlobalMatrixKey &gkey,
2400  const Array<OneD,const NekDouble> &inarray,
2401  Array<OneD, NekDouble> &outarray,
2402  CoeffState coeffstate)
2403  {
2404  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
2405  }
2406 
2407  /**
2408  * The operation is evaluated locally by the elemental
2409  * function StdRegions#StdExpansion#GetCoords.
2410  *
2411  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
2412  * will be stored in this array.
2413  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
2414  * will be stored in this array.
2415  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
2416  * will be stored in this array.
2417  */
2418  void ExpList::v_GetCoords(Array<OneD, NekDouble> &coord_0,
2419  Array<OneD, NekDouble> &coord_1,
2420  Array<OneD, NekDouble> &coord_2)
2421  {
2422  int i;
2423  Array<OneD, NekDouble> e_coord_0;
2424  Array<OneD, NekDouble> e_coord_1;
2425  Array<OneD, NekDouble> e_coord_2;
2426 
2427  switch(GetExp(0)->GetCoordim())
2428  {
2429  case 1:
2430  for(i= 0; i < (*m_exp).size(); ++i)
2431  {
2432  e_coord_0 = coord_0 + m_phys_offset[i];
2433  (*m_exp)[i]->GetCoords(e_coord_0);
2434  }
2435  break;
2436  case 2:
2437  ASSERTL0(coord_1.num_elements() != 0,
2438  "output coord_1 is not defined");
2439 
2440  for(i= 0; i < (*m_exp).size(); ++i)
2441  {
2442  e_coord_0 = coord_0 + m_phys_offset[i];
2443  e_coord_1 = coord_1 + m_phys_offset[i];
2444  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1);
2445  }
2446  break;
2447  case 3:
2448  ASSERTL0(coord_1.num_elements() != 0,
2449  "output coord_1 is not defined");
2450  ASSERTL0(coord_2.num_elements() != 0,
2451  "output coord_2 is not defined");
2452 
2453  for(i= 0; i < (*m_exp).size(); ++i)
2454  {
2455  e_coord_0 = coord_0 + m_phys_offset[i];
2456  e_coord_1 = coord_1 + m_phys_offset[i];
2457  e_coord_2 = coord_2 + m_phys_offset[i];
2458  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1,e_coord_2);
2459  }
2460  break;
2461  }
2462  }
2463 
2464  /**
2465  */
2466  void ExpList::v_SetUpPhysNormals()
2467  {
2468  ASSERTL0(false,
2469  "This method is not defined or valid for this class type");
2470  }
2471 
2472  /**
2473  */
2474  void ExpList::v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
2475  Array<OneD,int> &EdgeID)
2476  {
2477  ASSERTL0(false,
2478  "This method is not defined or valid for this class type");
2479  }
2480 
2481  /**
2482  */
2483  void ExpList::v_ReadGlobalOptimizationParameters()
2484  {
2485  ASSERTL0(false,
2486  "This method is not defined or valid for this class type");
2487  }
2488 
2489  /**
2490  */
2491  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>
2492  &ExpList::v_GetBndConditions(void)
2493  {
2494  ASSERTL0(false,
2495  "This method is not defined or valid for this class type");
2496  static Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
2497  result;
2498  return result;
2499  }
2500 
2501  /**
2502  */
2503  Array<OneD,SpatialDomains::BoundaryConditionShPtr> &ExpList::v_UpdateBndConditions()
2504  {
2505  ASSERTL0(false,
2506  "This method is not defined or valid for this class type");
2507  static Array<OneD, SpatialDomains::BoundaryConditionShPtr> result;
2508  return result;
2509  }
2510 
2511  /**
2512  */
2513  void ExpList::v_EvaluateBoundaryConditions(
2514  const NekDouble time,
2515  const std::string varName,
2516  const NekDouble x2_in,
2517  const NekDouble x3_in)
2518  {
2519  ASSERTL0(false,
2520  "This method is not defined or valid for this class type");
2521  }
2522 
2523  /**
2524  */
2525  map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
2526  {
2527  ASSERTL0(false,
2528  "This method is not defined or valid for this class type");
2529  static map<int,RobinBCInfoSharedPtr> result;
2530  return result;
2531  }
2532 
2533  /**
2534  */
2535  void ExpList::v_GetPeriodicEntities(
2536  PeriodicMap &periodicVerts,
2537  PeriodicMap &periodicEdges,
2538  PeriodicMap &periodicFaces)
2539  {
2540  ASSERTL0(false,
2541  "This method is not defined or valid for this class type");
2542  }
2543 
2544  SpatialDomains::BoundaryConditionShPtr ExpList::GetBoundaryCondition(
2545  const SpatialDomains::BoundaryConditionCollection& collection,
2546  unsigned int regionId,
2547  const std::string& variable)
2548  {
2549  SpatialDomains::BoundaryConditionCollection::const_iterator collectionIter = collection.find(regionId);
2550  ASSERTL1(collectionIter != collection.end(), "Unable to locate collection "+boost::lexical_cast<string>(regionId));
2551  const SpatialDomains::BoundaryConditionMapShPtr boundaryConditionMap = (*collectionIter).second;
2552  SpatialDomains::BoundaryConditionMap::const_iterator conditionMapIter = boundaryConditionMap->find(variable);
2553  ASSERTL1(conditionMapIter != boundaryConditionMap->end(), "Unable to locate condition map.");
2554  const SpatialDomains::BoundaryConditionShPtr boundaryCondition = (*conditionMapIter).second;
2555  return boundaryCondition;
2556  }
2557 
2558  ExpListSharedPtr &ExpList::v_GetPlane(int n)
2559  {
2560  ASSERTL0(false,
2561  "This method is not defined or valid for this class type");
2562  return NullExpListSharedPtr;
2563  }
2564  } //end of namespace
2565 } //end of namespace
2566 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2