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, 0.0);
182  m_phys = Array<OneD, NekDouble>(m_npoints, 0.0);
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 dist = 0.0;
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,
1244  locCoords,
1245  tol, resid))
1246  {
1247  w.SetX(gloCoords[0]);
1248  w.SetY(gloCoords[1]);
1249  w.SetZ(gloCoords[2]);
1250 
1251  // Find closest vertex
1252  for (int j = 0; j < (*m_exp)[i]->GetNverts(); ++j) {
1253  v = m_graph->GetVertex(
1254  (*m_exp)[i]->GetGeom()->GetVid(j));
1255  if (j == 0 || dist > v->dist(w))
1256  {
1257  dist = v->dist(w);
1258  }
1259  }
1260  elmtIdDist.push_back(
1261  std::pair<int, NekDouble>(i, dist));
1262  }
1263  }
1264 
1265  // Find nearest element
1266  if (!elmtIdDist.empty())
1267  {
1268  int min_id = elmtIdDist[0].first;
1269  NekDouble min_d = elmtIdDist[0].second;
1270 
1271  for (int i = 1; i < elmtIdDist.size(); ++i)
1272  {
1273  if (elmtIdDist[i].second < min_d) {
1274  min_id = elmtIdDist[i].first;
1275  min_d = elmtIdDist[i].second;
1276  }
1277  }
1278 
1279  // retrieve local coordinate of point
1280  (*m_exp)[min_id]->GetGeom()->GetLocCoords(gloCoords,
1281  locCoords);
1282  return min_id;
1283  }
1284  else
1285  {
1286  return -1;
1287  }
1288  }
1289  // non-embedded mesh (point can only match one element)
1290  else
1291  {
1292  static int start = 0;
1293  int min_id = 0;
1294  NekDouble resid_min = 1e6;
1295  Array<OneD, NekDouble> savLocCoords(locCoords.num_elements());
1296 
1297  // restart search from last found value
1298  for (int i = start; i < (*m_exp).size(); ++i)
1299  {
1300  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
1301  tol, resid))
1302  {
1303  start = i;
1304  return i;
1305  }
1306  else
1307  {
1308  if(resid < resid_min)
1309  {
1310  min_id = i;
1311  resid_min = resid;
1312  Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
1313  }
1314  }
1315  }
1316 
1317  for (int i = 0; i < start; ++i)
1318  {
1319  if ((*m_exp)[i]->GetGeom()->ContainsPoint(gloCoords, locCoords,
1320  tol, resid))
1321  {
1322  start = i;
1323  return i;
1324  }
1325  else
1326  {
1327  if(resid < resid_min)
1328  {
1329  min_id = i;
1330  resid_min = resid;
1331  Vmath::Vcopy(locCoords.num_elements(),savLocCoords,1,locCoords,1);
1332  }
1333  }
1334  }
1335 
1336  std::string msg = "Failed to find point in element to tolerance of "
1337  + boost::lexical_cast<std::string>(resid)
1338  + " using nearest point found";
1339  WARNINGL0(true,msg.c_str());
1340 
1341  if(returnNearestElmt)
1342  {
1343  Vmath::Vcopy(locCoords.num_elements(),locCoords,1,savLocCoords,1);
1344  return min_id;
1345  }
1346  else
1347  {
1348  return -1;
1349  }
1350 
1351  }
1352  }
1353 
1354 
1355  /**
1356  * Configures geometric info, such as tangent direction, on each
1357  * expansion.
1358  * @param graph2D Mesh
1359  */
1360  void ExpList::ApplyGeomInfo()
1361  {
1362 
1363  }
1364 
1365  /**
1366  * Write Tecplot Files Header
1367  * @param outfile Output file name.
1368  * @param var variables names
1369  */
1370  void ExpList::v_WriteTecplotHeader(std::ostream &outfile,
1371  std::string var)
1372  {
1373  int coordim = GetExp(0)->GetCoordim();
1374  char vars[3] = { 'x', 'y', 'z' };
1375 
1376  if (m_expType == e3DH1D)
1377  {
1378  coordim += 1;
1379  }
1380  else if (m_expType == e3DH2D)
1381  {
1382  coordim += 2;
1383  }
1384 
1385  outfile << "Variables = x";
1386  for (int i = 1; i < coordim; ++i)
1387  {
1388  outfile << ", " << vars[i];
1389  }
1390 
1391  if (var.size() > 0)
1392  {
1393  outfile << ", " << var;
1394  }
1395 
1396  outfile << std::endl << std::endl;
1397  }
1398 
1399  /**
1400  * Write Tecplot Files Zone
1401  * @param outfile Output file name.
1402  * @param expansion Expansion that is considered
1403  */
1404  void ExpList::v_WriteTecplotZone(std::ostream &outfile, int expansion)
1405  {
1406  int i, j;
1407  int coordim = GetCoordim(0);
1408  int nPoints = GetTotPoints();
1409  int nBases = (*m_exp)[0]->GetNumBases();
1410  int numBlocks = 0;
1411 
1412  Array<OneD, Array<OneD, NekDouble> > coords(3);
1413 
1414  if (expansion == -1)
1415  {
1416  nPoints = GetTotPoints();
1417 
1418  coords[0] = Array<OneD, NekDouble>(nPoints);
1419  coords[1] = Array<OneD, NekDouble>(nPoints);
1420  coords[2] = Array<OneD, NekDouble>(nPoints);
1421 
1422  GetCoords(coords[0], coords[1], coords[2]);
1423 
1424  for (i = 0; i < m_exp->size(); ++i)
1425  {
1426  int numInt = 1;
1427 
1428  for (j = 0; j < nBases; ++j)
1429  {
1430  numInt *= (*m_exp)[i]->GetNumPoints(j)-1;
1431  }
1432 
1433  numBlocks += numInt;
1434  }
1435  }
1436  else
1437  {
1438  nPoints = (*m_exp)[expansion]->GetTotPoints();
1439 
1440  coords[0] = Array<OneD, NekDouble>(nPoints);
1441  coords[1] = Array<OneD, NekDouble>(nPoints);
1442  coords[2] = Array<OneD, NekDouble>(nPoints);
1443 
1444  (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
1445 
1446  numBlocks = 1;
1447  for (j = 0; j < nBases; ++j)
1448  {
1449  numBlocks *= (*m_exp)[expansion]->GetNumPoints(j)-1;
1450  }
1451  }
1452 
1453  if (m_expType == e3DH1D)
1454  {
1455  nBases += 1;
1456  coordim += 1;
1457  int nPlanes = GetZIDs().num_elements();
1458  NekDouble tmp = numBlocks * (nPlanes-1.0) / nPlanes;
1459  numBlocks = (int)tmp;
1460  }
1461  else if (m_expType == e3DH2D)
1462  {
1463  nBases += 2;
1464  coordim += 1;
1465  }
1466 
1467  outfile << "Zone, N=" << nPoints << ", E="
1468  << numBlocks << ", F=FEBlock" ;
1469 
1470  switch(nBases)
1471  {
1472  case 2:
1473  outfile << ", ET=QUADRILATERAL" << std::endl;
1474  break;
1475  case 3:
1476  outfile << ", ET=BRICK" << std::endl;
1477  break;
1478  default:
1479  ASSERTL0(false,"Not set up for this type of output");
1480  break;
1481  }
1482 
1483  // Write out coordinates
1484  for (j = 0; j < coordim; ++j)
1485  {
1486  for (i = 0; i < nPoints; ++i)
1487  {
1488  outfile << coords[j][i] << " ";
1489  if (i % 1000 == 0 && i)
1490  {
1491  outfile << std::endl;
1492  }
1493  }
1494  outfile << std::endl;
1495  }
1496  }
1497 
1498  void ExpList::v_WriteTecplotConnectivity(std::ostream &outfile,
1499  int expansion)
1500  {
1501  int i,j,k,l;
1502  int nbase = (*m_exp)[0]->GetNumBases();
1503  int cnt = 0;
1504 
1505  boost::shared_ptr<LocalRegions::ExpansionVector> exp = m_exp;
1506 
1507  if (expansion != -1)
1508  {
1509  exp = boost::shared_ptr<LocalRegions::ExpansionVector>(
1510  new LocalRegions::ExpansionVector(1));
1511  (*exp)[0] = (*m_exp)[expansion];
1512  }
1513 
1514  if (nbase == 2)
1515  {
1516  for(i = 0; i < (*exp).size(); ++i)
1517  {
1518  const int np0 = (*exp)[i]->GetNumPoints(0);
1519  const int np1 = (*exp)[i]->GetNumPoints(1);
1520 
1521  for(j = 1; j < np1; ++j)
1522  {
1523  for(k = 1; k < np0; ++k)
1524  {
1525  outfile << cnt + (j-1)*np0 + k << " ";
1526  outfile << cnt + (j-1)*np0 + k+1 << " ";
1527  outfile << cnt + j *np0 + k+1 << " ";
1528  outfile << cnt + j *np0 + k << endl;
1529  }
1530  }
1531 
1532  cnt += np0*np1;
1533  }
1534  }
1535  else if (nbase == 3)
1536  {
1537  for(i = 0; i < (*exp).size(); ++i)
1538  {
1539  const int np0 = (*exp)[i]->GetNumPoints(0);
1540  const int np1 = (*exp)[i]->GetNumPoints(1);
1541  const int np2 = (*exp)[i]->GetNumPoints(2);
1542  const int np01 = np0*np1;
1543 
1544  for(j = 1; j < np2; ++j)
1545  {
1546  for(k = 1; k < np1; ++k)
1547  {
1548  for(l = 1; l < np0; ++l)
1549  {
1550  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l << " ";
1551  outfile << cnt + (j-1)*np01 + (k-1)*np0 + l+1 << " ";
1552  outfile << cnt + (j-1)*np01 + k *np0 + l+1 << " ";
1553  outfile << cnt + (j-1)*np01 + k *np0 + l << " ";
1554  outfile << cnt + j *np01 + (k-1)*np0 + l << " ";
1555  outfile << cnt + j *np01 + (k-1)*np0 + l+1 << " ";
1556  outfile << cnt + j *np01 + k *np0 + l+1 << " ";
1557  outfile << cnt + j *np01 + k *np0 + l << endl;
1558  }
1559  }
1560  }
1561  cnt += np0*np1*np2;
1562  }
1563  }
1564  else
1565  {
1566  ASSERTL0(false,"Not set up for this dimension");
1567  }
1568  }
1569 
1570  /**
1571  * Write Tecplot Files Field
1572  * @param outfile Output file name.
1573  * @param expansion Expansion that is considered
1574  */
1575  void ExpList::v_WriteTecplotField(std::ostream &outfile, int expansion)
1576  {
1577  if (expansion == -1)
1578  {
1579  int totpoints = GetTotPoints();
1580  if(m_physState == false)
1581  {
1582  BwdTrans(m_coeffs,m_phys);
1583  }
1584 
1585  for(int i = 0; i < totpoints; ++i)
1586  {
1587  outfile << m_phys[i] << " ";
1588  if(i % 1000 == 0 && i)
1589  {
1590  outfile << std::endl;
1591  }
1592  }
1593  outfile << std::endl;
1594 
1595  }
1596  else
1597  {
1598  int nPoints = (*m_exp)[expansion]->GetTotPoints();
1599 
1600  for (int i = 0; i < nPoints; ++i)
1601  {
1602  outfile << m_phys[i + m_phys_offset[expansion]] << " ";
1603  }
1604 
1605  outfile << std::endl;
1606  }
1607  }
1608 
1609  void ExpList::WriteVtkHeader(std::ostream &outfile)
1610  {
1611  outfile << "<?xml version=\"1.0\"?>" << endl;
1612  outfile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
1613  << "byte_order=\"LittleEndian\">" << endl;
1614  outfile << " <UnstructuredGrid>" << endl;
1615  }
1616 
1617  void ExpList::WriteVtkFooter(std::ostream &outfile)
1618  {
1619  outfile << " </UnstructuredGrid>" << endl;
1620  outfile << "</VTKFile>" << endl;
1621  }
1622 
1623  void ExpList::v_WriteVtkPieceHeader(std::ostream &outfile, int expansion)
1624  {
1625  ASSERTL0(false, "Routine not implemented for this expansion.");
1626  }
1627 
1628  void ExpList::WriteVtkPieceFooter(std::ostream &outfile, int expansion)
1629  {
1630  outfile << " </PointData>" << endl;
1631  outfile << " </Piece>" << endl;
1632  }
1633 
1634  void ExpList::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
1635  std::string var)
1636  {
1637  int i;
1638  int nq = (*m_exp)[expansion]->GetTotPoints();
1639 
1640  // printing the fields of that zone
1641  outfile << " <DataArray type=\"Float64\" Name=\""
1642  << var << "\">" << endl;
1643  outfile << " ";
1644  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion];
1645  for(i = 0; i < nq; ++i)
1646  {
1647  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
1648  }
1649  outfile << endl;
1650  outfile << " </DataArray>" << endl;
1651  }
1652 
1653  /**
1654  * Given a spectral/hp approximation
1655  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1656  * (which should be contained in #m_phys), this function calculates the
1657  * \f$L_\infty\f$ error of this approximation with respect to an exact
1658  * solution. The local distribution of the quadrature points allows an
1659  * elemental evaluation of this operation through the functions
1660  * StdRegions#StdExpansion#Linf.
1661  *
1662  * The exact solution, also evaluated at the quadrature
1663  * points, should be contained in the variable #m_phys of
1664  * the ExpList object \a Sol.
1665  *
1666  * @param soln A 1D array, containing the discrete
1667  * evaluation of the exact solution at the
1668  * quadrature points in its array #m_phys.
1669  * @return The \f$L_\infty\f$ error of the approximation.
1670  */
1671  NekDouble ExpList::Linf(
1672  const Array<OneD, const NekDouble> &inarray,
1673  const Array<OneD, const NekDouble> &soln)
1674  {
1675  NekDouble err = 0.0;
1676 
1677  if (soln == NullNekDouble1DArray)
1678  {
1679  err = Vmath::Vmax(m_npoints, inarray, 1);
1680  }
1681  else
1682  {
1683  for (int i = 0; i < m_npoints; ++i)
1684  {
1685  err = max(err, abs(inarray[i] - soln[i]));
1686  }
1687  }
1688 
1689  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceMax);
1690 
1691  return err;
1692  }
1693 
1694  /**
1695  * Given a spectral/hp approximation \f$u^{\delta}(\boldsymbol{x})\f$
1696  * evaluated at the quadrature points (which should be contained in
1697  * #m_phys), this function calculates the \f$L_2\f$ error of this
1698  * approximation with respect to an exact solution. The local
1699  * distribution of the quadrature points allows an elemental evaluation
1700  * of this operation through the functions StdRegions#StdExpansion#L2.
1701  *
1702  * The exact solution, also evaluated at the quadrature points, should
1703  * be contained in the variable #m_phys of the ExpList object \a Sol.
1704  *
1705  * @param Sol An ExpList, containing the discrete
1706  * evaluation of the exact solution at the
1707  * quadrature points in its array #m_phys.
1708  * @return The \f$L_2\f$ error of the approximation.
1709  */
1710  NekDouble ExpList::v_L2(
1711  const Array<OneD, const NekDouble> &inarray,
1712  const Array<OneD, const NekDouble> &soln)
1713  {
1714  NekDouble err = 0.0, errl2;
1715  int i;
1716 
1717  if (soln == NullNekDouble1DArray)
1718  {
1719  for (i = 0; i < (*m_exp).size(); ++i)
1720  {
1721  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i]);
1722  err += errl2*errl2;
1723  }
1724  }
1725  else
1726  {
1727  for (i = 0; i < (*m_exp).size(); ++i)
1728  {
1729  errl2 = (*m_exp)[i]->L2(inarray + m_phys_offset[i],
1730  soln + m_phys_offset[i]);
1731  err += errl2*errl2;
1732  }
1733  }
1734 
1735  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1736 
1737  return sqrt(err);
1738  }
1739 
1740  NekDouble ExpList::v_Integral(const Array<OneD, const NekDouble> &inarray)
1741  {
1742  NekDouble err = 0.0;
1743  int i = 0;
1744 
1745  for (i = 0; i < (*m_exp).size(); ++i)
1746  {
1747  err += (*m_exp)[m_offset_elmt_id[i]]->Integral(inarray + m_phys_offset[i]);
1748  }
1749  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1750 
1751  return err;
1752  }
1753 
1754  Array<OneD, const NekDouble> ExpList::v_HomogeneousEnergy (void)
1755  {
1756  ASSERTL0(false,
1757  "This method is not defined or valid for this class type");
1758  Array<OneD, NekDouble> NoEnergy(1,0.0);
1759  return NoEnergy;
1760  }
1761 
1762  LibUtilities::TranspositionSharedPtr ExpList::v_GetTransposition(void)
1763  {
1764  ASSERTL0(false,
1765  "This method is not defined or valid for this class type");
1766  LibUtilities::TranspositionSharedPtr trans;
1767 
1768  return trans;
1769  }
1770 
1771  NekDouble ExpList::v_GetHomoLen(void)
1772  {
1773  ASSERTL0(false,
1774  "This method is not defined or valid for this class type");
1775  NekDouble len = 0.0;
1776  return len;
1777  }
1778 
1779  Array<OneD, const unsigned int> ExpList::v_GetZIDs(void)
1780  {
1781  ASSERTL0(false,
1782  "This method is not defined or valid for this class type");
1783  Array<OneD, unsigned int> NoModes(1);
1784 
1785  return NoModes;
1786  }
1787 
1788  Array<OneD, const unsigned int> ExpList::v_GetYIDs(void)
1789  {
1790  ASSERTL0(false,
1791  "This method is not defined or valid for this class type");
1792  Array<OneD, unsigned int> NoModes(1);
1793 
1794  return NoModes;
1795  }
1796 
1797 
1798  void ExpList::v_PhysInterp1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
1799  {
1800  ASSERTL0(false,
1801  "This method is not defined or valid for this class type");
1802  }
1803 
1804  void ExpList::v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray) {
1805  ASSERTL0(false,
1806  "This method is not defined or valid for this class type");
1807  }
1808 
1809  void ExpList::ExtractFileBCs(
1810  const std::string &fileName,
1811  const std::string &varName,
1812  const boost::shared_ptr<ExpList> locExpList)
1813  {
1814  string varString = fileName.substr(0, fileName.find_last_of("."));
1815  int j, k, len = varString.length();
1816  varString = varString.substr(len-1, len);
1817 
1818  std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
1819  std::vector<std::vector<NekDouble> > FieldData;
1820 
1821  LibUtilities::FieldIO f(m_session->GetComm());
1822  f.Import(fileName, FieldDef, FieldData);
1823 
1824  bool found = false;
1825  for (j = 0; j < FieldDef.size(); ++j)
1826  {
1827  for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
1828  {
1829  if (FieldDef[j]->m_fields[k] == varName)
1830  {
1831  // Copy FieldData into locExpList
1832  locExpList->ExtractDataToCoeffs(
1833  FieldDef[j], FieldData[j],
1834  FieldDef[j]->m_fields[k],
1835  locExpList->UpdateCoeffs());
1836  found = true;
1837  }
1838  }
1839  }
1840 
1841  ASSERTL0(found, "Could not find variable '"+varName+
1842  "' in file boundary condition "+fileName);
1843  locExpList->BwdTrans_IterPerExp(
1844  locExpList->GetCoeffs(),
1845  locExpList->UpdatePhys());
1846  }
1847 
1848  /**
1849  * Given a spectral/hp approximation
1850  * \f$u^{\delta}(\boldsymbol{x})\f$ evaluated at the quadrature points
1851  * (which should be contained in #m_phys), this function calculates the
1852  * \f$H^1_2\f$ error of this approximation with respect to an exact
1853  * solution. The local distribution of the quadrature points allows an
1854  * elemental evaluation of this operation through the functions
1855  * StdRegions#StdExpansion#H1.
1856  *
1857  * The exact solution, also evaluated at the quadrature points, should
1858  * be contained in the variable #m_phys of the ExpList object \a Sol.
1859  *
1860  * @param soln An 1D array, containing the discrete evaluation
1861  * of the exact solution at the quadrature points.
1862  *
1863  * @return The \f$H^1_2\f$ error of the approximation.
1864  */
1865  NekDouble ExpList::H1(
1866  const Array<OneD, const NekDouble> &inarray,
1867  const Array<OneD, const NekDouble> &soln)
1868  {
1869  NekDouble err = 0.0, errh1;
1870  int i;
1871 
1872  for (i = 0; i < (*m_exp).size(); ++i)
1873  {
1874  errh1 = (*m_exp)[i]->H1(inarray + m_phys_offset[i],
1875  soln + m_phys_offset[i]);
1876  err += errh1*errh1;
1877  }
1878 
1879  m_comm->GetRowComm()->AllReduce(err, LibUtilities::ReduceSum);
1880 
1881  return sqrt(err);
1882  }
1883 
1884  void ExpList::GeneralGetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
1885  int NumHomoDir,
1886  Array<OneD, LibUtilities::BasisSharedPtr> &HomoBasis,
1887  std::vector<NekDouble> &HomoLen,
1888  std::vector<unsigned int> &HomoZIDs,
1889  std::vector<unsigned int> &HomoYIDs)
1890  {
1891  int startenum = (int) LibUtilities::eSegment;
1892  int endenum = (int) LibUtilities::eHexahedron;
1893  int s = 0;
1894  LibUtilities::ShapeType shape;
1895 
1896  ASSERTL1(NumHomoDir == HomoBasis.num_elements(),"Homogeneous basis is not the same length as NumHomoDir");
1897  ASSERTL1(NumHomoDir == HomoLen.size(),"Homogeneous length vector is not the same length as NumHomDir");
1898 
1899  // count number of shapes
1900  switch((*m_exp)[0]->GetShapeDimension())
1901  {
1902  case 1:
1903  startenum = (int) LibUtilities::eSegment;
1904  endenum = (int) LibUtilities::eSegment;
1905  break;
1906  case 2:
1907  startenum = (int) LibUtilities::eTriangle;
1908  endenum = (int) LibUtilities::eQuadrilateral;
1909  break;
1910  case 3:
1911  startenum = (int) LibUtilities::eTetrahedron;
1912  endenum = (int) LibUtilities::eHexahedron;
1913  break;
1914  }
1915 
1916  for(s = startenum; s <= endenum; ++s)
1917  {
1918  std::vector<unsigned int> elementIDs;
1919  std::vector<LibUtilities::BasisType> basis;
1920  std::vector<unsigned int> numModes;
1921  std::vector<std::string> fields;
1922 
1923  bool first = true;
1924  bool UniOrder = true;
1925  int n;
1926 
1927  shape = (LibUtilities::ShapeType) s;
1928 
1929  for(int i = 0; i < (*m_exp).size(); ++i)
1930  {
1931  if((*m_exp)[i]->GetGeom()->GetShapeType() == shape)
1932  {
1933  elementIDs.push_back((*m_exp)[i]->GetGeom()->GetGlobalID());
1934  if(first)
1935  {
1936  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
1937  {
1938  basis.push_back((*m_exp)[i]->GetBasis(j)->GetBasisType());
1939  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
1940  }
1941 
1942  // add homogeneous direction details if defined
1943  for(n = 0 ; n < NumHomoDir; ++n)
1944  {
1945  basis.push_back(HomoBasis[n]->GetBasisType());
1946  numModes.push_back(HomoBasis[n]->GetNumModes());
1947  }
1948 
1949  first = false;
1950  }
1951  else
1952  {
1953  ASSERTL0((*m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],"Routine is not set up for multiple bases definitions");
1954 
1955  for(int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
1956  {
1957  numModes.push_back((*m_exp)[i]->GetBasis(j)->GetNumModes());
1958  if(numModes[j] != (*m_exp)[i]->GetBasis(j)->GetNumModes())
1959  {
1960  UniOrder = false;
1961  }
1962  }
1963  // add homogeneous direction details if defined
1964  for(n = 0 ; n < NumHomoDir; ++n)
1965  {
1966  numModes.push_back(HomoBasis[n]->GetNumModes());
1967  }
1968  }
1969  }
1970  }
1971 
1972 
1973  if(elementIDs.size() > 0)
1974  {
1975  LibUtilities::FieldDefinitionsSharedPtr fdef = MemoryManager<LibUtilities::FieldDefinitions>::AllocateSharedPtr(shape, elementIDs, basis, UniOrder, numModes,fields, NumHomoDir, HomoLen, HomoZIDs, HomoYIDs);
1976  fielddef.push_back(fdef);
1977  }
1978  }
1979  }
1980 
1981 
1982  //
1983  // Virtual functions
1984  //
1985  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpList::v_GetFieldDefinitions()
1986  {
1987  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
1988  v_GetFieldDefinitions(returnval);
1989  return returnval;
1990  }
1991 
1992  void ExpList::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
1993  {
1994  GeneralGetFieldDefinitions(fielddef);
1995  }
1996 
1997  //Append the element data listed in elements
1998  //fielddef->m_ElementIDs onto fielddata
1999  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
2000  {
2001  v_AppendFieldData(fielddef,fielddata,m_coeffs);
2002  }
2003 
2004  void ExpList::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
2005  {
2006  int i;
2007  // Determine mapping from element ids to location in
2008  // expansion list
2009  // Determine mapping from element ids to location in
2010  // expansion list
2011  map<int, int> ElmtID_to_ExpID;
2012  for(i = 0; i < (*m_exp).size(); ++i)
2013  {
2014  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2015  }
2016 
2017  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
2018  {
2019  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
2020  int datalen = (*m_exp)[eid]->GetNcoeffs();
2021  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]],&coeffs[m_coeff_offset[eid]]+datalen);
2022  }
2023 
2024  }
2025 
2026  /// Extract the data in fielddata into the coeffs
2027  void ExpList::ExtractDataToCoeffs(
2028  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
2029  std::vector<NekDouble> &fielddata,
2030  std::string &field,
2031  Array<OneD, NekDouble> &coeffs)
2032  {
2033  v_ExtractDataToCoeffs(fielddef,fielddata,field,coeffs);
2034  }
2035 
2036  void ExpList::ExtractCoeffsToCoeffs(const boost::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2037  {
2038  v_ExtractCoeffsToCoeffs(fromExpList,fromCoeffs,toCoeffs);
2039  }
2040 
2041  /**
2042  * @brief Extract data from raw field data into expansion list.
2043  *
2044  * @param fielddef Field definitions.
2045  * @param fielddata Data for associated field.
2046  * @param field Field variable name.
2047  * @param coeffs Resulting coefficient array.
2048  */
2049  void ExpList::v_ExtractDataToCoeffs(
2050  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
2051  std::vector<NekDouble> &fielddata,
2052  std::string &field,
2053  Array<OneD, NekDouble> &coeffs)
2054  {
2055  int i, expId;
2056  int offset = 0;
2057  int modes_offset = 0;
2058  int datalen = fielddata.size()/fielddef->m_fields.size();
2059 
2060  // Find data location according to field definition
2061  for(i = 0; i < fielddef->m_fields.size(); ++i)
2062  {
2063  if(fielddef->m_fields[i] == field)
2064  {
2065  break;
2066  }
2067  offset += datalen;
2068  }
2069 
2070  ASSERTL0(i != fielddef->m_fields.size(),
2071  "Field (" + field + ") not found in file.");
2072 
2073  // Determine mapping from element ids to location in expansion list
2074  map<int, int> elmtToExpId;
2075 
2076  // Loop in reverse order so that in case where using a Homogeneous
2077  // expansion it sets geometry ids to first part of m_exp
2078  // list. Otherwise will set to second (complex) expansion
2079  for(i = (*m_exp).size()-1; i >= 0; --i)
2080  {
2081  elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
2082  }
2083 
2084  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
2085  {
2086  // Reset modes_offset in the case where all expansions of
2087  // the same order.
2088  if (fielddef->m_uniOrder == true)
2089  {
2090  modes_offset = 0;
2091  }
2092 
2093  datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType,
2094  fielddef->m_numModes, modes_offset);
2095 
2096  const int elmtId = fielddef->m_elementIDs[i];
2097  if (elmtToExpId.count(elmtId) == 0)
2098  {
2099  offset += datalen;
2100  continue;
2101  }
2102 
2103  expId = elmtToExpId[elmtId];
2104 
2105  if (datalen == (*m_exp)[expId]->GetNcoeffs())
2106  {
2107  Vmath::Vcopy(datalen, &fielddata[offset], 1,
2108  &coeffs[m_coeff_offset[expId]], 1);
2109  }
2110  else
2111  {
2112  (*m_exp)[expId]->ExtractDataToCoeffs(
2113  &fielddata[offset], fielddef->m_numModes,
2114  modes_offset, &coeffs[m_coeff_offset[expId]]);
2115  }
2116 
2117  offset += datalen;
2118  modes_offset += (*m_exp)[0]->GetNumBases();
2119  }
2120 
2121  return;
2122  }
2123 
2124  void ExpList::v_ExtractCoeffsToCoeffs(const boost::shared_ptr<ExpList> &fromExpList, const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
2125  {
2126  int i;
2127  int offset = 0;
2128 
2129  // check if the same and if so just copy over coeffs
2130  if(fromExpList->GetNcoeffs() == m_ncoeffs)
2131  {
2132  Vmath::Vcopy(m_ncoeffs,fromCoeffs,1,toCoeffs,1);
2133  }
2134  else
2135  {
2136  std::vector<unsigned int> nummodes;
2137  for(i = 0; i < (*m_exp).size(); ++i)
2138  {
2139  int eid = m_offset_elmt_id[i];
2140  for(int j= 0; j < fromExpList->GetExp(eid)->GetNumBases(); ++j)
2141  {
2142  nummodes.push_back(fromExpList->GetExp(eid)->GetBasisNumModes(j));
2143  }
2144 
2145  (*m_exp)[eid]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes,0,
2146  &toCoeffs[m_coeff_offset[eid]]);
2147 
2148  offset += fromExpList->GetExp(eid)->GetNcoeffs();
2149  }
2150  }
2151  }
2152 
2153 
2154  const Array<OneD,const boost::shared_ptr<ExpList> >
2155  &ExpList::v_GetBndCondExpansions(void)
2156  {
2157  ASSERTL0(false,
2158  "This method is not defined or valid for this class type");
2159  static Array<OneD,const boost::shared_ptr<ExpList> > result;
2160  return result;
2161  }
2162 
2163  boost::shared_ptr<ExpList> &ExpList::v_UpdateBndCondExpansion(int i)
2164  {
2165  ASSERTL0(false,
2166  "This method is not defined or valid for this class type");
2167  static boost::shared_ptr<ExpList> result;
2168  return result;
2169  }
2170 
2171  void ExpList::v_Upwind(
2172  const Array<OneD, const Array<OneD, NekDouble> > &Vec,
2173  const Array<OneD, const NekDouble> &Fwd,
2174  const Array<OneD, const NekDouble> &Bwd,
2175  Array<OneD, NekDouble> &Upwind)
2176  {
2177  ASSERTL0(false,
2178  "This method is not defined or valid for this class type");
2179  }
2180 
2181  void ExpList::v_Upwind(
2182  const Array<OneD, const NekDouble> &Vn,
2183  const Array<OneD, const NekDouble> &Fwd,
2184  const Array<OneD, const NekDouble> &Bwd,
2185  Array<OneD, NekDouble> &Upwind)
2186  {
2187  ASSERTL0(false,
2188  "This method is not defined or valid for this class type");
2189  }
2190 
2191  boost::shared_ptr<ExpList> &ExpList::v_GetTrace()
2192  {
2193  ASSERTL0(false,
2194  "This method is not defined or valid for this class type");
2195  static boost::shared_ptr<ExpList> returnVal;
2196  return returnVal;
2197  }
2198 
2199  boost::shared_ptr<AssemblyMapDG> &ExpList::v_GetTraceMap()
2200  {
2201  ASSERTL0(false,
2202  "This method is not defined or valid for this class type");
2203  static boost::shared_ptr<AssemblyMapDG> result;
2204  return result;
2205  }
2206 
2207  const Array<OneD, const int> &ExpList::v_GetTraceBndMap()
2208  {
2209  return GetTraceMap()->GetBndCondTraceToGlobalTraceMap();
2210  }
2211 
2212  void ExpList::v_GetNormals(
2213  Array<OneD, Array<OneD, NekDouble> > &normals)
2214  {
2215  ASSERTL0(false,
2216  "This method is not defined or valid for this class type");
2217  }
2218 
2219  void ExpList::v_AddTraceIntegral(
2220  const Array<OneD, const NekDouble> &Fx,
2221  const Array<OneD, const NekDouble> &Fy,
2222  Array<OneD, NekDouble> &outarray)
2223  {
2224  ASSERTL0(false,
2225  "This method is not defined or valid for this class type");
2226  }
2227 
2228  void ExpList::v_AddTraceIntegral(
2229  const Array<OneD, const NekDouble> &Fn,
2230  Array<OneD, NekDouble> &outarray)
2231  {
2232  ASSERTL0(false,
2233  "This method is not defined or valid for this class type");
2234  }
2235 
2236  void ExpList::v_AddFwdBwdTraceIntegral(
2237  const Array<OneD, const NekDouble> &Fwd,
2238  const Array<OneD, const NekDouble> &Bwd,
2239  Array<OneD, NekDouble> &outarray)
2240  {
2241  ASSERTL0(false,
2242  "This method is not defined or valid for this class type");
2243  }
2244 
2245  void ExpList::v_GetFwdBwdTracePhys(Array<OneD,NekDouble> &Fwd,
2246  Array<OneD,NekDouble> &Bwd)
2247  {
2248  ASSERTL0(false,
2249  "This method is not defined or valid for this class type");
2250  }
2251 
2252  void ExpList::v_GetFwdBwdTracePhys(
2253  const Array<OneD,const NekDouble> &field,
2254  Array<OneD,NekDouble> &Fwd,
2255  Array<OneD,NekDouble> &Bwd)
2256  {
2257  ASSERTL0(false,
2258  "This method is not defined or valid for this class type");
2259  }
2260 
2261  void ExpList::v_ExtractTracePhys(Array<OneD,NekDouble> &outarray)
2262  {
2263  ASSERTL0(false,
2264  "This method is not defined or valid for this class type");
2265  }
2266 
2267  void ExpList::v_ExtractTracePhys(
2268  const Array<OneD, const NekDouble> &inarray,
2269  Array<OneD,NekDouble> &outarray)
2270  {
2271  ASSERTL0(false,
2272  "This method is not defined or valid for this class type");
2273  }
2274 
2275  void ExpList::v_MultiplyByInvMassMatrix(
2276  const Array<OneD,const NekDouble> &inarray,
2277  Array<OneD, NekDouble> &outarray,
2278  CoeffState coeffstate)
2279  {
2280  ASSERTL0(false,
2281  "This method is not defined or valid for this class type");
2282  }
2283 
2284  void ExpList::v_HelmSolve(
2285  const Array<OneD, const NekDouble> &inarray,
2286  Array<OneD, NekDouble> &outarray,
2287  const FlagList &flags,
2288  const StdRegions::ConstFactorMap &factors,
2289  const StdRegions::VarCoeffMap &varcoeff,
2290  const Array<OneD, const NekDouble> &dirForcing)
2291  {
2292  ASSERTL0(false, "HelmSolve not implemented.");
2293  }
2294 
2295  void ExpList::v_LinearAdvectionDiffusionReactionSolve(
2296  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2297  const Array<OneD, const NekDouble> &inarray,
2298  Array<OneD, NekDouble> &outarray,
2299  const NekDouble lambda,
2300  CoeffState coeffstate,
2301  const Array<OneD, const NekDouble>& dirForcing)
2302  {
2303  ASSERTL0(false,
2304  "This method is not defined or valid for this class type");
2305  }
2306 
2307  void ExpList::v_LinearAdvectionReactionSolve(
2308  const Array<OneD, Array<OneD, NekDouble> > &velocity,
2309  const Array<OneD, const NekDouble> &inarray,
2310  Array<OneD, NekDouble> &outarray,
2311  const NekDouble lambda,
2312  CoeffState coeffstate,
2313  const Array<OneD, const NekDouble>& dirForcing)
2314  {
2315  ASSERTL0(false,
2316  "This method is not defined or valid for this class type");
2317  }
2318 
2319  void ExpList::v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
2320  Array<OneD, NekDouble> &outarray,
2321  CoeffState coeffstate,
2322  bool Shuff,
2323  bool UnShuff)
2324  {
2325  ASSERTL0(false,
2326  "This method is not defined or valid for this class type");
2327  }
2328 
2329  void ExpList::v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
2330  Array<OneD, NekDouble> &outarray,
2331  CoeffState coeffstate,
2332  bool Shuff,
2333  bool UnShuff)
2334  {
2335  ASSERTL0(false,
2336  "This method is not defined or valid for this class type");
2337  }
2338 
2339  void ExpList::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,const Array<OneD, NekDouble> &inarray2,Array<OneD, NekDouble> &outarray,CoeffState coeffstate)
2340  {
2341  ASSERTL0(false,
2342  "This method is not defined or valid for this class type");
2343  }
2344 
2345 
2346  void ExpList::v_GetBCValues(Array<OneD, NekDouble> &BndVals,
2347  const Array<OneD, NekDouble> &TotField,
2348  int BndID)
2349  {
2350  ASSERTL0(false,
2351  "This method is not defined or valid for this class type");
2352  }
2353 
2354  void ExpList::v_NormVectorIProductWRTBase(Array<OneD, const NekDouble> &V1,
2355  Array<OneD, const NekDouble> &V2,
2356  Array<OneD, NekDouble> &outarray,
2357  int BndID)
2358  {
2359  ASSERTL0(false,
2360  "This method is not defined or valid for this class type");
2361  }
2362 
2363  void ExpList::v_ImposeDirichletConditions(Array<OneD,NekDouble>& outarray)
2364  {
2365  ASSERTL0(false,
2366  "This method is not defined or valid for this class type");
2367  }
2368 
2369  /**
2370  */
2371  void ExpList::v_FillBndCondFromField()
2372  {
2373  ASSERTL0(false,
2374  "This method is not defined or valid for this class type");
2375  }
2376 
2377  void ExpList::v_LocalToGlobal(void)
2378  {
2379  ASSERTL0(false,
2380  "This method is not defined or valid for this class type");
2381  }
2382 
2383  void ExpList::v_GlobalToLocal(void)
2384  {
2385  ASSERTL0(false,
2386  "This method is not defined or valid for this class type");
2387  }
2388 
2389 
2390  void ExpList::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
2391  Array<OneD, NekDouble> &outarray,
2392  CoeffState coeffstate)
2393  {
2394  v_BwdTrans_IterPerExp(inarray,outarray);
2395  }
2396 
2397  void ExpList::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
2398  Array<OneD, NekDouble> &outarray,
2399  CoeffState coeffstate)
2400  {
2401  v_FwdTrans_IterPerExp(inarray,outarray);
2402  }
2403 
2404  void ExpList::v_IProductWRTBase(
2405  const Array<OneD, const NekDouble> &inarray,
2406  Array<OneD, NekDouble> &outarray,
2407  CoeffState coeffstate)
2408  {
2409  v_IProductWRTBase_IterPerExp(inarray,outarray);
2410  }
2411 
2412  void ExpList::v_GeneralMatrixOp(
2413  const GlobalMatrixKey &gkey,
2414  const Array<OneD,const NekDouble> &inarray,
2415  Array<OneD, NekDouble> &outarray,
2416  CoeffState coeffstate)
2417  {
2418  GeneralMatrixOp_IterPerExp(gkey,inarray,outarray);
2419  }
2420 
2421  /**
2422  * The operation is evaluated locally by the elemental
2423  * function StdRegions#StdExpansion#GetCoords.
2424  *
2425  * @param coord_0 After calculation, the \f$x_1\f$ coordinate
2426  * will be stored in this array.
2427  * @param coord_1 After calculation, the \f$x_2\f$ coordinate
2428  * will be stored in this array.
2429  * @param coord_2 After calculation, the \f$x_3\f$ coordinate
2430  * will be stored in this array.
2431  */
2432  void ExpList::v_GetCoords(Array<OneD, NekDouble> &coord_0,
2433  Array<OneD, NekDouble> &coord_1,
2434  Array<OneD, NekDouble> &coord_2)
2435  {
2436  int i;
2437  Array<OneD, NekDouble> e_coord_0;
2438  Array<OneD, NekDouble> e_coord_1;
2439  Array<OneD, NekDouble> e_coord_2;
2440 
2441  switch(GetExp(0)->GetCoordim())
2442  {
2443  case 1:
2444  for(i= 0; i < (*m_exp).size(); ++i)
2445  {
2446  e_coord_0 = coord_0 + m_phys_offset[i];
2447  (*m_exp)[i]->GetCoords(e_coord_0);
2448  }
2449  break;
2450  case 2:
2451  ASSERTL0(coord_1.num_elements() != 0,
2452  "output coord_1 is not defined");
2453 
2454  for(i= 0; i < (*m_exp).size(); ++i)
2455  {
2456  e_coord_0 = coord_0 + m_phys_offset[i];
2457  e_coord_1 = coord_1 + m_phys_offset[i];
2458  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1);
2459  }
2460  break;
2461  case 3:
2462  ASSERTL0(coord_1.num_elements() != 0,
2463  "output coord_1 is not defined");
2464  ASSERTL0(coord_2.num_elements() != 0,
2465  "output coord_2 is not defined");
2466 
2467  for(i= 0; i < (*m_exp).size(); ++i)
2468  {
2469  e_coord_0 = coord_0 + m_phys_offset[i];
2470  e_coord_1 = coord_1 + m_phys_offset[i];
2471  e_coord_2 = coord_2 + m_phys_offset[i];
2472  (*m_exp)[i]->GetCoords(e_coord_0,e_coord_1,e_coord_2);
2473  }
2474  break;
2475  }
2476  }
2477 
2478  /**
2479  */
2480  void ExpList::v_SetUpPhysNormals()
2481  {
2482  ASSERTL0(false,
2483  "This method is not defined or valid for this class type");
2484  }
2485 
2486  /**
2487  */
2488  void ExpList::v_GetBoundaryToElmtMap(Array<OneD, int> &ElmtID,
2489  Array<OneD,int> &EdgeID)
2490  {
2491  ASSERTL0(false,
2492  "This method is not defined or valid for this class type");
2493  }
2494 
2495  /**
2496  */
2497  void ExpList::v_ReadGlobalOptimizationParameters()
2498  {
2499  ASSERTL0(false,
2500  "This method is not defined or valid for this class type");
2501  }
2502 
2503  /**
2504  */
2505  const Array<OneD,const SpatialDomains::BoundaryConditionShPtr>
2506  &ExpList::v_GetBndConditions(void)
2507  {
2508  ASSERTL0(false,
2509  "This method is not defined or valid for this class type");
2510  static Array<OneD, const SpatialDomains::BoundaryConditionShPtr>
2511  result;
2512  return result;
2513  }
2514 
2515  /**
2516  */
2517  Array<OneD,SpatialDomains::BoundaryConditionShPtr> &ExpList::v_UpdateBndConditions()
2518  {
2519  ASSERTL0(false,
2520  "This method is not defined or valid for this class type");
2521  static Array<OneD, SpatialDomains::BoundaryConditionShPtr> result;
2522  return result;
2523  }
2524 
2525  /**
2526  */
2527  void ExpList::v_EvaluateBoundaryConditions(
2528  const NekDouble time,
2529  const std::string varName,
2530  const NekDouble x2_in,
2531  const NekDouble x3_in)
2532  {
2533  ASSERTL0(false,
2534  "This method is not defined or valid for this class type");
2535  }
2536 
2537  /**
2538  */
2539  map<int, RobinBCInfoSharedPtr> ExpList::v_GetRobinBCInfo(void)
2540  {
2541  ASSERTL0(false,
2542  "This method is not defined or valid for this class type");
2543  static map<int,RobinBCInfoSharedPtr> result;
2544  return result;
2545  }
2546 
2547  /**
2548  */
2549  void ExpList::v_GetPeriodicEntities(
2550  PeriodicMap &periodicVerts,
2551  PeriodicMap &periodicEdges,
2552  PeriodicMap &periodicFaces)
2553  {
2554  ASSERTL0(false,
2555  "This method is not defined or valid for this class type");
2556  }
2557 
2558  SpatialDomains::BoundaryConditionShPtr ExpList::GetBoundaryCondition(
2559  const SpatialDomains::BoundaryConditionCollection& collection,
2560  unsigned int regionId,
2561  const std::string& variable)
2562  {
2563  SpatialDomains::BoundaryConditionCollection::const_iterator collectionIter = collection.find(regionId);
2564  ASSERTL1(collectionIter != collection.end(), "Unable to locate collection "+boost::lexical_cast<string>(regionId));
2565  const SpatialDomains::BoundaryConditionMapShPtr boundaryConditionMap = (*collectionIter).second;
2566  SpatialDomains::BoundaryConditionMap::const_iterator conditionMapIter = boundaryConditionMap->find(variable);
2567  ASSERTL1(conditionMapIter != boundaryConditionMap->end(), "Unable to locate condition map.");
2568  const SpatialDomains::BoundaryConditionShPtr boundaryCondition = (*conditionMapIter).second;
2569  return boundaryCondition;
2570  }
2571 
2572  ExpListSharedPtr &ExpList::v_GetPlane(int n)
2573  {
2574  ASSERTL0(false,
2575  "This method is not defined or valid for this class type");
2576  return NullExpListSharedPtr;
2577  }
2578  } //end of namespace
2579 } //end of namespace
2580 
Generated on Wed Mar 4 2015 16:33:14 for Nektar++ by   doxygen 1.8.1.2