Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpListHomogeneous2D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpListHomogeneous2D.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 2-directions
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 #include <StdRegions/StdSegExp.h>
39 #include <StdRegions/StdQuadExp.h>
40 #include <LocalRegions/Expansion.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46  namespace MultiRegions
47  {
48  // Forward declaration for typedefs
49  ExpListHomogeneous2D::ExpListHomogeneous2D():
50  ExpList(),
51  m_homogeneousBasis_y(LibUtilities::NullBasisSharedPtr),
52  m_homogeneousBasis_z(LibUtilities::NullBasisSharedPtr),
53  m_lhom_y(1),
54  m_lhom_z(1),
55  m_homogeneous2DBlockMat(MemoryManager<Homo2DBlockMatrixMap>::AllocateSharedPtr())
56  {
57  }
58 
60  const LibUtilities::BasisKey &HomoBasis_y,
61  const LibUtilities::BasisKey &HomoBasis_z,
62  const NekDouble lhom_y,
63  const NekDouble lhom_z,
64  const bool useFFT,
65  const bool dealiasing):
66  ExpList(pSession),
67  m_useFFT(useFFT),
68  m_lhom_y(lhom_y),
69  m_lhom_z(lhom_z),
70  m_homogeneous2DBlockMat(MemoryManager<Homo2DBlockMatrixMap>::AllocateSharedPtr()),
71  m_dealiasing(dealiasing)
72  {
73  ASSERTL2(HomoBasis_y != LibUtilities::NullBasisKey,
74  "Homogeneous Basis in y direction is a null basis");
75  ASSERTL2(HomoBasis_z != LibUtilities::NullBasisKey,
76  "Homogeneous Basis in z direction is a null basis");
77 
80 
82 
83  m_Ycomm = m_comm->GetColumnComm()->GetRowComm();
84  m_Zcomm = m_comm->GetColumnComm()->GetRowComm();
85 
86  m_ny = m_homogeneousBasis_y->GetNumPoints()/m_Ycomm->GetSize();
87  m_nz = m_homogeneousBasis_z->GetNumPoints()/m_Zcomm->GetSize();
88 
90 
91  if(m_useFFT)
92  {
95  }
96 
97  if(m_dealiasing)
98  {
99  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Remove dealiasing if you want to run in parallel");
100  SetPaddingBase();
101  }
102  }
103 
104 
105  /**
106  * @param In ExpListHomogeneous2D object to copy.
107  */
109  ExpList(In,false),
110  m_useFFT(In.m_useFFT),
111  m_FFT_y(In.m_FFT_y),
112  m_FFT_z(In.m_FFT_z),
113  m_transposition(In.m_transposition),
114  m_Ycomm(In.m_Ycomm),
115  m_Zcomm(In.m_Ycomm),
116  m_homogeneousBasis_y(In.m_homogeneousBasis_y),
117  m_homogeneousBasis_z(In.m_homogeneousBasis_z),
118  m_lhom_y(In.m_lhom_y),
119  m_lhom_z(In.m_lhom_z),
120  m_homogeneous2DBlockMat(In.m_homogeneous2DBlockMat),
121  m_ny(In.m_ny),
122  m_nz(In.m_nz),
123  m_dealiasing(In.m_dealiasing),
124  m_padsize_y(In.m_padsize_y),
125  m_padsize_z(In.m_padsize_z),
126  MatFwdPAD(In.MatFwdPAD),
127  MatBwdPAD(In.MatBwdPAD)
128  {
129  m_lines = Array<OneD, ExpListSharedPtr>(In.m_lines.num_elements());
130  }
131 
132  /**
133  * Destructor
134  */
136  {
137  }
138 
140  Array<OneD, NekDouble> &outarray,
141  CoeffState coeffstate,
142  bool Shuff,
143  bool UnShuff)
144  {
145  // Forwards trans
146  Homogeneous2DTrans(inarray,outarray,true,coeffstate,Shuff,UnShuff);
147  }
148 
150  Array<OneD, NekDouble> &outarray,
151  CoeffState coeffstate,
152  bool Shuff,
153  bool UnShuff)
154  {
155  // Backwards trans
156  Homogeneous2DTrans(inarray,outarray,false,coeffstate,Shuff,UnShuff);
157  }
158 
160  const Array<OneD, NekDouble> &inarray2,
161  Array<OneD, NekDouble> &outarray,
162  CoeffState coeffstate)
163  {
164  int npoints = outarray.num_elements(); // number of total physical points
165  int nlines = m_lines.num_elements(); // number of lines == number of Fourier modes = number of Fourier coeff = number of points per slab
166  int nslabs = npoints/nlines; // number of slabs = numebr of physical points per line
167 
168  Array<OneD, NekDouble> V1(npoints);
169  Array<OneD, NekDouble> V2(npoints);
170  Array<OneD, NekDouble> V1V2(npoints);
171  Array<OneD, NekDouble> ShufV1(npoints);
172  Array<OneD, NekDouble> ShufV2(npoints);
173  Array<OneD, NekDouble> ShufV1V2(npoints);
174 
175  if(m_WaveSpace)
176  {
177  V1 = inarray1;
178  V2 = inarray2;
179  }
180  else
181  {
182  HomogeneousFwdTrans(inarray1,V1,coeffstate);
183  HomogeneousFwdTrans(inarray2,V2,coeffstate);
184  }
185 
186  m_transposition->Transpose(V1,ShufV1,false,LibUtilities::eXtoYZ);
187  m_transposition->Transpose(V2,ShufV2,false,LibUtilities::eXtoYZ);
188 
189  Array<OneD, NekDouble> PadV1_slab_coeff(m_padsize_y*m_padsize_z,0.0);
190  Array<OneD, NekDouble> PadV2_slab_coeff(m_padsize_y*m_padsize_z,0.0);
191  Array<OneD, NekDouble> PadRe_slab_coeff(m_padsize_y*m_padsize_z,0.0);
192 
193  Array<OneD, NekDouble> PadV1_slab_phys(m_padsize_y*m_padsize_z,0.0);
194  Array<OneD, NekDouble> PadV2_slab_phys(m_padsize_y*m_padsize_z,0.0);
195  Array<OneD, NekDouble> PadRe_slab_phys(m_padsize_y*m_padsize_z,0.0);
196 
197  NekVector<NekDouble> PadIN_V1(m_padsize_y*m_padsize_z,PadV1_slab_coeff,eWrapper);
198  NekVector<NekDouble> PadOUT_V1(m_padsize_y*m_padsize_z,PadV1_slab_phys,eWrapper);
199 
200  NekVector<NekDouble> PadIN_V2(m_padsize_y*m_padsize_z,PadV2_slab_coeff,eWrapper);
201  NekVector<NekDouble> PadOUT_V2(m_padsize_y*m_padsize_z,PadV2_slab_phys,eWrapper);
202 
203  NekVector<NekDouble> PadIN_Re(m_padsize_y*m_padsize_z,PadRe_slab_phys,eWrapper);
204  NekVector<NekDouble> PadOUT_Re(m_padsize_y*m_padsize_z,PadRe_slab_coeff,eWrapper);
205 
206  //Looping on the slabs
207  for(int j = 0 ; j< nslabs ; j++)
208  {
209  //Copying the j-th slab of size N*M into a bigger slab of lenght 2*N*M
210  //We are in Fourier space
211  for(int i = 0 ; i< m_nz ; i++)
212  {
213  Vmath::Vcopy(m_ny,&(ShufV1[i*m_ny + j*nlines]),1,&(PadV1_slab_coeff[i*2*m_ny]),1);
214  Vmath::Vcopy(m_ny,&(ShufV2[i*m_ny + j*nlines]),1,&(PadV2_slab_coeff[i*2*m_ny]),1);
215  }
216 
217  //Moving to physical space using the padded system
218  PadOUT_V1 = (*MatBwdPAD)*PadIN_V1;
219  PadOUT_V2 = (*MatBwdPAD)*PadIN_V2;
220 
221  //Perfroming the vectors multiplication in physical
222  //space on the padded system
223  Vmath::Vmul(m_padsize_y*m_padsize_z,PadV1_slab_phys,1,PadV2_slab_phys,1,PadRe_slab_phys,1);
224 
225  //Moving back the result (V1*V2)_phys in Fourier
226  //space, padded system
227  PadOUT_Re = (*MatFwdPAD)*PadIN_Re;
228 
229  //Copying the first half of the padded pencil in the
230  //full vector (Fourier space)
231  for (int i = 0; i < m_nz; i++)
232  {
233  Vmath::Vcopy(m_ny,&(PadRe_slab_coeff[i*2*m_ny]),1,&(ShufV1V2[i*m_ny + j*nlines]),1);
234  }
235  }
236 
237  if(m_WaveSpace)
238  {
239  m_transposition->Transpose(ShufV1V2,outarray,false,LibUtilities::eYZtoX);
240  }
241  else
242  {
243  m_transposition->Transpose(ShufV1V2,V1V2,false,LibUtilities::eYZtoX);
244 
245  //Moving the results in physical space for the output
246  HomogeneousBwdTrans(V1V2,outarray,coeffstate);
247  }
248  }
249 
251  {
252  int cnt = 0, cnt1 = 0;
253  Array<OneD, NekDouble> tmparray;
254  int nlines = m_lines.num_elements();
255 
256  for(int n = 0; n < nlines; ++n)
257  {
258  m_lines[n]->FwdTrans(inarray+cnt, tmparray = outarray + cnt1,
259  coeffstate);
260  cnt += m_lines[n]->GetTotPoints();
261  cnt1 += m_lines[n]->GetNcoeffs();
262  }
263  if(!m_WaveSpace)
264  {
265  HomogeneousFwdTrans(outarray,outarray,coeffstate);
266  }
267  }
268 
270  {
271  int cnt = 0, cnt1 = 0;
272  Array<OneD, NekDouble> tmparray;
273  int nlines = m_lines.num_elements();
274 
275  for(int n = 0; n < nlines; ++n)
276  {
277  m_lines[n]->FwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
278 
279  cnt += m_lines[n]->GetTotPoints();
280  cnt1 += m_lines[n]->GetNcoeffs();
281  }
282  if(!m_WaveSpace)
283  {
284  HomogeneousFwdTrans(outarray,outarray);
285  }
286  }
287 
289  {
290  int cnt = 0, cnt1 = 0;
291  Array<OneD, NekDouble> tmparray;
292  int nlines = m_lines.num_elements();
293 
294  for(int n = 0; n < nlines; ++n)
295  {
296  m_lines[n]->BwdTrans(inarray+cnt, tmparray = outarray + cnt1,
297  coeffstate);
298  cnt += m_lines[n]->GetNcoeffs();
299  cnt1 += m_lines[n]->GetTotPoints();
300  }
301  if(!m_WaveSpace)
302  {
303  HomogeneousBwdTrans(outarray,outarray);
304  }
305  }
306 
308  {
309  int cnt = 0, cnt1 = 0;
310  Array<OneD, NekDouble> tmparray;
311  int nlines = m_lines.num_elements();
312 
313  for(int n = 0; n < nlines; ++n)
314  {
315  m_lines[n]->BwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
316 
317  cnt += m_lines[n]->GetNcoeffs();
318  cnt1 += m_lines[n]->GetTotPoints();
319  }
320  if(!m_WaveSpace)
321  {
322  HomogeneousBwdTrans(outarray,outarray);
323  }
324  }
325 
326 
328  {
329  int cnt = 0, cnt1 = 0;
330  Array<OneD, NekDouble> tmparray;
331  int nlines = m_lines.num_elements();
332 
333  for(int n = 0; n < nlines; ++n)
334  {
335  m_lines[n]->IProductWRTBase(inarray+cnt, tmparray = outarray + cnt1,coeffstate);
336 
337  cnt += m_lines[n]->GetNcoeffs();
338  cnt1 += m_lines[n]->GetTotPoints();
339  }
340  }
341 
343  {
344  int cnt = 0, cnt1 = 0;
345  Array<OneD, NekDouble> tmparray;
346  int nlines = m_lines.num_elements();
347 
348  for(int n = 0; n < nlines; ++n)
349  {
350  m_lines[n]->IProductWRTBase_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
351 
352  cnt += m_lines[n]->GetNcoeffs();
353  cnt1 += m_lines[n]->GetTotPoints();
354  }
355  }
356 
358  Array<OneD, NekDouble> &outarray,
359  bool IsForwards,
360  CoeffState coeffstate,
361  bool Shuff,
362  bool UnShuff)
363  {
364  if(m_useFFT)
365  {
366 
367  int n = m_lines.num_elements(); //number of Fourier points in the Fourier directions (x-z grid)
368  int s = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per line
369  int p = s/n; //number of points per line = n of Fourier transform required
370 
371  Array<OneD, NekDouble> fft_in(s);
372  Array<OneD, NekDouble> fft_out(s);
373 
374  m_transposition->Transpose(inarray,fft_in,false,LibUtilities::eXtoYZ);
375 
376  if(IsForwards)
377  {
378  for(int i=0;i<(p*m_nz);i++)
379  {
380  m_FFT_y->FFTFwdTrans(m_tmpIN = fft_in + i*m_ny, m_tmpOUT = fft_out + i*m_ny);
381  }
382 
383  }
384  else
385  {
386  for(int i=0;i<(p*m_nz);i++)
387  {
388  m_FFT_y->FFTBwdTrans(m_tmpIN = fft_in + i*m_ny, m_tmpOUT = fft_out + i*m_ny);
389  }
390  }
391 
392  m_transposition->Transpose(fft_out,fft_in,false,LibUtilities::eYZtoZY);
393 
394  if(IsForwards)
395  {
396  for(int i=0;i<(p*m_ny);i++)
397  {
398  m_FFT_z->FFTFwdTrans(m_tmpIN = fft_in + i*m_nz, m_tmpOUT = fft_out + i*m_nz);
399  }
400 
401  }
402  else
403  {
404  for(int i=0;i<(p*m_ny);i++)
405  {
406  m_FFT_z->FFTBwdTrans(m_tmpIN = fft_in + i*m_nz, m_tmpOUT = fft_out + i*m_nz);
407  }
408  }
409 
410  //TODO: required ZYtoX routine
411  m_transposition->Transpose(fft_out,fft_in,false,LibUtilities::eZYtoYZ);
412 
413  m_transposition->Transpose(fft_in,outarray,false,LibUtilities::eYZtoX);
414 
415  }
416  else
417  {
418  DNekBlkMatSharedPtr blkmatY;
419  DNekBlkMatSharedPtr blkmatZ;
420 
421  if(inarray.num_elements() == m_npoints) //transform phys space
422  {
423  if(IsForwards)
424  {
427  }
428  else
429  {
432  }
433  }
434  else
435  {
436  if(IsForwards)
437  {
440  }
441  else
442  {
445  }
446  }
447 
448  int nrowsY = blkmatY->GetRows();
449  int ncolsY = blkmatY->GetColumns();
450 
451  Array<OneD, NekDouble> sortedinarrayY(ncolsY);
452  Array<OneD, NekDouble> sortedoutarrayY(nrowsY);
453 
454  int nrowsZ = blkmatZ->GetRows();
455  int ncolsZ = blkmatZ->GetColumns();
456 
457  Array<OneD, NekDouble> sortedinarrayZ(ncolsZ);
458  Array<OneD, NekDouble> sortedoutarrayZ(nrowsZ);
459 
460  NekVector<NekDouble> inY (ncolsY,sortedinarrayY,eWrapper);
461  NekVector<NekDouble> outY(nrowsY,sortedoutarrayY,eWrapper);
462 
463  NekVector<NekDouble> inZ (ncolsZ,sortedinarrayZ,eWrapper);
464  NekVector<NekDouble> outZ(nrowsZ,sortedoutarrayZ,eWrapper);
465 
466  m_transposition->Transpose(inarray,sortedinarrayY,!IsForwards,LibUtilities::eXtoYZ);
467 
468  outY = (*blkmatY)*inY;
469 
470  m_transposition->Transpose(sortedoutarrayY,sortedinarrayZ,false,LibUtilities::eYZtoZY);
471 
472  outZ = (*blkmatZ)*inZ;
473 
474  m_transposition->Transpose(sortedoutarrayZ,sortedoutarrayY,false,LibUtilities::eZYtoYZ);
475 
476  m_transposition->Transpose(sortedoutarrayY,outarray,false,LibUtilities::eYZtoX);
477  }
478  }
479 
481  {
482  Homo2DBlockMatrixMap::iterator matrixIter = m_homogeneous2DBlockMat->find(mattype);
483 
484  if(matrixIter == m_homogeneous2DBlockMat->end())
485  {
486  return ((*m_homogeneous2DBlockMat)[mattype] =
487  GenHomogeneous2DBlockMatrix(mattype,coeffstate));
488  }
489  else
490  {
491  return matrixIter->second;
492  }
493  }
494 
496  {
497  int i;
498  int n_exp = 0;
499 
500  DNekMatSharedPtr loc_mat;
501  DNekBlkMatSharedPtr BlkMatrix;
502 
504 
505  int NumPoints = 0;
506  int NumModes = 0;
507  int NumPencils = 0;
508 
509  if((mattype == eForwardsCoeffSpaceY1D) || (mattype == eBackwardsCoeffSpaceY1D)
510  ||(mattype == eForwardsPhysSpaceY1D) || (mattype == eBackwardsPhysSpaceY1D))
511  {
512  Basis = m_homogeneousBasis_y;
513  NumPoints = m_homogeneousBasis_y->GetNumModes();
514  NumModes = m_homogeneousBasis_y->GetNumPoints();
515  NumPencils = m_homogeneousBasis_z->GetNumPoints();
516  }
517  else
518  {
519  Basis = m_homogeneousBasis_z;
520  NumPoints = m_homogeneousBasis_z->GetNumModes();
521  NumModes = m_homogeneousBasis_z->GetNumPoints();
522  NumPencils = m_homogeneousBasis_y->GetNumPoints();
523  }
524 
525  if((mattype == eForwardsCoeffSpaceY1D) || (mattype == eForwardsCoeffSpaceZ1D)
526  ||(mattype == eBackwardsCoeffSpaceY1D)||(mattype == eBackwardsCoeffSpaceZ1D))
527  {
528  n_exp = m_lines[0]->GetNcoeffs();
529  }
530  else
531  {
532  n_exp = m_lines[0]->GetTotPoints(); // will operatore on m_phys
533  }
534 
535  Array<OneD,unsigned int> nrows(n_exp);
536  Array<OneD,unsigned int> ncols(n_exp);
537 
538  if((mattype == eForwardsCoeffSpaceY1D)||(mattype == eForwardsPhysSpaceY1D) ||
539  (mattype == eForwardsCoeffSpaceZ1D)||(mattype == eForwardsPhysSpaceZ1D))
540  {
541  nrows = Array<OneD, unsigned int>(n_exp*NumPencils,NumModes);
542  ncols = Array<OneD, unsigned int>(n_exp*NumPencils,NumPoints);
543  }
544  else
545  {
546  nrows = Array<OneD, unsigned int>(n_exp*NumPencils,NumPoints);
547  ncols = Array<OneD, unsigned int>(n_exp*NumPencils,NumModes);
548  }
549 
550  MatrixStorage blkmatStorage = eDIAGONAL;
551  BlkMatrix = MemoryManager<DNekBlkMat>::AllocateSharedPtr(nrows,ncols,blkmatStorage);
552 
553  StdRegions::StdSegExp StdSeg(Basis->GetBasisKey());
554 
555  if((mattype == eForwardsCoeffSpaceY1D)||(mattype == eForwardsPhysSpaceY1D) ||
556  (mattype == eForwardsCoeffSpaceZ1D)||(mattype == eForwardsPhysSpaceZ1D))
557  {
559  StdSeg.DetShapeType(),
560  StdSeg);
561 
562  loc_mat = StdSeg.GetStdMatrix(matkey);
563  }
564  else
565  {
567  StdSeg.DetShapeType(),
568  StdSeg);
569 
570  loc_mat = StdSeg.GetStdMatrix(matkey);
571  }
572 
573  // set up array of block matrices.
574  for(i = 0; i < (n_exp*NumPencils); ++i)
575  {
576  BlkMatrix->SetBlock(i,i,loc_mat);
577  }
578 
579  return BlkMatrix;
580  }
581 
582  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous2D::v_GetFieldDefinitions()
583  {
584  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
585  // Set up Homogeneous length details.
587  HomoBasis[0] = m_homogeneousBasis_y;
588  HomoBasis[1] = m_homogeneousBasis_z;
589 
590  std::vector<NekDouble> HomoLen(2);
591  HomoLen[0] = m_lhom_y;
592  HomoLen[1] = m_lhom_z;
593 
594  int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
595  int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
596 
597  std::vector<unsigned int> sIDs
599 
600  std::vector<unsigned int> yIDs;
601  std::vector<unsigned int> zIDs;
602 
603  for(int n = 0; n < nhom_modes_z; ++n)
604  {
605  for(int m = 0; m < nhom_modes_y; ++m)
606  {
607  zIDs.push_back(n);
608  yIDs.push_back(m);
609  }
610  }
611 
612  m_lines[0]->GeneralGetFieldDefinitions(returnval, 2, HomoBasis,
613  HomoLen, false,
614  sIDs, zIDs, yIDs);
615  return returnval;
616  }
617 
618  void ExpListHomogeneous2D::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
619  {
620  // Set up Homogeneous length details.
622  HomoBasis[0] = m_homogeneousBasis_y;
623  HomoBasis[1] = m_homogeneousBasis_z;
624  std::vector<NekDouble> HomoLen(2);
625  HomoLen[0] = m_lhom_y;
626  HomoLen[1] = m_lhom_z;
627 
628  int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
629  int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
630 
631  std::vector<unsigned int> sIDs
633 
634  std::vector<unsigned int> yIDs;
635  std::vector<unsigned int> zIDs;
636 
637  for(int n = 0; n < nhom_modes_z; ++n)
638  {
639  for(int m = 0; m < nhom_modes_y; ++m)
640  {
641  zIDs.push_back(n);
642  yIDs.push_back(m);
643  }
644  }
645 
646  // enforce NumHomoDir == 1 by direct call
647  m_lines[0]->GeneralGetFieldDefinitions(fielddef, 2, HomoBasis,
648  HomoLen, false,
649  sIDs, zIDs, yIDs);
650  }
651 
653  {
654  int i,k;
655 
656  int NumMod_y = m_homogeneousBasis_y->GetNumModes();
657  int NumMod_z = m_homogeneousBasis_z->GetNumModes();
658 
659  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
660 
661  // Determine mapping from element ids to location in
662  // expansion list
663  map<int, int> ElmtID_to_ExpID;
664  for(i = 0; i < m_lines[0]->GetExpSize(); ++i)
665  {
666  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
667  }
668 
669  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
670  {
671  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
672  int datalen = (*m_exp)[eid]->GetNcoeffs();
673 
674  for(k = 0; k < (NumMod_y*NumMod_z); ++k)
675  {
676  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]+k*ncoeffs_per_line],&coeffs[m_coeff_offset[eid]+k*ncoeffs_per_line]+datalen);
677  }
678  }
679  }
680 
681  void ExpListHomogeneous2D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
682  {
683  v_AppendFieldData(fielddef,fielddata,m_coeffs);
684  }
685 
686  //Extract the data in fielddata into the m_coeff list
687  void ExpListHomogeneous2D::v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, std::string &field, Array<OneD, NekDouble> &coeffs)
688  {
689  int i,k;
690  int offset = 0;
691  int datalen = fielddata.size()/fielddef->m_fields.size();
692  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
693  int NumMod_y = m_homogeneousBasis_y->GetNumModes();
694  int NumMod_z = m_homogeneousBasis_z->GetNumModes();
695 
696  // Find data location according to field definition
697  for(i = 0; i < fielddef->m_fields.size(); ++i)
698  {
699  if(fielddef->m_fields[i] == field)
700  {
701  break;
702  }
703  offset += datalen;
704  }
705 
706  ASSERTL0(i!= fielddef->m_fields.size(),"Field not found in data file");
707 
708  // Determine mapping from element ids to location in
709  // expansion list
710  map<int, int> ElmtID_to_ExpID;
711  for(i = 0; i < m_lines[0]->GetExpSize(); ++i)
712  {
713  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
714  }
715 
716  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
717  {
718  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
719  int datalen = (*m_exp)[eid]->GetNcoeffs();
720 
721  for(k = 0; k < (NumMod_y*NumMod_z); ++k)
722  {
723  Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid] + k*ncoeffs_per_line],1);
724  offset += datalen;
725  }
726  }
727  }
728 
729  void ExpListHomogeneous2D::v_WriteVtkPieceData(std::ostream &outfile, int expansion,
730  std::string var)
731  {
732  int i;
733  int nq = (*m_exp)[expansion]->GetTotPoints();
734  int npoints_per_line = m_lines[0]->GetTotPoints();
735 
736  // printing the fields of that zone
737  outfile << " <DataArray type=\"Float32\" Name=\""
738  << var << "\">" << endl;
739  outfile << " ";
740  for (int n = 0; n < m_lines.num_elements(); ++n)
741  {
742  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_line;
743  for(i = 0; i < nq; ++i)
744  {
745  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
746  }
747  }
748  outfile << endl;
749  outfile << " </DataArray>" << endl;
750  }
751 
753  Array<OneD, NekDouble> &out_d0,
754  Array<OneD, NekDouble> &out_d1,
755  Array<OneD, NekDouble> &out_d2)
756 
757  {
758  int nyzlines = m_lines.num_elements(); //number of Fourier points in the Fourier directions (nF_pts)
759  int npoints = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per line (nT_pts)
760  int n_points_line = npoints/nyzlines; //number of points per line
761 
762  Array<OneD, NekDouble> temparray(npoints);
763  Array<OneD, NekDouble> temparray1(npoints);
764  Array<OneD, NekDouble> temparray2(npoints);
768 
769  for( int i=0 ; i<nyzlines ; i++ )
770  {
771  m_lines[i]->PhysDeriv( tmp1 = inarray + i*n_points_line ,tmp2 = out_d0 + i*n_points_line);
772  }
773 
775  {
776  if(m_WaveSpace)
777  {
778  temparray = inarray;
779  }
780  else
781  {
782  HomogeneousFwdTrans(inarray,temparray);
783  }
784  NekDouble sign = -1.0;
785  NekDouble beta;
786 
787  //along y
788  for(int i = 0; i < m_ny; i++)
789  {
790  beta = -sign*2*M_PI*(i/2)/m_lhom_y;
791 
792  for(int j = 0; j < m_nz; j++)
793  {
794  Vmath::Smul(n_points_line,beta,tmp1 = temparray + n_points_line*(i+j*m_ny),1, tmp2 = temparray1 + n_points_line*((i-int(sign))+j*m_ny),1);
795  }
796 
797  sign = -1.0*sign;
798  }
799 
800  //along z
801  sign = -1.0;
802  for(int i = 0; i < m_nz; i++)
803  {
804  beta = -sign*2*M_PI*(i/2)/m_lhom_z;
805  Vmath::Smul(m_ny*n_points_line,beta,tmp1 = temparray + i*m_ny*n_points_line,1,tmp2 = temparray2 + (i-int(sign))*m_ny*n_points_line,1);
806  sign = -1.0*sign;
807  }
808  if(m_WaveSpace)
809  {
810  out_d1 = temparray1;
811  out_d2 = temparray2;
812  }
813  else
814  {
815  HomogeneousBwdTrans(temparray1,out_d1);
816  HomogeneousBwdTrans(temparray2,out_d2);
817  }
818  }
819  else
820  {
821  if(m_WaveSpace)
822  {
823  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis")
824  }
825  else
826  {
827  StdRegions::StdQuadExp StdQuad(m_homogeneousBasis_y->GetBasisKey(),m_homogeneousBasis_z->GetBasisKey());
828 
829  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXtoYZ);
830 
831  for(int i = 0; i < n_points_line; i++)
832  {
833  StdQuad.PhysDeriv(tmp1 = temparray + i*nyzlines, tmp2 = temparray1 + i*nyzlines, tmp3 = temparray2 + i*nyzlines);
834  }
835 
836  m_transposition->Transpose(temparray1,out_d1,false,LibUtilities::eYZtoX);
837  m_transposition->Transpose(temparray2,out_d2,false,LibUtilities::eYZtoX);
838  Vmath::Smul(npoints,2.0/m_lhom_y,out_d1,1,out_d1,1);
839  Vmath::Smul(npoints,2.0/m_lhom_z,out_d2,1,out_d2,1);
840  }
841  }
842  }
843 
845  const Array<OneD, const NekDouble> &inarray,
846  Array<OneD, NekDouble> &out_d)
847 
848  {
849  int nyzlines = m_lines.num_elements(); //number of Fourier points in the Fourier directions (nF_pts)
850  int npoints = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per line (nT_pts)
851  int n_points_line = npoints/nyzlines; //number of points per line
852  //convert enum into int
853  int dir = (int)edir;
854 
855  Array<OneD, NekDouble> temparray(npoints);
856  Array<OneD, NekDouble> temparray1(npoints);
857  Array<OneD, NekDouble> temparray2(npoints);
861 
862  if (dir < 1)
863  {
864  for( int i=0 ; i<nyzlines ; i++)
865  {
866  m_lines[i]->PhysDeriv( tmp1 = inarray + i*n_points_line ,tmp2 = out_d + i*n_points_line);
867  }
868  }
869  else
870  {
872  {
873  if(m_WaveSpace)
874  {
875  temparray = inarray;
876  }
877  else
878  {
879  HomogeneousFwdTrans(inarray,temparray);
880  }
881  NekDouble sign = -1.0;
882  NekDouble beta;
883 
884  if (dir == 1)
885  {
886  //along y
887  for(int i = 0; i < m_ny; i++)
888  {
889  beta = -sign*2*M_PI*(i/2)/m_lhom_y;
890 
891  for(int j = 0; j < m_nz; j++)
892  {
893  Vmath::Smul(n_points_line,beta,tmp1 = temparray + n_points_line*(i+j*m_ny),1, tmp2 = temparray1 + n_points_line*((i-int(sign))+j*m_ny),1);
894  }
895  sign = -1.0*sign;
896  }
897  if(m_WaveSpace)
898  {
899  out_d = temparray1;
900  }
901  else
902  {
903  HomogeneousBwdTrans(temparray1,out_d);
904  }
905  }
906  else
907  {
908  //along z
909  for(int i = 0; i < m_nz; i++)
910  {
911  beta = -sign*2*M_PI*(i/2)/m_lhom_z;
912  Vmath::Smul(m_ny*n_points_line,beta,tmp1 = temparray + i*m_ny*n_points_line,1,tmp2 = temparray2 + (i-int(sign))*m_ny*n_points_line,1);
913  sign = -1.0*sign;
914  }
915  if(m_WaveSpace)
916  {
917  out_d = temparray2;
918  }
919  else
920  {
921  HomogeneousBwdTrans(temparray2,out_d);
922  }
923  }
924  }
925  else
926  {
927  if(m_WaveSpace)
928  {
929  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis")
930  }
931  else
932  {
933  StdRegions::StdQuadExp StdQuad(m_homogeneousBasis_y->GetBasisKey(),m_homogeneousBasis_z->GetBasisKey());
934 
935  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXtoYZ);
936 
937  for(int i = 0; i < n_points_line; i++)
938  {
939  StdQuad.PhysDeriv(tmp1 = temparray + i*nyzlines, tmp2 = temparray1 + i*nyzlines, tmp3 = temparray2 + i*nyzlines);
940  }
941 
942  if (dir == 1)
943  {
944  m_transposition->Transpose(temparray1,out_d,false,LibUtilities::eYZtoX);
945  Vmath::Smul(npoints,2.0/m_lhom_y,out_d,1,out_d,1);
946  }
947  else
948  {
949  m_transposition->Transpose(temparray2,out_d,false,LibUtilities::eYZtoX);
950  Vmath::Smul(npoints,2.0/m_lhom_z,out_d,1,out_d,1);
951  }
952  }
953  }
954  }
955  }
956 
958  Array<OneD, NekDouble> &out_d0,
959  Array<OneD, NekDouble> &out_d1,
960  Array<OneD, NekDouble> &out_d2)
961 
962  {
963  v_PhysDeriv(inarray,out_d0,out_d1,out_d2);
964  }
965 
967  const Array<OneD, const NekDouble> &inarray,
968  Array<OneD, NekDouble> &out_d)
969  {
970  //convert int into enum
971  v_PhysDeriv(edir,inarray,out_d);
972  }
973 
975  {
976  NekDouble size_y = 1.5*m_ny;
977  NekDouble size_z = 1.5*m_nz;
978  m_padsize_y = int(size_y);
979  m_padsize_z = int(size_z);
980 
983 
986 
989 
990  StdRegions::StdQuadExp StdQuad(m_paddingBasis_y->GetBasisKey(),m_paddingBasis_z->GetBasisKey());
991 
992  StdRegions::StdMatrixKey matkey1(StdRegions::eFwdTrans,StdQuad.DetShapeType(),StdQuad);
993  StdRegions::StdMatrixKey matkey2(StdRegions::eBwdTrans,StdQuad.DetShapeType(),StdQuad);
994 
995  MatFwdPAD = StdQuad.GetStdMatrix(matkey1);
996  MatBwdPAD = StdQuad.GetStdMatrix(matkey2);
997  }
998  } //end of namespace
999 } //end of namespace
1000 
1001 
1002 /**
1003 * $Log: v $
1004 *
1005 **/
1006 
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
LibUtilities::TranspositionSharedPtr m_transposition
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
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
DNekBlkMatSharedPtr GenHomogeneous2DBlockMatrix(Homogeneous2DMatType mattype, CoeffState coeffstate=eLocal) const
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:22
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
LibUtilities::BasisSharedPtr m_paddingBasis_z
Base expansion in z direction.
boost::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:118
static BasisSharedPtr NullBasisSharedPtr
Definition: Basis.h:358
STL namespace.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:956
Fourier Expansion .
Definition: BasisType.h:52
LibUtilities::NektarFFTSharedPtr m_FFT_z
NekDouble m_lhom_z
Width of homogeneous direction z.
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal)
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:939
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
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 void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
Array< OneD, ExpListSharedPtr > m_lines
Vector of ExpList, will be filled with ExpList1D.
BasisManagerT & BasisManager(void)
virtual void v_IProductWRTBase_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:988
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:64
static const NekDouble kNekZeroTol
LibUtilities::BasisSharedPtr m_homogeneousBasis_z
Base expansion in z direction.
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:199
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:991
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
NekDouble m_lhom_y
Width of homogeneous direction y.
static std::vector< unsigned int > NullUnsignedIntVector
Definition: FieldIO.h:51
Defines a specification for a set of points.
Definition: Points.h:58
double NekDouble
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
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)
int m_ny
Number of modes = number of poitns in y direction.
LibUtilities::BasisSharedPtr m_homogeneousBasis_y
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
void Homogeneous2DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
std::map< Homogeneous2DMatType, DNekBlkMatSharedPtr > Homo2DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
int m_nz
Number of modes = number of poitns in z direction.
LibUtilities::BasisSharedPtr m_paddingBasis_y
Base expansion in y direction.
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_BwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:907
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed t...
Definition: ErrorUtil.hpp:213
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
LibUtilities::NektarFFTSharedPtr m_FFT_y
virtual void v_FwdTrans_IterPerExp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
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< Basis > BasisSharedPtr
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1047
Describes the specification for a Basis.
Definition: Basis.h:50
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate)
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:169
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, CoeffState coeffstate=eLocal, bool Shuff=true, bool UnShuff=true)
DNekBlkMatSharedPtr GetHomogeneous2DBlockMatrix(Homogeneous2DMatType mattype, CoeffState coeffstate=eLocal) const