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