Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ExpListHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpListHomogeneous1D.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: An ExpList which is homogeneous in 1-direction
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 #include <StdRegions/StdSegExp.h>
39 #include <StdRegions/StdPointExp.h>
40 #include <LocalRegions/Expansion.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47  namespace MultiRegions
48  {
49  // Forward declaration for typedefs
50  ExpListHomogeneous1D::ExpListHomogeneous1D():
51  ExpList(),
52  m_homogeneousBasis(LibUtilities::NullBasisSharedPtr),
53  m_lhom(1),
54  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr())
55  {
56  }
57 
60  &pSession,const LibUtilities::BasisKey &HomoBasis,
61  const NekDouble lhom,
62  const bool useFFT,
63  const bool dealiasing):
64  ExpList(pSession),
65  m_useFFT(useFFT),
66  m_lhom(lhom),
67  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
68  m_dealiasing(dealiasing)
69  {
71  "Homogeneous Basis is a null basis");
72 
74 
75  m_StripZcomm = m_session->DefinesSolverInfo("HomoStrip") ?
76  m_comm->GetColumnComm()->GetColumnComm() :
77  m_comm->GetColumnComm();
78 
79  ASSERTL0( m_homogeneousBasis->GetNumPoints() %
80  m_StripZcomm->GetSize() == 0,
81  "HomModesZ should be a multiple of npz.");
82 
83  if ( (m_homogeneousBasis->GetBasisType() !=
85  && (m_homogeneousBasis->GetBasisType() !=
87  {
88  ASSERTL0(
89  (m_homogeneousBasis->GetNumPoints() /
90  m_StripZcomm->GetSize()) % 2 == 0,
91  "HomModesZ/npz should be an even integer.");
92  }
93 
95  ::AllocateSharedPtr(HomoBasis, m_comm,
96  m_StripZcomm);
97 
99  m_homogeneousBasis->GetNumPoints() /
100  m_StripZcomm->GetSize());
101 
102  if(m_useFFT)
103  {
105  "NekFFTW", m_homogeneousBasis->GetNumPoints());
106  }
107 
108  if(m_dealiasing)
109  {
110  if(m_useFFT)
111  {
112  int size = m_homogeneousBasis->GetNumPoints() +
113  m_homogeneousBasis->GetNumPoints() / 2;
114  m_padsize = size + (size % 2);
116  .CreateInstance("NekFFTW", m_padsize);
117  }
118  else
119  {
120  ASSERTL0(false, "Dealiasing available just in combination "
121  "with FFTW");
122  }
123  }
124  }
125 
126 
127  /**
128  * @param In ExpListHomogeneous1D object to copy.
129  */
131  ExpList(In,false),
132  m_transposition(In.m_transposition),
133  m_StripZcomm(In.m_StripZcomm),
134  m_useFFT(In.m_useFFT),
135  m_FFT(In.m_FFT),
136  m_tmpIN(In.m_tmpIN),
137  m_tmpOUT(In.m_tmpOUT),
138  m_homogeneousBasis(In.m_homogeneousBasis),
139  m_lhom(In.m_lhom),
140  m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
141  m_dealiasing(In.m_dealiasing),
142  m_padsize(In.m_padsize)
143  {
144  m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
145  }
146 
148  const std::vector<unsigned int> &eIDs):
149  ExpList(In,eIDs,false),
150  m_transposition(In.m_transposition),
151  m_useFFT(In.m_useFFT),
152  m_FFT(In.m_FFT),
153  m_tmpIN(In.m_tmpIN),
154  m_tmpOUT(In.m_tmpOUT),
155  m_homogeneousBasis(In.m_homogeneousBasis),
156  m_lhom(In.m_lhom),
157  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
158  m_dealiasing(In.m_dealiasing),
159  m_padsize(In.m_padsize)
160  {
161  m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
162  }
163 
164  /**
165  * Destructor
166  */
168  {
169  }
170 
172  Array<OneD, NekDouble> &outarray,
173  CoeffState coeffstate,
174  bool Shuff,
175  bool UnShuff)
176  {
177  // Forwards trans
178  Homogeneous1DTrans(inarray,outarray,true,coeffstate,Shuff,UnShuff);
179  }
180 
182  Array<OneD, NekDouble> &outarray,
183  CoeffState coeffstate,
184  bool Shuff,
185  bool UnShuff)
186  {
187  // Backwards trans
188  Homogeneous1DTrans(inarray,outarray,false,coeffstate,Shuff,UnShuff);
189  }
190 
191  /**
192  * Dealiasing routine
193  *
194  * @param inarray1 First term of the product
195  * @param inarray2 Second term of the product
196  * @param outarray Dealiased product
197  */
199  const Array<OneD, NekDouble> &inarray2,
200  Array<OneD, NekDouble> &outarray,
201  CoeffState coeffstate)
202  {
203  int num_dofs = inarray1.num_elements();
204  int N = m_homogeneousBasis->GetNumPoints();
205 
206  Array<OneD, NekDouble> V1(num_dofs);
207  Array<OneD, NekDouble> V2(num_dofs);
208  Array<OneD, NekDouble> V1V2(num_dofs);
209 
210  if(m_WaveSpace)
211  {
212  V1 = inarray1;
213  V2 = inarray2;
214  }
215  else
216  {
217  HomogeneousFwdTrans(inarray1,V1,coeffstate);
218  HomogeneousFwdTrans(inarray2,V2,coeffstate);
219  }
220 
221  int num_points_per_plane = num_dofs/m_planes.num_elements();
222  int num_proc;
223  if(!m_session->DefinesSolverInfo("HomoStrip"))
224  {
225  num_proc = m_comm->GetColumnComm()->GetSize();
226  }
227  else
228  {
229  num_proc = m_StripZcomm->GetSize();
230  }
231  int num_dfts_per_proc = num_points_per_plane / num_proc
232  + (num_points_per_plane % num_proc > 0);
233 
234  Array<OneD, NekDouble> ShufV1(num_dfts_per_proc*N,0.0);
235  Array<OneD, NekDouble> ShufV2(num_dfts_per_proc*N,0.0);
236  Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc*N,0.0);
237 
238  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
239  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
240  Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize,0.0);
241  Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize,0.0);
242 
243  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
244  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
245 
246  m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXYtoZ);
247  m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXYtoZ);
248 
249  // Looping on the pencils
250  for(int i = 0 ; i < num_dfts_per_proc ; i++)
251  {
252  // Copying the i-th pencil pf lenght N into a bigger
253  // pencil of lenght 2N We are in Fourier space
254  Vmath::Vcopy(N, &(ShufV1[i*N]), 1, &(ShufV1_PAD_coef[0]), 1);
255  Vmath::Vcopy(N, &(ShufV2[i*N]), 1, &(ShufV2_PAD_coef[0]), 1);
256 
257  // Moving to physical space using the padded system
258  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
259  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);
260 
261  // Perfroming the vectors multiplication in physical space on
262  // the padded system
263  Vmath::Vmul(m_padsize, ShufV1_PAD_phys, 1,
264  ShufV2_PAD_phys, 1,
265  ShufV1V2_PAD_phys, 1);
266 
267  // Moving back the result (V1*V2)_phys in Fourier space, padded
268  // system
269  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
270 
271  // Copying the first half of the padded pencil in the full
272  // vector (Fourier space)
273  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1,
274  &(ShufV1V2[i*N]), 1);
275  }
276 
277  // Moving the results to the output
278  if (m_WaveSpace)
279  {
280  m_transposition->Transpose(ShufV1V2, outarray, false,
282  }
283  else
284  {
285  m_transposition->Transpose(ShufV1V2, V1V2, false,
287  HomogeneousBwdTrans(V1V2, outarray, coeffstate);
288  }
289  }
290 
291  /**
292  * Dealiasing routine for dot product
293  *
294  * @param inarray1 First term of the product with dimension ndim
295  * (e.g. U)
296  * @param inarray2 Second term of the product with dimension ndim*nvec
297  * (e.g. grad U)
298  * @param outarray Dealiased product with dimension nvec
299  */
301  const Array<OneD, Array<OneD, NekDouble> > &inarray1,
302  const Array<OneD, Array<OneD, NekDouble> > &inarray2,
303  Array<OneD, Array<OneD, NekDouble> > &outarray,
304  CoeffState coeffstate)
305  {
306  int ndim = inarray1.num_elements();
307  ASSERTL1( inarray2.num_elements() % ndim == 0,
308  "Wrong dimensions for DealiasedDotProd.");
309  int nvec = inarray2.num_elements() / ndim;
310 
311  int num_dofs = inarray1[0].num_elements();
312  int N = m_homogeneousBasis->GetNumPoints();
313 
314  int num_points_per_plane = num_dofs/m_planes.num_elements();
315  int num_proc;
316  if(!m_session->DefinesSolverInfo("HomoStrip"))
317  {
318  num_proc = m_comm->GetColumnComm()->GetSize();
319  }
320  else
321  {
322  num_proc = m_StripZcomm->GetSize();
323  }
324  int num_dfts_per_proc = num_points_per_plane / num_proc
325  + (num_points_per_plane % num_proc > 0);
326 
327  // Get inputs in Fourier space
329  Array<OneD, Array<OneD, NekDouble> > V2(ndim*nvec);
330  if(m_WaveSpace)
331  {
332  for (int i = 0; i < ndim; i++)
333  {
334  V1[i] = inarray1[i];
335  }
336  for (int i = 0; i < ndim*nvec; i++)
337  {
338  V2[i] = inarray2[i];
339  }
340  }
341  else
342  {
343  for (int i = 0; i < ndim; i++)
344  {
345  V1[i] = Array<OneD, NekDouble> (num_dofs);
346  HomogeneousFwdTrans(inarray1[i],V1[i],coeffstate);
347  }
348  for (int i = 0; i < ndim*nvec; i++)
349  {
350  V2[i] = Array<OneD, NekDouble> (num_dofs);
351  HomogeneousFwdTrans(inarray2[i],V2[i],coeffstate);
352  }
353  }
354 
355  // Allocate variables for ffts
356  Array<OneD, Array<OneD, NekDouble> > ShufV1(ndim);
357  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
358  Array<OneD, Array<OneD, NekDouble> > ShufV1_PAD_phys(ndim);
359  for (int i = 0; i < ndim; i++)
360  {
361  ShufV1[i] = Array<OneD, NekDouble>
362  (num_dfts_per_proc*N,0.0);
363  ShufV1_PAD_phys[i] = Array<OneD, NekDouble>
364  (m_padsize,0.0);
365  }
366 
367  Array<OneD, Array<OneD, NekDouble> > ShufV2(ndim*nvec);
368  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
369  Array<OneD, Array<OneD, NekDouble> > ShufV2_PAD_phys(ndim*nvec);
370  for (int i = 0; i < ndim*nvec; i++)
371  {
372  ShufV2[i] = Array<OneD, NekDouble>
373  (num_dfts_per_proc*N,0.0);
374  ShufV2_PAD_phys[i] = Array<OneD, NekDouble>
375  (m_padsize,0.0);
376  }
377 
378  Array<OneD, Array<OneD, NekDouble> > ShufV1V2(nvec);
379  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
380  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
381  for (int i = 0; i < nvec; i++)
382  {
383  ShufV1V2[i] = Array<OneD, NekDouble>
384  (num_dfts_per_proc*N,0.0);
385  }
386 
387  // Transpose inputs
388  for (int i = 0; i < ndim; i++)
389  {
390  m_transposition->Transpose(V1[i], ShufV1[i], false,
392  }
393  for (int i = 0; i < ndim*nvec; i++)
394  {
395  m_transposition->Transpose(V2[i], ShufV2[i], false,
397  }
398 
399  // Looping on the pencils
400  for(int i = 0 ; i < num_dfts_per_proc ; i++)
401  {
402  for (int j = 0; j < ndim; j++)
403  {
404  // Copying the i-th pencil pf lenght N into a bigger
405  // pencil of lenght 1.5N We are in Fourier space
406  Vmath::Vcopy(N, &(ShufV1[j][i*N]), 1,
407  &(ShufV1_PAD_coef[0]), 1);
408  // Moving to physical space using the padded system
409  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys[j]);
410  }
411  for (int j = 0; j < ndim*nvec; j++)
412  {
413  Vmath::Vcopy(N, &(ShufV2[j][i*N]), 1,
414  &(ShufV2_PAD_coef[0]), 1);
415  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys[j]);
416  }
417 
418  // Performing the vectors multiplication in physical space on
419  // the padded system
420  for (int j = 0; j < nvec; j++)
421  {
422  Vmath::Zero(m_padsize, ShufV1V2_PAD_phys, 1);
423  for (int k = 0; k < ndim; k++)
424  {
425  Vmath::Vvtvp(m_padsize, ShufV1_PAD_phys[k], 1,
426  ShufV2_PAD_phys[j*ndim+k], 1,
427  ShufV1V2_PAD_phys, 1,
428  ShufV1V2_PAD_phys, 1);
429  }
430  // Moving back the result (V1*V2)_phys in Fourier space,
431  // padded system
432  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
433  // Copying the first half of the padded pencil in the full
434  // vector (Fourier space)
435  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1,
436  &(ShufV1V2[j][i*N]), 1);
437  }
438  }
439 
440  // Moving the results to the output
441  if (m_WaveSpace)
442  {
443  for (int j = 0; j < nvec; j++)
444  {
445  m_transposition->Transpose(ShufV1V2[j], outarray[j],
446  false,
448  }
449  }
450  else
451  {
452  Array<OneD, NekDouble> V1V2(num_dofs);
453  for (int j = 0; j < nvec; j++)
454  {
455  m_transposition->Transpose(ShufV1V2[j], V1V2, false,
457  HomogeneousBwdTrans(V1V2, outarray[j], coeffstate);
458  }
459  }
460  }
461 
462  /**
463  * Forward transform
464  */
466  {
467  int cnt = 0, cnt1 = 0;
468  Array<OneD, NekDouble> tmparray;
469 
470  for(int n = 0; n < m_planes.num_elements(); ++n)
471  {
472  m_planes[n]->FwdTrans(inarray+cnt, tmparray = outarray + cnt1,
473  coeffstate);
474  cnt += m_planes[n]->GetTotPoints();
475 
476  cnt1 += m_planes[n]->GetNcoeffs(); // need to skip ncoeffs
477  }
478  if(!m_WaveSpace)
479  {
480  HomogeneousFwdTrans(outarray,outarray,coeffstate);
481  }
482  }
483 
484  /**
485  * Forward transform element by element
486  */
488  {
489  int cnt = 0, cnt1 = 0;
490  Array<OneD, NekDouble> tmparray;
491 
492  //spectral element FwdTrans plane by plane
493  for(int n = 0; n < m_planes.num_elements(); ++n)
494  {
495  m_planes[n]->FwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
496  cnt += m_planes[n]->GetTotPoints();
497  cnt1 += m_planes[n]->GetNcoeffs();
498  }
499  if(!m_WaveSpace)
500  {
501  HomogeneousFwdTrans(outarray,outarray);
502  }
503  }
504 
505  /**
506  * Backward transform
507  */
509  {
510  int cnt = 0, cnt1 = 0;
511  Array<OneD, NekDouble> tmparray;
512 
513  for(int n = 0; n < m_planes.num_elements(); ++n)
514  {
515  m_planes[n]->BwdTrans(inarray+cnt, tmparray = outarray + cnt1,
516  coeffstate);
517  cnt += m_planes[n]->GetNcoeffs();
518  cnt1 += m_planes[n]->GetTotPoints();
519  }
520  if(!m_WaveSpace)
521  {
522  HomogeneousBwdTrans(outarray,outarray);
523  }
524  }
525 
526  /**
527  * Backward transform element by element
528  */
530  {
531  int cnt = 0, cnt1 = 0;
532  Array<OneD, NekDouble> tmparray;
533 
534  for(int n = 0; n < m_planes.num_elements(); ++n)
535  {
536  m_planes[n]->BwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
537 
538  cnt += m_planes[n]->GetNcoeffs();
539  cnt1 += m_planes[n]->GetTotPoints();
540  }
541  if(!m_WaveSpace)
542  {
543  HomogeneousBwdTrans(outarray,outarray);
544  }
545  }
546 
547  /**
548  * Inner product
549  */
551  {
552  int cnt = 0, cnt1 = 0;
553  Array<OneD, NekDouble> tmparray, tmpIn;
554 
555  if(m_WaveSpace)
556  {
557  tmpIn = inarray;
558  }
559  else
560  {
561  tmpIn = Array<OneD, NekDouble> (inarray.num_elements(), 0.0);
562  HomogeneousFwdTrans(inarray,tmpIn);
563  }
564 
565  for(int n = 0; n < m_planes.num_elements(); ++n)
566  {
567  m_planes[n]->IProductWRTBase(tmpIn+cnt, tmparray = outarray + cnt1,coeffstate);
568 
569  cnt1 += m_planes[n]->GetNcoeffs();
570  cnt += m_planes[n]->GetTotPoints();
571  }
572  }
573 
574  /**
575  * Inner product element by element
576  */
578  {
579  int cnt = 0, cnt1 = 0;
580  Array<OneD, NekDouble> tmparray, tmpIn;
581 
582  if(m_WaveSpace)
583  {
584  tmpIn = inarray;
585  }
586  else
587  {
588  tmpIn = Array<OneD, NekDouble> (inarray.num_elements(), 0.0);
589  HomogeneousFwdTrans(inarray,tmpIn);
590  }
591 
592  for(int n = 0; n < m_planes.num_elements(); ++n)
593  {
594  m_planes[n]->IProductWRTBase_IterPerExp(tmpIn+cnt, tmparray = outarray + cnt1);
595 
596  cnt1 += m_planes[n]->GetNcoeffs();
597  cnt += m_planes[n]->GetTotPoints();
598  }
599  }
600 
601  /**
602  * Homogeneous transform Bwd/Fwd (MVM and FFT)
603  */
605  bool IsForwards,
606  CoeffState coeffstate,
607  bool Shuff,
608  bool UnShuff)
609  {
610  int num_dofs;
611 
612  if(IsForwards)
613  {
614  num_dofs = inarray.num_elements();
615  }
616  else
617  {
618  num_dofs = outarray.num_elements();
619  }
620 
621  if(m_useFFT)
622  {
623  int num_points_per_plane = num_dofs/m_planes.num_elements();
624  int num_dfts_per_proc;
625  if(!m_session->DefinesSolverInfo("HomoStrip"))
626  {
627  int nP = m_comm->GetColumnComm()->GetSize();
628  num_dfts_per_proc = num_points_per_plane / nP
629  + (num_points_per_plane % nP > 0);
630  }
631  else
632  {
633  int nP = m_StripZcomm->GetSize();
634  num_dfts_per_proc = num_points_per_plane / nP
635  + (num_points_per_plane % nP > 0);
636  }
637 
638  Array<OneD, NekDouble> fft_in (num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
639  Array<OneD, NekDouble> fft_out(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
640 
641  if(Shuff)
642  {
643  m_transposition->Transpose(inarray,fft_in,false,LibUtilities::eXYtoZ);
644  }
645  else
646  {
647  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
648  inarray,1,fft_in,1);
649  }
650 
651  if(IsForwards)
652  {
653  for(int i = 0 ; i < num_dfts_per_proc ; i++)
654  {
655  m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
656  }
657  }
658  else
659  {
660  for(int i = 0 ; i < num_dfts_per_proc ; i++)
661  {
662  m_FFT->FFTBwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
663  }
664  }
665 
666  if(UnShuff)
667  {
668  m_transposition->Transpose(fft_out,outarray,false,LibUtilities::eZtoXY);
669  }
670  else
671  {
672  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
673  fft_out,1,outarray,1);
674  }
675  }
676  else
677  {
678  DNekBlkMatSharedPtr blkmat;
679 
680  if(num_dofs == m_npoints) //transform phys space
681  {
682  if(IsForwards)
683  {
685  }
686  else
687  {
689  }
690  }
691  else
692  {
693  if(IsForwards)
694  {
696  }
697  else
698  {
700  }
701  }
702 
703  int nrows = blkmat->GetRows();
704  int ncols = blkmat->GetColumns();
705 
706  Array<OneD, NekDouble> sortedinarray(ncols,0.0);
707  Array<OneD, NekDouble> sortedoutarray(nrows,0.0);
708 
709  if(Shuff)
710  {
711  m_transposition->Transpose(inarray,sortedinarray,!IsForwards,LibUtilities::eXYtoZ);
712  }
713  else
714  {
715  Vmath::Vcopy(ncols,inarray,1,sortedinarray,1);
716  }
717 
718  // Create NekVectors from the given data arrays
719  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
720  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
721 
722  // Perform matrix-vector multiply.
723  out = (*blkmat)*in;
724 
725  if(UnShuff)
726  {
727  m_transposition->Transpose(sortedoutarray,outarray,IsForwards,LibUtilities::eZtoXY);
728  }
729  else
730  {
731  Vmath::Vcopy(nrows,sortedoutarray,1,outarray,1);
732  }
733 
734  }
735  }
736 
738  {
739  Homo1DBlockMatrixMap::iterator matrixIter = m_homogeneous1DBlockMat->find(mattype);
740 
741  if(matrixIter == m_homogeneous1DBlockMat->end())
742  {
743  return ((*m_homogeneous1DBlockMat)[mattype] =
744  GenHomogeneous1DBlockMatrix(mattype,coeffstate));
745  }
746  else
747  {
748  return matrixIter->second;
749  }
750  }
751 
752 
754  {
755  DNekMatSharedPtr loc_mat;
756  DNekBlkMatSharedPtr BlkMatrix;
757  int n_exp = 0;
758  int num_trans_per_proc = 0;
759 
760  if((mattype == eForwardsCoeffSpace1D)
761  ||(mattype == eBackwardsCoeffSpace1D)) // will operate on m_coeffs
762  {
763  n_exp = m_planes[0]->GetNcoeffs();
764  }
765  else
766  {
767  n_exp = m_planes[0]->GetTotPoints(); // will operatore on m_phys
768  }
769 
770  num_trans_per_proc = n_exp/m_comm->GetColumnComm()->GetSize() + (n_exp%m_comm->GetColumnComm()->GetSize() > 0);
771 
772  Array<OneD,unsigned int> nrows(num_trans_per_proc);
773  Array<OneD,unsigned int> ncols(num_trans_per_proc);
774 
775  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
776  {
777  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
778  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
779  }
780  else
781  {
782  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
783  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
784  }
785 
786  MatrixStorage blkmatStorage = eDIAGONAL;
787  BlkMatrix = MemoryManager<DNekBlkMat>
788  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
789 
790  //Half Mode
792  {
793  StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
794 
795  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
796  {
798  StdPoint.DetShapeType(),
799  StdPoint);
800 
801  loc_mat = StdPoint.GetStdMatrix(matkey);
802  }
803  else
804  {
806  StdPoint.DetShapeType(),
807  StdPoint);
808 
809  loc_mat = StdPoint.GetStdMatrix(matkey);
810  }
811  }
812  //other cases
813  else
814  {
815  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
816 
817  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
818  {
820  StdSeg.DetShapeType(),
821  StdSeg);
822 
823  loc_mat = StdSeg.GetStdMatrix(matkey);
824  }
825  else
826  {
828  StdSeg.DetShapeType(),
829  StdSeg);
830 
831  loc_mat = StdSeg.GetStdMatrix(matkey);
832  }
833 
834  }
835 
836  // set up array of block matrices.
837  for(int i = 0; i < num_trans_per_proc; ++i)
838  {
839  BlkMatrix->SetBlock(i,i,loc_mat);
840  }
841 
842  return BlkMatrix;
843  }
844 
845  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous1D::v_GetFieldDefinitions()
846  {
847  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
848 
849  // Set up Homogeneous length details.
851 
852  std::vector<NekDouble> HomoLen;
853  HomoLen.push_back(m_lhom);
854 
855  std::vector<unsigned int> StripsIDs;
856 
857  bool strips;
858  m_session->MatchSolverInfo("HomoStrip","True",strips,false);
859  if (strips)
860  {
861  StripsIDs.push_back(m_transposition->GetStripID());
862  }
863 
864  std::vector<unsigned int> PlanesIDs;
865  int IDoffset = 0;
866 
867  // introduce a 2 plane offset for single mode case so can
868  // be post-processed or used in MultiMode expansion.
870  {
871  IDoffset = 2;
872  }
873 
874  for(int i = 0; i < m_planes.num_elements(); i++)
875  {
876  PlanesIDs.push_back(m_transposition->GetPlaneID(i)+IDoffset);
877  }
878 
879  m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, HomoBasis,
880  HomoLen, strips, StripsIDs, PlanesIDs);
881  return returnval;
882  }
883 
884  void ExpListHomogeneous1D::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
885  {
886  // Set up Homogeneous length details.
888 
889  std::vector<NekDouble> HomoLen;
890  HomoLen.push_back(m_lhom);
891 
892  std::vector<unsigned int> StripsIDs;
893 
894  bool strips;
895  m_session->MatchSolverInfo("HomoStrip","True",strips,false);
896  if (strips)
897  {
898  StripsIDs.push_back(m_transposition->GetStripID());
899  }
900 
901  std::vector<unsigned int> PlanesIDs;
902  int IDoffset = 0;
903 
905  {
906  IDoffset = 2;
907  }
908 
909  for(int i = 0; i < m_planes.num_elements(); i++)
910  {
911  PlanesIDs.push_back(m_transposition->GetPlaneID(i)+IDoffset);
912  }
913 
914  // enforce NumHomoDir == 1 by direct call
915  m_planes[0]->GeneralGetFieldDefinitions(fielddef, 1, HomoBasis,
916  HomoLen, strips, StripsIDs, PlanesIDs);
917  }
918 
919 
920  /** This routine appends the data from the expansion list into
921  the output format where each element is given by looping
922  over its Fourier modes where as data in the expandion is
923  stored with all consecutive elements and then the Fourier
924  modes
925  */
927  {
928  int i,n;
929  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
930 
931  // Determine mapping from element ids to location in
932  // expansion list
933  if (m_elmtToExpId.size() == 0)
934  {
935  for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
936  {
937  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
938  }
939  }
940 
941  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
942  {
943  int eid = m_elmtToExpId[fielddef->m_elementIDs[i]];
944  int datalen = (*m_exp)[eid]->GetNcoeffs();
945 
946  for(n = 0; n < m_planes.num_elements(); ++n)
947  {
948  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane],&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane]+datalen);
949  }
950  }
951  }
952 
953  void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
954  {
955  v_AppendFieldData(fielddef,fielddata,m_coeffs);
956  }
957 
958  //Extract the data in fielddata into the m_coeff list
961  std::vector<NekDouble> &fielddata,
962  std::string &field,
963  Array<OneD, NekDouble> &coeffs)
964  {
965  int i,n;
966  int offset = 0;
967  int nzmodes = 1;
968  int datalen = fielddata.size()/fielddef->m_fields.size();
969  std::vector<unsigned int> fieldDefHomoZids;
970 
971 
972  // Find data location according to field definition
973  for(i = 0; i < fielddef->m_fields.size(); ++i)
974  {
975  if(fielddef->m_fields[i] == field)
976  {
977  break;
978  }
979  offset += datalen;
980  }
981 
982  if(i == fielddef->m_fields.size())
983  {
984  cout << "Field "<< field<< "not found in data file. " << endl;
985  }
986  else
987  {
988 
989  int modes_offset = 0;
990  int planes_offset = 0;
991  Array<OneD, NekDouble> coeff_tmp;
993 
994  // Build map of plane IDs lying on this processor and determine
995  // mapping from element ids to location in expansion list.
996  if (m_zIdToPlane.size() == 0)
997  {
998  for (i = 0; i < m_planes.num_elements(); ++i)
999  {
1000  m_zIdToPlane[m_transposition->GetPlaneID(i)] = i;
1001  }
1002 
1003  for (i = 0; i < m_planes[0]->GetExpSize(); ++i)
1004  {
1005  m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
1006  }
1007  }
1008 
1009  if(fielddef->m_numHomogeneousDir)
1010  {
1011  nzmodes = fielddef->m_homogeneousZIDs.size();
1012  fieldDefHomoZids = fielddef->m_homogeneousZIDs;
1013  }
1014  else // input file is 2D and so set nzmodes to 1
1015  {
1016  nzmodes = 1;
1017  fieldDefHomoZids.push_back(0);
1018  }
1019 
1020  // calculate number of modes in the current partition
1021  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
1022 
1023  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
1024  {
1025  if(fielddef->m_uniOrder == true) // reset modes_offset to zero
1026  {
1027  modes_offset = 0;
1028  }
1029 
1030  int datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType,
1031  fielddef->m_numModes,
1032  modes_offset);
1033 
1034  it = m_elmtToExpId.find(fielddef->m_elementIDs[i]);
1035 
1036  // ensure element is on this partition for parallel case.
1037  if(it == m_elmtToExpId.end())
1038  {
1039  // increase offset for correct FieldData access
1040  offset += datalen*nzmodes;
1041  modes_offset += (*m_exp)[0]->GetNumBases() +
1042  fielddef->m_numHomogeneousDir;
1043  continue;
1044  }
1045 
1046  int eid = it->second;
1047 
1048 
1049  for(n = 0; n < nzmodes; ++n, offset += datalen)
1050  {
1051 
1052  it = m_zIdToPlane.find(fieldDefHomoZids[n]);
1053 
1054  // Check to make sure this mode number lies in this field.
1055  if (it == m_zIdToPlane.end())
1056  {
1057  continue;
1058  }
1059 
1060  planes_offset = it->second;
1061  if(datalen == (*m_exp)[eid]->GetNcoeffs())
1062  {
1063  Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid]+planes_offset*ncoeffs_per_plane],1);
1064  }
1065  else // unpack data to new order
1066  {
1067  (*m_exp)[eid]->ExtractDataToCoeffs(&fielddata[offset], fielddef->m_numModes,modes_offset,&coeffs[m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane], fielddef->m_basis);
1068  }
1069  }
1070  modes_offset += (*m_exp)[0]->GetNumBases() + fielddef->m_numHomogeneousDir;
1071  }
1072  }
1073  }
1074 
1075  //Extract the data in fielddata into the m_coeff list
1077  const boost::shared_ptr<ExpList> &fromExpList,const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
1078  {
1079  int i;
1080  int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
1081  int toNcoeffs_per_plane = m_planes[0]->GetNcoeffs();
1082  Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp;
1083 
1084  for(i = 0; i < m_planes.num_elements(); ++i)
1085  {
1086  m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + toNcoeffs_per_plane*i);
1087  }
1088  }
1089 
1090  void ExpListHomogeneous1D::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
1091  std::string var)
1092  {
1093  // If there is only one plane (e.g. HalfMode), we write a 2D plane.
1094  if (m_planes.num_elements() == 1)
1095  {
1096  m_planes[0]->WriteVtkPieceData(outfile, expansion, var);
1097  return;
1098  }
1099 
1100  int i;
1101  int nq = (*m_exp)[expansion]->GetTotPoints();
1102  int npoints_per_plane = m_planes[0]->GetTotPoints();
1103 
1104  // If we are using Fourier points, output extra plane to fill domain
1105  int outputExtraPlane = 0;
1106  Array<OneD, NekDouble> extraPlane;
1107  if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
1108  && m_homogeneousBasis->GetPointsType() ==
1110  {
1111  outputExtraPlane = 1;
1112  // Get extra plane data
1113  if (m_StripZcomm->GetSize() == 1)
1114  {
1115  extraPlane = m_phys + m_phys_offset[expansion];
1116  }
1117  else
1118  {
1119  // Determine to and from rank for communication
1120  int size = m_StripZcomm->GetSize();
1121  int rank = m_StripZcomm->GetRank();
1122  int fromRank = (rank+1) % size;
1123  int toRank = (rank == 0) ? size-1 : rank-1;
1124  // Communicate using SendRecv
1125  extraPlane = Array<OneD, NekDouble>(nq);
1126  Array<OneD, NekDouble> send (nq,
1127  m_phys + m_phys_offset[expansion]);
1128  m_StripZcomm->SendRecv(toRank, send,
1129  fromRank, extraPlane);
1130  }
1131  }
1132 
1133  // printing the fields of that zone
1134  outfile << " <DataArray type=\"Float64\" Name=\""
1135  << var << "\">" << endl;
1136  outfile << " ";
1137  for (int n = 0; n < m_planes.num_elements(); ++n)
1138  {
1139  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_plane;
1140  for(i = 0; i < nq; ++i)
1141  {
1142  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
1143  }
1144  }
1145  if (outputExtraPlane)
1146  {
1147  for(i = 0; i < nq; ++i)
1148  {
1149  outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol ?
1150  0 : extraPlane[i]) << " ";
1151  }
1152  }
1153  outfile << endl;
1154  outfile << " </DataArray>" << endl;
1155  }
1156 
1158  {
1159  int cnt,cnt1;
1160  Array<OneD, NekDouble> tmparray;
1161  cnt = m_planes[0]->GetTotPoints();
1162  cnt1 = m_planes[0]->Get1DScaledTotPoints(scale);
1163 
1164  ASSERTL1(m_planes.num_elements()*cnt1 <= outarray.num_elements(),"size of outarray does not match internal estimage");
1165 
1166 
1167  for(int i = 0; i < m_planes.num_elements(); i++)
1168  {
1169 
1170  m_planes[i]->PhysInterp1DScaled(scale,inarray+i*cnt,
1171  tmparray = outarray+i*cnt1);
1172  }
1173  }
1174 
1175 
1177  {
1178  int cnt,cnt1;
1179  Array<OneD, NekDouble> tmparray;
1180  cnt = m_planes[0]->Get1DScaledTotPoints(scale);
1181  cnt1 = m_planes[0]->GetTotPoints();
1182 
1183  ASSERTL1(m_planes.num_elements()*cnt <= inarray.num_elements(),"size of outarray does not match internal estimage");
1184 
1185 
1186  for(int i = 0; i < m_planes.num_elements(); i++)
1187  {
1188  m_planes[i]->PhysGalerkinProjection1DScaled(scale,inarray+i*cnt,
1189  tmparray = outarray+i*cnt1);
1190  }
1191 
1192  }
1194  Array<OneD, NekDouble> &out_d0,
1195  Array<OneD, NekDouble> &out_d1,
1196  Array<OneD, NekDouble> &out_d2)
1197  {
1198  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
1199  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
1200 
1201  Array<OneD, NekDouble> temparray(nT_pts);
1202  Array<OneD, NekDouble> outarray(nT_pts);
1205  Array<OneD, NekDouble> tmp3;
1206 
1207  for(int i = 0; i < m_planes.num_elements(); i++)
1208  {
1209  m_planes[i]->PhysDeriv(inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
1210  }
1211 
1212  if(out_d2 != NullNekDouble1DArray)
1213  {
1216  {
1217  if(m_WaveSpace)
1218  {
1219  temparray = inarray;
1220  }
1221  else
1222  {
1223  HomogeneousFwdTrans(inarray,temparray);
1224  }
1225 
1226  NekDouble sign = -1.0;
1227  NekDouble beta;
1228 
1229  //Half Mode
1231  {
1232  beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
1233 
1234  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1235  }
1236  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
1237  {
1238  beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
1239 
1240  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1241  }
1242 
1243  //Fully complex
1244  else
1245  {
1246  for(int i = 0; i < m_planes.num_elements(); i++)
1247  {
1248  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
1249 
1250  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
1251 
1252  sign = -1.0*sign;
1253  }
1254  }
1255 
1256  if(m_WaveSpace)
1257  {
1258  out_d2 = outarray;
1259  }
1260  else
1261  {
1262  HomogeneousBwdTrans(outarray,out_d2);
1263  }
1264  }
1265  else
1266  {
1267  if(!m_session->DefinesSolverInfo("HomoStrip"))
1268  {
1269  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
1270  "Parallelisation in the homogeneous direction "
1271  "implemented just for Fourier basis");
1272  }
1273  else
1274  {
1275  ASSERTL0(m_StripZcomm->GetSize() == 1,
1276  "Parallelisation in the homogeneous direction "
1277  "implemented just for Fourier basis");
1278  }
1279 
1280  if(m_WaveSpace)
1281  {
1282  ASSERTL0(false, "Semi-phyisical time-stepping not "
1283  "implemented yet for non-Fourier "
1284  "basis");
1285  }
1286  else
1287  {
1288  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1289 
1290  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
1291 
1292  for(int i = 0; i < nP_pts; i++)
1293  {
1294  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
1295  }
1296 
1297  m_transposition->Transpose(outarray,out_d2,false,LibUtilities::eZtoXY);
1298 
1299  Vmath::Smul(nT_pts,2.0/m_lhom,out_d2,1,out_d2,1);
1300  }
1301  }
1302  }
1303  }
1304 
1307 
1308  {
1309  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
1310  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
1311 
1312  int dir= (int)edir;
1313 
1314  Array<OneD, NekDouble> temparray(nT_pts);
1315  Array<OneD, NekDouble> outarray(nT_pts);
1318 
1319  if (dir < 2)
1320  {
1321  for(int i=0; i < m_planes.num_elements(); i++)
1322  {
1323  m_planes[i]->PhysDeriv(edir, inarray + i*nP_pts ,tmp2 = out_d + i*nP_pts);
1324  }
1325  }
1326  else
1327  {
1330  {
1331  if(m_WaveSpace)
1332  {
1333  temparray = inarray;
1334  }
1335  else
1336  {
1337  HomogeneousFwdTrans(inarray,temparray);
1338  }
1339 
1340  NekDouble sign = -1.0;
1341  NekDouble beta;
1342 
1343  //HalfMode
1345  {
1346  beta = 2*sign*M_PI*(m_transposition->GetK(0))/m_lhom;
1347 
1348  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1349  }
1350  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
1351  {
1352  beta = -2*sign*M_PI*(m_transposition->GetK(0))/m_lhom;
1353 
1354  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
1355  }
1356  //Fully complex
1357  else
1358  {
1359  for(int i = 0; i < m_planes.num_elements(); i++)
1360  {
1361  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
1362 
1363  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
1364 
1365  sign = -1.0*sign;
1366  }
1367  }
1368  if(m_WaveSpace)
1369  {
1370  out_d = outarray;
1371  }
1372  else
1373  {
1374  HomogeneousBwdTrans(outarray,out_d);
1375  }
1376  }
1377  else
1378  {
1379  if(!m_session->DefinesSolverInfo("HomoStrip"))
1380  {
1381  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
1382  "Parallelisation in the homogeneous direction "
1383  "implemented just for Fourier basis");
1384  }
1385  else
1386  {
1387  ASSERTL0(m_StripZcomm->GetSize() == 1,
1388  "Parallelisation in the homogeneous direction "
1389  "implemented just for Fourier basis");
1390  }
1391 
1392  if(m_WaveSpace)
1393  {
1394  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
1395  }
1396  else
1397  {
1398  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1399 
1400  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
1401 
1402  for(int i = 0; i < nP_pts; i++)
1403  {
1404  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
1405  }
1406 
1407  m_transposition->Transpose(outarray,out_d,false,LibUtilities::eZtoXY);
1408 
1409  Vmath::Smul(nT_pts,2.0/m_lhom,out_d,1,out_d,1);
1410  }
1411  }
1412  }
1413  }
1414 
1416  Array<OneD, NekDouble> &out_d0,
1417  Array<OneD, NekDouble> &out_d1,
1418  Array<OneD, NekDouble> &out_d2)
1419 
1420  {
1421  v_PhysDeriv(inarray,out_d0,out_d1,out_d2);
1422  }
1423 
1425  const Array<OneD, const NekDouble> &inarray,
1426  Array<OneD, NekDouble> &out_d)
1427  {
1428  v_PhysDeriv(edir,inarray,out_d);
1429  }
1430 
1432  {
1433  return m_transposition;
1434  }
1435 
1437  {
1438  return m_lhom;
1439  }
1440 
1442  {
1443  return m_transposition->GetPlanesIDs();
1444  }
1445  } //end of namespace
1446 } //end of namespace
1447 
1448 
1449 /**
1450 * $Log: v $
1451 *
1452 **/
1453 
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
DNekBlkMatSharedPtr GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate=eLocal) const
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
boost::shared_ptr< Transposition > TranspositionSharedPtr
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static Array< OneD, NekDouble > NullNekDouble1DArray
boost::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
Definition: ExpList.h:1062
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:27
LibUtilities::TranspositionSharedPtr m_transposition
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual void v_ExtractCoeffsToCoeffs(const boost::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
boost::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:181
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
LibUtilities::NektarFFTSharedPtr m_FFT_deal
static BasisSharedPtr NullBasisSharedPtr
Definition: Basis.h:358
STL namespace.
NekDouble m_lhom
Width of homogeneous direction.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1015
Fourier Expansion .
Definition: BasisType.h:52
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:998
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset)
Definition: ShapeType.hpp:312
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
BasisManagerT & BasisManager(void)
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1047
LibUtilities::BasisSharedPtr m_homogeneousBasis
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
boost::unordered_map< int, int > m_zIdToPlane
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:66
static const NekDouble kNekZeroTol
Fourier Modified expansions with just the real part of the first mode .
Definition: BasisType.h:59
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:213
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
Definition: ExpList.h:1036
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1050
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
virtual void v_BwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:969
Array< OneD, ExpListSharedPtr > m_planes
double NekDouble
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual void v_FwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
Fourier Modified expansions with just the imaginary part of the first mode .
Definition: BasisType.h:60
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble > > &inarray1, const Array< OneD, Array< OneD, NekDouble > > &inarray2, Array< OneD, Array< OneD, NekDouble > > &outarray, CoeffState coeffstate=eLocal)
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
boost::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:72
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:966
void Homogeneous1DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
DNekBlkMatSharedPtr GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate=eLocal) const
virtual void v_IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Fourier ModifiedExpansion with just the first mode .
Definition: BasisType.h:58
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:250
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
std::map< Homogeneous1DMatType, DNekBlkMatSharedPtr > Homo1DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Describes the specification for a Basis.
Definition: Basis.h:50
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:183
LibUtilities::NektarFFTSharedPtr m_FFT