Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: An ExpList which is homogeneous in 2-directions
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 #include <LocalRegions/Expansion.h>
40 #include <StdRegions/StdQuadExp.h>
41 #include <StdRegions/StdSegExp.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace MultiRegions
48 {
49 // Forward declaration for typedefs
50 ExpListHomogeneous2D::ExpListHomogeneous2D(const ExpansionType type)
51  : ExpList(type), m_homogeneousBasis_y(LibUtilities::NullBasisSharedPtr),
52  m_homogeneousBasis_z(LibUtilities::NullBasisSharedPtr), m_lhom_y(1),
53  m_lhom_z(1), m_homogeneous2DBlockMat(
54  MemoryManager<Homo2DBlockMatrixMap>::AllocateSharedPtr())
55 {
56 }
57 
59  const ExpansionType type,
61  const LibUtilities::BasisKey &HomoBasis_y,
62  const LibUtilities::BasisKey &HomoBasis_z, const NekDouble lhom_y,
63  const NekDouble lhom_z, const bool useFFT, const bool dealiasing)
64  : ExpList(type), m_useFFT(useFFT), m_lhom_y(lhom_y), m_lhom_z(lhom_z),
65  m_homogeneous2DBlockMat(
66  MemoryManager<Homo2DBlockMatrixMap>::AllocateSharedPtr()),
67  m_dealiasing(dealiasing)
68 {
69  m_session = pSession;
70  m_comm = pSession->GetComm();
71 
72  ASSERTL2(HomoBasis_y != LibUtilities::NullBasisKey,
73  "Homogeneous Basis in y direction is a null basis");
74  ASSERTL2(HomoBasis_z != LibUtilities::NullBasisKey,
75  "Homogeneous Basis in z direction is a null basis");
76 
79 
82  HomoBasis_y, HomoBasis_z, m_comm->GetColumnComm());
83 
84  m_Ycomm = m_comm->GetColumnComm()->GetRowComm();
85  m_Zcomm = m_comm->GetColumnComm()->GetRowComm();
86 
87  m_ny = m_homogeneousBasis_y->GetNumPoints() / m_Ycomm->GetSize();
88  m_nz = m_homogeneousBasis_z->GetNumPoints() / m_Zcomm->GetSize();
89 
91 
92  if (m_useFFT)
93  {
94  m_FFT_y =
96  m_FFT_z =
98  }
99 
100  if (m_dealiasing)
101  {
102  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,
103  "Remove dealiasing if you want to run in parallel");
104  SetPaddingBase();
105  }
106 }
107 
108 /**
109  * @param In ExpListHomogeneous2D object to copy.
110  */
112  : ExpList(In, false), m_useFFT(In.m_useFFT), m_FFT_y(In.m_FFT_y),
113  m_FFT_z(In.m_FFT_z), m_transposition(In.m_transposition),
114  m_Ycomm(In.m_Ycomm), m_Zcomm(In.m_Ycomm),
115  m_homogeneousBasis_y(In.m_homogeneousBasis_y),
116  m_homogeneousBasis_z(In.m_homogeneousBasis_z), m_lhom_y(In.m_lhom_y),
117  m_lhom_z(In.m_lhom_z),
118  m_homogeneous2DBlockMat(In.m_homogeneous2DBlockMat), m_ny(In.m_ny),
119  m_nz(In.m_nz), m_dealiasing(In.m_dealiasing), m_padsize_y(In.m_padsize_y),
120  m_padsize_z(In.m_padsize_z), MatFwdPAD(In.MatFwdPAD),
121  MatBwdPAD(In.MatBwdPAD)
122 {
124 }
125 
127  const ExpListHomogeneous2D &In, const std::vector<unsigned int> &eIDs)
128  : ExpList(In, eIDs, false), m_useFFT(In.m_useFFT), m_FFT_y(In.m_FFT_y),
129  m_FFT_z(In.m_FFT_z), m_transposition(In.m_transposition),
130  m_Ycomm(In.m_Ycomm), m_Zcomm(In.m_Ycomm),
131  m_homogeneousBasis_y(In.m_homogeneousBasis_y),
132  m_homogeneousBasis_z(In.m_homogeneousBasis_z), m_lhom_y(In.m_lhom_y),
133  m_lhom_z(In.m_lhom_z),
134  m_homogeneous2DBlockMat(
135  MemoryManager<Homo2DBlockMatrixMap>::AllocateSharedPtr()),
136  m_ny(In.m_ny), m_nz(In.m_nz), m_dealiasing(In.m_dealiasing),
137  m_padsize_y(In.m_padsize_y), m_padsize_z(In.m_padsize_z),
138  MatFwdPAD(In.MatFwdPAD), MatBwdPAD(In.MatBwdPAD)
139 {
141 }
142 
143 /**
144  * Destructor
145  */
147 {
148 }
149 
151  const Array<OneD, const NekDouble> &inarray,
152  Array<OneD, NekDouble> &outarray, bool Shuff, bool UnShuff)
153 {
154  // Forwards trans
155  Homogeneous2DTrans(inarray, outarray, true, Shuff, UnShuff);
156 }
157 
159  const Array<OneD, const NekDouble> &inarray,
160  Array<OneD, NekDouble> &outarray, bool Shuff, bool UnShuff)
161 {
162  // Backwards trans
163  Homogeneous2DTrans(inarray, outarray, false, Shuff, UnShuff);
164 }
165 
167  const Array<OneD, NekDouble> &inarray1,
168  const Array<OneD, NekDouble> &inarray2, Array<OneD, NekDouble> &outarray)
169 {
170  int npoints = outarray.size(); // number of total physical points
171  int nlines =
172  m_lines.size(); // number of lines == number of Fourier modes = number
173  // of Fourier coeff = number of points per slab
174  int nslabs = npoints /
175  nlines; // number of slabs = numebr of physical points per line
176 
177  Array<OneD, NekDouble> V1(npoints);
178  Array<OneD, NekDouble> V2(npoints);
179  Array<OneD, NekDouble> V1V2(npoints);
180  Array<OneD, NekDouble> ShufV1(npoints);
181  Array<OneD, NekDouble> ShufV2(npoints);
182  Array<OneD, NekDouble> ShufV1V2(npoints);
183 
184  if (m_WaveSpace)
185  {
186  V1 = inarray1;
187  V2 = inarray2;
188  }
189  else
190  {
191  HomogeneousFwdTrans(inarray1, V1);
192  HomogeneousFwdTrans(inarray2, V2);
193  }
194 
195  m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXtoYZ);
196  m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXtoYZ);
197 
198  Array<OneD, NekDouble> PadV1_slab_coeff(m_padsize_y * m_padsize_z, 0.0);
199  Array<OneD, NekDouble> PadV2_slab_coeff(m_padsize_y * m_padsize_z, 0.0);
200  Array<OneD, NekDouble> PadRe_slab_coeff(m_padsize_y * m_padsize_z, 0.0);
201 
202  Array<OneD, NekDouble> PadV1_slab_phys(m_padsize_y * m_padsize_z, 0.0);
203  Array<OneD, NekDouble> PadV2_slab_phys(m_padsize_y * m_padsize_z, 0.0);
204  Array<OneD, NekDouble> PadRe_slab_phys(m_padsize_y * m_padsize_z, 0.0);
205 
206  NekVector<NekDouble> PadIN_V1(m_padsize_y * m_padsize_z, PadV1_slab_coeff,
207  eWrapper);
208  NekVector<NekDouble> PadOUT_V1(m_padsize_y * m_padsize_z, PadV1_slab_phys,
209  eWrapper);
210 
211  NekVector<NekDouble> PadIN_V2(m_padsize_y * m_padsize_z, PadV2_slab_coeff,
212  eWrapper);
213  NekVector<NekDouble> PadOUT_V2(m_padsize_y * m_padsize_z, PadV2_slab_phys,
214  eWrapper);
215 
216  NekVector<NekDouble> PadIN_Re(m_padsize_y * m_padsize_z, PadRe_slab_phys,
217  eWrapper);
218  NekVector<NekDouble> PadOUT_Re(m_padsize_y * m_padsize_z, PadRe_slab_coeff,
219  eWrapper);
220 
221  // Looping on the slabs
222  for (int j = 0; j < nslabs; j++)
223  {
224  // Copying the j-th slab of size N*M into a bigger slab of lenght 2*N*M
225  // We are in Fourier space
226  for (int i = 0; i < m_nz; i++)
227  {
228  Vmath::Vcopy(m_ny, &(ShufV1[i * m_ny + j * nlines]), 1,
229  &(PadV1_slab_coeff[i * 2 * m_ny]), 1);
230  Vmath::Vcopy(m_ny, &(ShufV2[i * m_ny + j * nlines]), 1,
231  &(PadV2_slab_coeff[i * 2 * m_ny]), 1);
232  }
233 
234  // Moving to physical space using the padded system
235  PadOUT_V1 = (*MatBwdPAD) * PadIN_V1;
236  PadOUT_V2 = (*MatBwdPAD) * PadIN_V2;
237 
238  // Perfroming the vectors multiplication in physical
239  // space on the padded system
240  Vmath::Vmul(m_padsize_y * m_padsize_z, PadV1_slab_phys, 1,
241  PadV2_slab_phys, 1, PadRe_slab_phys, 1);
242 
243  // Moving back the result (V1*V2)_phys in Fourier
244  // space, padded system
245  PadOUT_Re = (*MatFwdPAD) * PadIN_Re;
246 
247  // Copying the first half of the padded pencil in the
248  // full vector (Fourier space)
249  for (int i = 0; i < m_nz; i++)
250  {
251  Vmath::Vcopy(m_ny, &(PadRe_slab_coeff[i * 2 * m_ny]), 1,
252  &(ShufV1V2[i * m_ny + j * nlines]), 1);
253  }
254  }
255 
256  if (m_WaveSpace)
257  {
258  m_transposition->Transpose(ShufV1V2, outarray, false,
260  }
261  else
262  {
263  m_transposition->Transpose(ShufV1V2, V1V2, false, LibUtilities::eYZtoX);
264 
265  // Moving the results in physical space for the output
266  HomogeneousBwdTrans(V1V2, outarray);
267  }
268 }
269 
271  const Array<OneD, Array<OneD, NekDouble>> &inarray1,
272  const Array<OneD, Array<OneD, NekDouble>> &inarray2,
273  Array<OneD, Array<OneD, NekDouble>> &outarray)
274 {
275  // TODO Proper implementation of this
276  int ndim = inarray1.size();
277  ASSERTL1(inarray2.size() % ndim == 0,
278  "Wrong dimensions for DealiasedDotProd.");
279  int nvec = inarray2.size() % ndim;
280  int npts = inarray1[0].size();
281 
282  Array<OneD, NekDouble> out(npts);
283  for (int i = 0; i < nvec; i++)
284  {
285  Vmath::Zero(npts, outarray[i], 1);
286  for (int j = 0; j < ndim; j++)
287  {
288  DealiasedProd(inarray1[j], inarray2[i * ndim + j], out);
289  Vmath::Vadd(npts, outarray[i], 1, out, 1, outarray[i], 1);
290  }
291  }
292 }
293 
295  const Array<OneD, const NekDouble> &inarray,
296  Array<OneD, NekDouble> &outarray)
297 {
298  int cnt = 0, cnt1 = 0;
299  Array<OneD, NekDouble> tmparray;
300  int nlines = m_lines.size();
301 
302  for (int n = 0; n < nlines; ++n)
303  {
304  m_lines[n]->FwdTrans(inarray + cnt, tmparray = outarray + cnt1);
305  cnt += m_lines[n]->GetTotPoints();
306  cnt1 += m_lines[n]->GetNcoeffs();
307  }
308  if (!m_WaveSpace)
309  {
310  HomogeneousFwdTrans(outarray, outarray);
311  }
312 }
313 
315  const Array<OneD, const NekDouble> &inarray,
316  Array<OneD, NekDouble> &outarray)
317 {
318  int cnt = 0, cnt1 = 0;
319  Array<OneD, NekDouble> tmparray;
320  int nlines = m_lines.size();
321 
322  for (int n = 0; n < nlines; ++n)
323  {
324  m_lines[n]->FwdTransLocalElmt(inarray + cnt,
325  tmparray = outarray + cnt1);
326 
327  cnt += m_lines[n]->GetTotPoints();
328  cnt1 += m_lines[n]->GetNcoeffs();
329  }
330  if (!m_WaveSpace)
331  {
332  HomogeneousFwdTrans(outarray, outarray);
333  }
334 }
335 
337  const Array<OneD, const NekDouble> &inarray,
338  Array<OneD, NekDouble> &outarray)
339 {
340  int cnt = 0, cnt1 = 0;
341  Array<OneD, NekDouble> tmparray;
342  int nlines = m_lines.size();
343 
344  for (int n = 0; n < nlines; ++n)
345  {
346  m_lines[n]->BwdTrans(inarray + cnt, tmparray = outarray + cnt1);
347  cnt += m_lines[n]->GetNcoeffs();
348  cnt1 += m_lines[n]->GetTotPoints();
349  }
350  if (!m_WaveSpace)
351  {
352  HomogeneousBwdTrans(outarray, outarray);
353  }
354 }
355 
357  const Array<OneD, const NekDouble> &inarray,
358  Array<OneD, NekDouble> &outarray)
359 {
360  int cnt = 0, cnt1 = 0;
361  Array<OneD, NekDouble> tmparray;
362  int nlines = m_lines.size();
363 
364  for (int n = 0; n < nlines; ++n)
365  {
366  m_lines[n]->IProductWRTBase(inarray + cnt, tmparray = outarray + cnt1);
367 
368  cnt += m_lines[n]->GetNcoeffs();
369  cnt1 += m_lines[n]->GetTotPoints();
370  }
371 }
372 
374  const Array<OneD, const NekDouble> &inarray,
375  Array<OneD, NekDouble> &outarray, bool IsForwards, bool Shuff, bool UnShuff)
376 {
377  boost::ignore_unused(Shuff, UnShuff);
378 
379  if (m_useFFT)
380  {
381 
382  int n = m_lines.size(); // number of Fourier points in the Fourier
383  // directions (x-z grid)
384  int s = inarray.size(); // number of total points = n. of Fourier points
385  // * n. of points per line
386  int p =
387  s /
388  n; // number of points per line = n of Fourier transform required
389 
390  Array<OneD, NekDouble> fft_in(s);
391  Array<OneD, NekDouble> fft_out(s);
392 
393  m_transposition->Transpose(inarray, fft_in, false,
395 
396  if (IsForwards)
397  {
398  for (int i = 0; i < (p * m_nz); i++)
399  {
400  m_FFT_y->FFTFwdTrans(m_tmpIN = fft_in + i * m_ny,
401  m_tmpOUT = fft_out + i * m_ny);
402  }
403  }
404  else
405  {
406  for (int i = 0; i < (p * m_nz); i++)
407  {
408  m_FFT_y->FFTBwdTrans(m_tmpIN = fft_in + i * m_ny,
409  m_tmpOUT = fft_out + i * m_ny);
410  }
411  }
412 
413  m_transposition->Transpose(fft_out, fft_in, false,
415 
416  if (IsForwards)
417  {
418  for (int i = 0; i < (p * m_ny); i++)
419  {
420  m_FFT_z->FFTFwdTrans(m_tmpIN = fft_in + i * m_nz,
421  m_tmpOUT = fft_out + i * m_nz);
422  }
423  }
424  else
425  {
426  for (int i = 0; i < (p * m_ny); i++)
427  {
428  m_FFT_z->FFTBwdTrans(m_tmpIN = fft_in + i * m_nz,
429  m_tmpOUT = fft_out + i * m_nz);
430  }
431  }
432 
433  // TODO: required ZYtoX routine
434  m_transposition->Transpose(fft_out, fft_in, false,
436 
437  m_transposition->Transpose(fft_in, outarray, false,
439  }
440  else
441  {
442  DNekBlkMatSharedPtr blkmatY;
443  DNekBlkMatSharedPtr blkmatZ;
444 
445  if (inarray.size() == m_npoints) // transform phys space
446  {
447  if (IsForwards)
448  {
451  }
452  else
453  {
456  }
457  }
458  else
459  {
460  if (IsForwards)
461  {
464  }
465  else
466  {
469  }
470  }
471 
472  int nrowsY = blkmatY->GetRows();
473  int ncolsY = blkmatY->GetColumns();
474 
475  Array<OneD, NekDouble> sortedinarrayY(ncolsY);
476  Array<OneD, NekDouble> sortedoutarrayY(nrowsY);
477 
478  int nrowsZ = blkmatZ->GetRows();
479  int ncolsZ = blkmatZ->GetColumns();
480 
481  Array<OneD, NekDouble> sortedinarrayZ(ncolsZ);
482  Array<OneD, NekDouble> sortedoutarrayZ(nrowsZ);
483 
484  NekVector<NekDouble> inY(ncolsY, sortedinarrayY, eWrapper);
485  NekVector<NekDouble> outY(nrowsY, sortedoutarrayY, eWrapper);
486 
487  NekVector<NekDouble> inZ(ncolsZ, sortedinarrayZ, eWrapper);
488  NekVector<NekDouble> outZ(nrowsZ, sortedoutarrayZ, eWrapper);
489 
490  m_transposition->Transpose(inarray, sortedinarrayY, !IsForwards,
492 
493  outY = (*blkmatY) * inY;
494 
495  m_transposition->Transpose(sortedoutarrayY, sortedinarrayZ, false,
497 
498  outZ = (*blkmatZ) * inZ;
499 
500  m_transposition->Transpose(sortedoutarrayZ, sortedoutarrayY, false,
502 
503  m_transposition->Transpose(sortedoutarrayY, outarray, false,
505  }
506 }
507 
509  Homogeneous2DMatType mattype) const
510 {
511  auto matrixIter = m_homogeneous2DBlockMat->find(mattype);
512 
513  if (matrixIter == m_homogeneous2DBlockMat->end())
514  {
515  return ((*m_homogeneous2DBlockMat)[mattype] =
516  GenHomogeneous2DBlockMatrix(mattype));
517  }
518  else
519  {
520  return matrixIter->second;
521  }
522 }
523 
525  Homogeneous2DMatType mattype) const
526 {
527  int i;
528  int n_exp = 0;
529 
530  DNekMatSharedPtr loc_mat;
531  DNekBlkMatSharedPtr BlkMatrix;
532 
534 
535  int NumPoints = 0;
536  int NumModes = 0;
537  int NumPencils = 0;
538 
539  if ((mattype == eForwardsCoeffSpaceY1D) ||
540  (mattype == eBackwardsCoeffSpaceY1D) ||
541  (mattype == eForwardsPhysSpaceY1D) ||
542  (mattype == eBackwardsPhysSpaceY1D))
543  {
545  NumPoints = m_homogeneousBasis_y->GetNumModes();
546  NumModes = m_homogeneousBasis_y->GetNumPoints();
547  NumPencils = m_homogeneousBasis_z->GetNumPoints();
548  }
549  else
550  {
552  NumPoints = m_homogeneousBasis_z->GetNumModes();
553  NumModes = m_homogeneousBasis_z->GetNumPoints();
554  NumPencils = m_homogeneousBasis_y->GetNumPoints();
555  }
556 
557  if ((mattype == eForwardsCoeffSpaceY1D) ||
558  (mattype == eForwardsCoeffSpaceZ1D) ||
559  (mattype == eBackwardsCoeffSpaceY1D) ||
560  (mattype == eBackwardsCoeffSpaceZ1D))
561  {
562  n_exp = m_lines[0]->GetNcoeffs();
563  }
564  else
565  {
566  n_exp = m_lines[0]->GetTotPoints(); // will operatore on m_phys
567  }
568 
569  Array<OneD, unsigned int> nrows(n_exp);
570  Array<OneD, unsigned int> ncols(n_exp);
571 
572  if ((mattype == eForwardsCoeffSpaceY1D) ||
573  (mattype == eForwardsPhysSpaceY1D) ||
574  (mattype == eForwardsCoeffSpaceZ1D) ||
575  (mattype == eForwardsPhysSpaceZ1D))
576  {
577  nrows = Array<OneD, unsigned int>(n_exp * NumPencils, NumModes);
578  ncols = Array<OneD, unsigned int>(n_exp * NumPencils, NumPoints);
579  }
580  else
581  {
582  nrows = Array<OneD, unsigned int>(n_exp * NumPencils, NumPoints);
583  ncols = Array<OneD, unsigned int>(n_exp * NumPencils, NumModes);
584  }
585 
586  MatrixStorage blkmatStorage = eDIAGONAL;
587  BlkMatrix = MemoryManager<DNekBlkMat>::AllocateSharedPtr(nrows, ncols,
588  blkmatStorage);
589 
591 
592  if ((mattype == eForwardsCoeffSpaceY1D) ||
593  (mattype == eForwardsPhysSpaceY1D) ||
594  (mattype == eForwardsCoeffSpaceZ1D) ||
595  (mattype == eForwardsPhysSpaceZ1D))
596  {
598  StdSeg.DetShapeType(), StdSeg);
599 
600  loc_mat = StdSeg.GetStdMatrix(matkey);
601  }
602  else
603  {
605  StdSeg.DetShapeType(), StdSeg);
606 
607  loc_mat = StdSeg.GetStdMatrix(matkey);
608  }
609 
610  // set up array of block matrices.
611  for (i = 0; i < (n_exp * NumPencils); ++i)
612  {
613  BlkMatrix->SetBlock(i, i, loc_mat);
614  }
615 
616  return BlkMatrix;
617 }
618 
619 std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous2D::
621 {
622  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
623  // Set up Homogeneous length details.
625  HomoBasis[0] = m_homogeneousBasis_y;
626  HomoBasis[1] = m_homogeneousBasis_z;
627 
628  std::vector<NekDouble> HomoLen(2);
629  HomoLen[0] = m_lhom_y;
630  HomoLen[1] = m_lhom_z;
631 
632  int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
633  int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
634 
635  std::vector<unsigned int> sIDs = LibUtilities::NullUnsignedIntVector;
636 
637  std::vector<unsigned int> yIDs;
638  std::vector<unsigned int> zIDs;
639 
640  for (int n = 0; n < nhom_modes_z; ++n)
641  {
642  for (int m = 0; m < nhom_modes_y; ++m)
643  {
644  zIDs.push_back(n);
645  yIDs.push_back(m);
646  }
647  }
648 
649  m_lines[0]->GeneralGetFieldDefinitions(returnval, 2, HomoBasis, HomoLen,
650  false, sIDs, zIDs, yIDs);
651  return returnval;
652 }
653 
655  std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
656 {
657  // Set up Homogeneous length details.
659  HomoBasis[0] = m_homogeneousBasis_y;
660  HomoBasis[1] = m_homogeneousBasis_z;
661  std::vector<NekDouble> HomoLen(2);
662  HomoLen[0] = m_lhom_y;
663  HomoLen[1] = m_lhom_z;
664 
665  int nhom_modes_y = m_homogeneousBasis_y->GetNumModes();
666  int nhom_modes_z = m_homogeneousBasis_z->GetNumModes();
667 
668  std::vector<unsigned int> sIDs = LibUtilities::NullUnsignedIntVector;
669 
670  std::vector<unsigned int> yIDs;
671  std::vector<unsigned int> zIDs;
672 
673  for (int n = 0; n < nhom_modes_z; ++n)
674  {
675  for (int m = 0; m < nhom_modes_y; ++m)
676  {
677  zIDs.push_back(n);
678  yIDs.push_back(m);
679  }
680  }
681 
682  // enforce NumHomoDir == 1 by direct call
683  m_lines[0]->GeneralGetFieldDefinitions(fielddef, 2, HomoBasis, HomoLen,
684  false, sIDs, zIDs, yIDs);
685 }
686 
689  std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
690 {
691  int i, k;
692 
693  int NumMod_y = m_homogeneousBasis_y->GetNumModes();
694  int NumMod_z = m_homogeneousBasis_z->GetNumModes();
695 
696  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
697 
698  // Determine mapping from element ids to location in
699  // expansion list
700  map<int, int> ElmtID_to_ExpID;
701  for (i = 0; i < m_lines[0]->GetExpSize(); ++i)
702  {
703  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
704  }
705 
706  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
707  {
708  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
709  int datalen = (*m_exp)[eid]->GetNcoeffs();
710 
711  for (k = 0; k < (NumMod_y * NumMod_z); ++k)
712  {
713  fielddata.insert(
714  fielddata.end(),
715  &coeffs[m_coeff_offset[eid] + k * ncoeffs_per_line],
716  &coeffs[m_coeff_offset[eid] + k * ncoeffs_per_line] + datalen);
717  }
718  }
719 }
720 
723  std::vector<NekDouble> &fielddata)
724 {
725  v_AppendFieldData(fielddef, fielddata, m_coeffs);
726 }
727 
728 // Extract the data in fielddata into the m_coeff list
731  std::vector<NekDouble> &fielddata, std::string &field,
732  Array<OneD, NekDouble> &coeffs)
733 {
734  int i, k;
735  int offset = 0;
736  int datalen = fielddata.size() / fielddef->m_fields.size();
737  int ncoeffs_per_line = m_lines[0]->GetNcoeffs();
738  int NumMod_y = m_homogeneousBasis_y->GetNumModes();
739  int NumMod_z = m_homogeneousBasis_z->GetNumModes();
740 
741  // Find data location according to field definition
742  for (i = 0; i < fielddef->m_fields.size(); ++i)
743  {
744  if (fielddef->m_fields[i] == field)
745  {
746  break;
747  }
748  offset += datalen;
749  }
750 
751  ASSERTL0(i != fielddef->m_fields.size(), "Field not found in data file");
752 
753  // Determine mapping from element ids to location in
754  // expansion list
755  map<int, int> ElmtID_to_ExpID;
756  for (i = 0; i < m_lines[0]->GetExpSize(); ++i)
757  {
758  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
759  }
760 
761  for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
762  {
763  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
764  int datalen = (*m_exp)[eid]->GetNcoeffs();
765 
766  for (k = 0; k < (NumMod_y * NumMod_z); ++k)
767  {
768  Vmath::Vcopy(datalen, &fielddata[offset], 1,
769  &coeffs[m_coeff_offset[eid] + k * ncoeffs_per_line],
770  1);
771  offset += datalen;
772  }
773  }
774 }
775 
776 void ExpListHomogeneous2D::v_WriteVtkPieceData(std::ostream &outfile,
777  int expansion, std::string var)
778 {
779  int i;
780  int nq = (*m_exp)[expansion]->GetTotPoints();
781  int npoints_per_line = m_lines[0]->GetTotPoints();
782 
783  // printing the fields of that zone
784  outfile << " <DataArray type=\"Float64\" Name=\"" << var << "\">"
785  << endl;
786  outfile << " ";
787  for (int n = 0; n < m_lines.size(); ++n)
788  {
789  const Array<OneD, NekDouble> phys =
790  m_phys + m_phys_offset[expansion] + n * npoints_per_line;
791  for (i = 0; i < nq; ++i)
792  {
793  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i])
794  << " ";
795  }
796  }
797  outfile << endl;
798  outfile << " </DataArray>" << endl;
799 }
800 
804 
805 {
806  int nyzlines = m_lines.size(); // number of Fourier points in the Fourier
807  // directions (nF_pts)
808  int npoints = inarray.size(); // number of total points = n. of Fourier
809  // points * n. of points per line (nT_pts)
810  int n_points_line = npoints / nyzlines; // number of points per line
811 
812  Array<OneD, NekDouble> temparray(npoints);
813  Array<OneD, NekDouble> temparray1(npoints);
814  Array<OneD, NekDouble> temparray2(npoints);
818 
819  for (int i = 0; i < nyzlines; i++)
820  {
821  m_lines[i]->PhysDeriv(tmp1 = inarray + i * n_points_line,
822  tmp2 = out_d0 + i * n_points_line);
823  }
824 
825  if (m_homogeneousBasis_y->GetBasisType() == LibUtilities::eFourier &&
826  m_homogeneousBasis_z->GetBasisType() == LibUtilities::eFourier)
827  {
828  if (m_WaveSpace)
829  {
830  temparray = inarray;
831  }
832  else
833  {
834  HomogeneousFwdTrans(inarray, temparray);
835  }
836  NekDouble sign = -1.0;
837  NekDouble beta;
838 
839  // along y
840  for (int i = 0; i < m_ny; i++)
841  {
842  beta = -sign * 2 * M_PI * (i / 2) / m_lhom_y;
843 
844  for (int j = 0; j < m_nz; j++)
845  {
846  Vmath::Smul(n_points_line, beta,
847  tmp1 = temparray + n_points_line * (i + j * m_ny),
848  1,
849  tmp2 = temparray1 +
850  n_points_line * ((i - int(sign)) + j * m_ny),
851  1);
852  }
853 
854  sign = -1.0 * sign;
855  }
856 
857  // along z
858  sign = -1.0;
859  for (int i = 0; i < m_nz; i++)
860  {
861  beta = -sign * 2 * M_PI * (i / 2) / m_lhom_z;
862  Vmath::Smul(
863  m_ny * n_points_line, beta,
864  tmp1 = temparray + i * m_ny * n_points_line, 1,
865  tmp2 = temparray2 + (i - int(sign)) * m_ny * n_points_line, 1);
866  sign = -1.0 * sign;
867  }
868  if (m_WaveSpace)
869  {
870  out_d1 = temparray1;
871  out_d2 = temparray2;
872  }
873  else
874  {
875  HomogeneousBwdTrans(temparray1, out_d1);
876  HomogeneousBwdTrans(temparray2, out_d2);
877  }
878  }
879  else
880  {
881  if (m_WaveSpace)
882  {
883  ASSERTL0(false, "Semi-phyisical time-stepping not implemented yet "
884  "for non-Fourier basis")
885  }
886  else
887  {
888  StdRegions::StdQuadExp StdQuad(m_homogeneousBasis_y->GetBasisKey(),
889  m_homogeneousBasis_z->GetBasisKey());
890 
891  m_transposition->Transpose(inarray, temparray, false,
893 
894  for (int i = 0; i < n_points_line; i++)
895  {
896  StdQuad.PhysDeriv(tmp1 = temparray + i * nyzlines,
897  tmp2 = temparray1 + i * nyzlines,
898  tmp3 = temparray2 + i * nyzlines);
899  }
900 
901  m_transposition->Transpose(temparray1, out_d1, false,
903  m_transposition->Transpose(temparray2, out_d2, false,
905  Vmath::Smul(npoints, 2.0 / m_lhom_y, out_d1, 1, out_d1, 1);
906  Vmath::Smul(npoints, 2.0 / m_lhom_z, out_d2, 1, out_d2, 1);
907  }
908  }
909 }
910 
912  Direction edir, const Array<OneD, const NekDouble> &inarray,
913  Array<OneD, NekDouble> &out_d)
914 
915 {
916  int nyzlines = m_lines.size(); // number of Fourier points in the Fourier
917  // directions (nF_pts)
918  int npoints = inarray.size(); // number of total points = n. of Fourier
919  // points * n. of points per line (nT_pts)
920  int n_points_line = npoints / nyzlines; // number of points per line
921  // convert enum into int
922  int dir = (int)edir;
923 
924  Array<OneD, NekDouble> temparray(npoints);
925  Array<OneD, NekDouble> temparray1(npoints);
926  Array<OneD, NekDouble> temparray2(npoints);
930 
931  if (dir < 1)
932  {
933  for (int i = 0; i < nyzlines; i++)
934  {
935  m_lines[i]->PhysDeriv(tmp1 = inarray + i * n_points_line,
936  tmp2 = out_d + i * n_points_line);
937  }
938  }
939  else
940  {
941  if (m_homogeneousBasis_y->GetBasisType() == LibUtilities::eFourier &&
942  m_homogeneousBasis_z->GetBasisType() == LibUtilities::eFourier)
943  {
944  if (m_WaveSpace)
945  {
946  temparray = inarray;
947  }
948  else
949  {
950  HomogeneousFwdTrans(inarray, temparray);
951  }
952  NekDouble sign = -1.0;
953  NekDouble beta;
954 
955  if (dir == 1)
956  {
957  // along y
958  for (int i = 0; i < m_ny; i++)
959  {
960  beta = -sign * 2 * M_PI * (i / 2) / m_lhom_y;
961 
962  for (int j = 0; j < m_nz; j++)
963  {
964  Vmath::Smul(
965  n_points_line, beta,
966  tmp1 = temparray + n_points_line * (i + j * m_ny),
967  1,
968  tmp2 = temparray1 +
969  n_points_line * ((i - int(sign)) + j * m_ny),
970  1);
971  }
972  sign = -1.0 * sign;
973  }
974  if (m_WaveSpace)
975  {
976  out_d = temparray1;
977  }
978  else
979  {
980  HomogeneousBwdTrans(temparray1, out_d);
981  }
982  }
983  else
984  {
985  // along z
986  for (int i = 0; i < m_nz; i++)
987  {
988  beta = -sign * 2 * M_PI * (i / 2) / m_lhom_z;
989  Vmath::Smul(m_ny * n_points_line, beta,
990  tmp1 = temparray + i * m_ny * n_points_line, 1,
991  tmp2 = temparray2 +
992  (i - int(sign)) * m_ny * n_points_line,
993  1);
994  sign = -1.0 * sign;
995  }
996  if (m_WaveSpace)
997  {
998  out_d = temparray2;
999  }
1000  else
1001  {
1002  HomogeneousBwdTrans(temparray2, out_d);
1003  }
1004  }
1005  }
1006  else
1007  {
1008  if (m_WaveSpace)
1009  {
1010  ASSERTL0(false, "Semi-phyisical time-stepping not implemented "
1011  "yet for non-Fourier basis")
1012  }
1013  else
1014  {
1015  StdRegions::StdQuadExp StdQuad(
1016  m_homogeneousBasis_y->GetBasisKey(),
1017  m_homogeneousBasis_z->GetBasisKey());
1018 
1019  m_transposition->Transpose(inarray, temparray, false,
1021 
1022  for (int i = 0; i < n_points_line; i++)
1023  {
1024  StdQuad.PhysDeriv(tmp1 = temparray + i * nyzlines,
1025  tmp2 = temparray1 + i * nyzlines,
1026  tmp3 = temparray2 + i * nyzlines);
1027  }
1028 
1029  if (dir == 1)
1030  {
1031  m_transposition->Transpose(temparray1, out_d, false,
1033  Vmath::Smul(npoints, 2.0 / m_lhom_y, out_d, 1, out_d, 1);
1034  }
1035  else
1036  {
1037  m_transposition->Transpose(temparray2, out_d, false,
1039  Vmath::Smul(npoints, 2.0 / m_lhom_z, out_d, 1, out_d, 1);
1040  }
1041  }
1042  }
1043  }
1044 }
1045 
1047  const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &out_d0,
1049 
1050 {
1051  v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
1052 }
1053 
1055  Direction edir, const Array<OneD, const NekDouble> &inarray,
1056  Array<OneD, NekDouble> &out_d)
1057 {
1058  // convert int into enum
1059  v_PhysDeriv(edir, inarray, out_d);
1060 }
1061 
1063 {
1064  NekDouble size_y = 1.5 * m_ny;
1065  NekDouble size_z = 1.5 * m_nz;
1066  m_padsize_y = int(size_y);
1067  m_padsize_z = int(size_z);
1068 
1069  const LibUtilities::PointsKey Ppad_y(m_padsize_y,
1072  Ppad_y);
1073 
1074  const LibUtilities::PointsKey Ppad_z(m_padsize_z,
1077  Ppad_z);
1078 
1081 
1082  StdRegions::StdQuadExp StdQuad(m_paddingBasis_y->GetBasisKey(),
1083  m_paddingBasis_z->GetBasisKey());
1084 
1086  StdQuad.DetShapeType(), StdQuad);
1088  StdQuad.DetShapeType(), StdQuad);
1089 
1090  MatFwdPAD = StdQuad.GetStdMatrix(matkey1);
1091  MatBwdPAD = StdQuad.GetStdMatrix(matkey2);
1092 }
1093 } // namespace MultiRegions
1094 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Represents a basis of a given type.
Definition: Basis.h:213
const BasisKey GetBasisKey() const
Returns the specification for the Basis.
Definition: Basis.h:328
Describes the specification for a Basis.
Definition: Basis.h:50
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:144
Defines a specification for a set of points.
Definition: Points.h:59
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Abstraction of a two-dimensional multi-elemental expansion which is merely a collection of local expa...
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
DNekBlkMatSharedPtr GenHomogeneous2DBlockMatrix(Homogeneous2DMatType mattype) const
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
LibUtilities::NektarFFTSharedPtr m_FFT_y
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
ExpListHomogeneous2D(const ExpansionType type)
Default constructor.
int m_nz
Number of modes = number of poitns in z direction.
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.
LibUtilities::NektarFFTSharedPtr m_FFT_z
LibUtilities::BasisSharedPtr m_homogeneousBasis_y
Definition of the total number of degrees of freedom and quadrature points. Sets up the storage for m...
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
NekDouble m_lhom_z
Width of homogeneous direction z.
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)
void HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
int m_ny
Number of modes = number of poitns in y direction.
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_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void Homogeneous2DTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool IsForwards, bool Shuff=true, bool UnShuff=true)
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
LibUtilities::TranspositionSharedPtr m_transposition
Array< OneD, ExpListSharedPtr > m_lines
Vector of ExpList, will be filled with ExpList1D.
LibUtilities::BasisSharedPtr m_homogeneousBasis_z
Base expansion in z direction.
LibUtilities::BasisSharedPtr m_paddingBasis_z
Base expansion in z direction.
DNekBlkMatSharedPtr GetHomogeneous2DBlockMatrix(Homogeneous2DMatType mattype) const
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble >> &inarray1, const Array< OneD, Array< OneD, NekDouble >> &inarray2, Array< OneD, Array< OneD, NekDouble >> &outarray)
LibUtilities::BasisSharedPtr m_paddingBasis_y
Base expansion in y direction.
NekDouble m_lhom_y
Width of homogeneous direction y.
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
Base class for all multi-elemental spectral/hp expansions.
Definition: ExpList.h:101
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
Definition: ExpList.h:1158
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
Definition: ExpList.h:1210
LibUtilities::CommSharedPtr m_comm
Communicator.
Definition: ExpList.h:1126
void DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
Definition: ExpList.h:2005
LibUtilities::SessionReaderSharedPtr m_session
Session.
Definition: ExpList.h:1129
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
Definition: ExpList.h:1213
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Definition: ExpList.h:1175
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:850
Class representing a segment element in reference space.
Definition: StdSegExp.h:54
static BasisSharedPtr NullBasisSharedPtr
Definition: Basis.h:369
static std::vector< unsigned int > NullUnsignedIntVector
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
NektarFFTFactory & GetNektarFFTFactory()
Definition: NektarFFT.cpp:69
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:177
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:61
@ eFourierEvenlySpaced
1D Evenly-spaced points using Fourier Fit
Definition: PointsType.h:76
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
std::map< Homogeneous2DMatType, DNekBlkMatSharedPtr > Homo2DBlockMatrixMap
A map between homo matrix keys and their associated block matrices.
static const NekDouble kNekZeroTol
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekBlkMat > DNekBlkMatSharedPtr
Definition: NekTypeDefs.hpp:77
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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:209
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255