Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ExpListHomogeneous1D.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File ExpListHomogeneous1D.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: An ExpList which is homogeneous in 1-direction
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <MultiRegions/ExpListHomogeneous1D.h>
37 #include <LibUtilities/Foundations/ManagerAccess.h> // for PointsManager, etc
38 #include <StdRegions/StdSegExp.h>
39 #include <StdRegions/StdPointExp.h>
40 #include <LocalRegions/Expansion.h>
41 #include <LocalRegions/Expansion2D.h>
42 
43 namespace Nektar
44 {
45  namespace MultiRegions
46  {
47  // Forward declaration for typedefs
48  ExpListHomogeneous1D::ExpListHomogeneous1D():
49  ExpList(),
50  m_homogeneousBasis(LibUtilities::NullBasisSharedPtr),
51  m_lhom(1),
52  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr())
53  {
54  }
55 
56  ExpListHomogeneous1D::ExpListHomogeneous1D(const LibUtilities::SessionReaderSharedPtr
57  &pSession,const LibUtilities::BasisKey &HomoBasis, const NekDouble lhom, const bool useFFT, const bool dealiasing):
58  ExpList(pSession),
59  m_useFFT(useFFT),
60  m_lhom(lhom),
61  m_homogeneous1DBlockMat(MemoryManager<Homo1DBlockMatrixMap>::AllocateSharedPtr()),
62  m_dealiasing(dealiasing)
63  {
64  ASSERTL2(HomoBasis != LibUtilities::NullBasisKey,"Homogeneous Basis is a null basis");
65 
66  m_homogeneousBasis = LibUtilities::BasisManager()[HomoBasis];
67 
68  m_transposition = MemoryManager<LibUtilities::Transposition>::AllocateSharedPtr(HomoBasis,m_comm->GetColumnComm());
69 
70  m_planes = Array<OneD,ExpListSharedPtr>(m_homogeneousBasis->GetNumPoints()/m_comm->GetColumnComm()->GetSize());
71 
72  if(m_useFFT)
73  {
74  m_FFT = LibUtilities::GetNektarFFTFactory().CreateInstance("NekFFTW", m_homogeneousBasis->GetNumPoints());
75  }
76 
77  if(m_dealiasing)
78  {
79  if(m_useFFT)
80  {
81  NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
82  m_padsize = int(size);
83  m_FFT_deal = LibUtilities::GetNektarFFTFactory()
84  .CreateInstance("NekFFTW", m_padsize);
85  }
86  else
87  {
88  ASSERTL0(false, "Dealiasing available just in combination "
89  "with FFTW");
90  }
91  }
92  }
93 
94 
95  /**
96  * @param In ExpListHomogeneous1D object to copy.
97  */
98  ExpListHomogeneous1D::ExpListHomogeneous1D(const ExpListHomogeneous1D &In):
99  ExpList(In,false),
100  m_transposition(In.m_transposition),
101  m_useFFT(In.m_useFFT),
102  m_FFT(In.m_FFT),
103  m_tmpIN(In.m_tmpIN),
104  m_tmpOUT(In.m_tmpOUT),
105  m_homogeneousBasis(In.m_homogeneousBasis),
106  m_lhom(In.m_lhom),
107  m_homogeneous1DBlockMat(In.m_homogeneous1DBlockMat),
108  m_dealiasing(In.m_dealiasing),
109  m_padsize(In.m_padsize)
110  {
111  m_planes = Array<OneD, ExpListSharedPtr>(In.m_planes.num_elements());
112  }
113 
114  /**
115  * Destructor
116  */
117  ExpListHomogeneous1D::~ExpListHomogeneous1D()
118  {
119  }
120 
121  void ExpListHomogeneous1D::v_HomogeneousFwdTrans(const Array<OneD, const NekDouble> &inarray,
122  Array<OneD, NekDouble> &outarray,
123  CoeffState coeffstate,
124  bool Shuff,
125  bool UnShuff)
126  {
127  // Forwards trans
128  Homogeneous1DTrans(inarray,outarray,true,coeffstate,Shuff,UnShuff);
129  }
130 
131  void ExpListHomogeneous1D::v_HomogeneousBwdTrans(const Array<OneD, const NekDouble> &inarray,
132  Array<OneD, NekDouble> &outarray,
133  CoeffState coeffstate,
134  bool Shuff,
135  bool UnShuff)
136  {
137  // Backwards trans
138  Homogeneous1DTrans(inarray,outarray,false,coeffstate,Shuff,UnShuff);
139  }
140 
141  /**
142  * Dealiasing routine
143  */
144  void ExpListHomogeneous1D::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,
145  const Array<OneD, NekDouble> &inarray2,
146  Array<OneD, NekDouble> &outarray,
147  CoeffState coeffstate)
148  {
149  // inarray1 = first term of the product in full physical space
150  // inarray2 = second term of the product in full physical space
151  // dealiased product stored in outarray
152 
153  int num_dofs = inarray1.num_elements();
154 
155  int N = m_homogeneousBasis->GetNumPoints();
156 
157  Array<OneD, NekDouble> V1(num_dofs);
158  Array<OneD, NekDouble> V2(num_dofs);
159  Array<OneD, NekDouble> V1V2(num_dofs);
160 
161  HomogeneousFwdTrans(inarray1,V1,coeffstate);
162  HomogeneousFwdTrans(inarray2,V2,coeffstate);
163 
164  int num_points_per_plane = num_dofs/m_planes.num_elements();
165  int num_proc = m_comm->GetColumnComm()->GetSize();
166  int num_dfts_per_proc = num_points_per_plane / num_proc
167  + (num_points_per_plane % num_proc > 0);
168 
169  Array<OneD, NekDouble> ShufV1(num_dfts_per_proc*N,0.0);
170  Array<OneD, NekDouble> ShufV2(num_dfts_per_proc*N,0.0);
171  Array<OneD, NekDouble> ShufV1V2(num_dfts_per_proc*N,0.0);
172 
173  Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
174  Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
175  Array<OneD, NekDouble> ShufV1_PAD_phys(m_padsize,0.0);
176  Array<OneD, NekDouble> ShufV2_PAD_phys(m_padsize,0.0);
177 
178  Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
179  Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
180 
181  m_transposition->Transpose(V1, ShufV1, false, LibUtilities::eXYtoZ);
182  m_transposition->Transpose(V2, ShufV2, false, LibUtilities::eXYtoZ);
183 
184  // Looping on the pencils
185  for(int i = 0 ; i < num_dfts_per_proc ; i++)
186  {
187  // Copying the i-th pencil pf lenght N into a bigger
188  // pencil of lenght 2N We are in Fourier space
189  Vmath::Vcopy(N, &(ShufV1[i*N]), 1, &(ShufV1_PAD_coef[0]), 1);
190  Vmath::Vcopy(N, &(ShufV2[i*N]), 1, &(ShufV2_PAD_coef[0]), 1);
191 
192  // Moving to physical space using the padded system
193  m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys);
194  m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys);
195 
196  // Perfroming the vectors multiplication in physical space on
197  // the padded system
198  Vmath::Vmul(m_padsize, ShufV1_PAD_phys, 1,
199  ShufV2_PAD_phys, 1,
200  ShufV1V2_PAD_phys, 1);
201 
202  // Moving back the result (V1*V2)_phys in Fourier space, padded
203  // system
204  m_FFT_deal->FFTFwdTrans(ShufV1V2_PAD_phys, ShufV1V2_PAD_coef);
205 
206  // Copying the first half of the padded pencil in the full
207  // vector (Fourier space)
208  Vmath::Vcopy(N, &(ShufV1V2_PAD_coef[0]), 1,
209  &(ShufV1V2[i*N]), 1);
210  }
211 
212  m_transposition->Transpose(ShufV1V2, V1V2, false,
213  LibUtilities::eZtoXY);
214 
215  // Moving the results in physical space for the output
216  HomogeneousBwdTrans(V1V2, outarray, coeffstate);
217  }
218 
219  /**
220  * Forward transform
221  */
222  void ExpListHomogeneous1D::v_FwdTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate )
223  {
224  int cnt = 0, cnt1 = 0;
225  Array<OneD, NekDouble> tmparray;
226 
227  for(int n = 0; n < m_planes.num_elements(); ++n)
228  {
229  m_planes[n]->FwdTrans(inarray+cnt, tmparray = outarray + cnt1,
230  coeffstate);
231  cnt += m_planes[n]->GetTotPoints();
232 
233  cnt1 += m_planes[n]->GetNcoeffs(); // need to skip ncoeffs
234  }
235  if(!m_WaveSpace)
236  {
237  HomogeneousFwdTrans(outarray,outarray,coeffstate);
238  }
239  }
240 
241  /**
242  * Forward transform element by element
243  */
244  void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
245  {
246  int cnt = 0, cnt1 = 0;
247  Array<OneD, NekDouble> tmparray;
248 
249  //spectral element FwdTrans plane by plane
250  for(int n = 0; n < m_planes.num_elements(); ++n)
251  {
252  m_planes[n]->FwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
253  cnt += m_planes[n]->GetTotPoints();
254  cnt1 += m_planes[n]->GetNcoeffs();
255  }
256  if(!m_WaveSpace)
257  {
258  HomogeneousFwdTrans(outarray,outarray);
259  }
260  }
261 
262  /**
263  * Backward transform
264  */
265  void ExpListHomogeneous1D::v_BwdTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate)
266  {
267  int cnt = 0, cnt1 = 0;
268  Array<OneD, NekDouble> tmparray;
269 
270  for(int n = 0; n < m_planes.num_elements(); ++n)
271  {
272  m_planes[n]->BwdTrans(inarray+cnt, tmparray = outarray + cnt1,
273  coeffstate);
274  cnt += m_planes[n]->GetNcoeffs();
275  cnt1 += m_planes[n]->GetTotPoints();
276  }
277  if(!m_WaveSpace)
278  {
279  HomogeneousBwdTrans(outarray,outarray);
280  }
281  }
282 
283  /**
284  * Backward transform element by element
285  */
286  void ExpListHomogeneous1D::v_BwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
287  {
288  int cnt = 0, cnt1 = 0;
289  Array<OneD, NekDouble> tmparray;
290 
291  for(int n = 0; n < m_planes.num_elements(); ++n)
292  {
293  m_planes[n]->BwdTrans_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
294 
295  cnt += m_planes[n]->GetNcoeffs();
296  cnt1 += m_planes[n]->GetTotPoints();
297  }
298  if(!m_WaveSpace)
299  {
300  HomogeneousBwdTrans(outarray,outarray);
301  }
302  }
303 
304  /**
305  * Inner product
306  */
307  void ExpListHomogeneous1D::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray, CoeffState coeffstate)
308  {
309  int cnt = 0, cnt1 = 0;
310  Array<OneD, NekDouble> tmparray;
311 
312  for(int n = 0; n < m_planes.num_elements(); ++n)
313  {
314  m_planes[n]->IProductWRTBase(inarray+cnt, tmparray = outarray + cnt1,coeffstate);
315 
316  cnt1 += m_planes[n]->GetNcoeffs();
317  cnt += m_planes[n]->GetTotPoints();
318  }
319  }
320 
321  /**
322  * Inner product element by element
323  */
324  void ExpListHomogeneous1D::v_IProductWRTBase_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
325  {
326  int cnt = 0, cnt1 = 0;
327  Array<OneD, NekDouble> tmparray;
328 
329  for(int n = 0; n < m_planes.num_elements(); ++n)
330  {
331  m_planes[n]->IProductWRTBase_IterPerExp(inarray+cnt, tmparray = outarray + cnt1);
332 
333  cnt1 += m_planes[n]->GetNcoeffs();
334  cnt += m_planes[n]->GetTotPoints();
335  }
336  }
337 
338  /**
339  * Homogeneous transform Bwd/Fwd (MVM and FFT)
340  */
341  void ExpListHomogeneous1D::Homogeneous1DTrans(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray,
342  bool IsForwards,
343  CoeffState coeffstate,
344  bool Shuff,
345  bool UnShuff)
346  {
347  int num_dofs;
348 
349  if(IsForwards)
350  {
351  num_dofs = inarray.num_elements();
352  }
353  else
354  {
355  num_dofs = outarray.num_elements();
356  }
357 
358  if(m_useFFT)
359  {
360 
361  int num_points_per_plane = num_dofs/m_planes.num_elements();
362  int num_dfts_per_proc = num_points_per_plane/m_comm->GetColumnComm()->GetSize() + (num_points_per_plane%m_comm->GetColumnComm()->GetSize() > 0);
363 
364  Array<OneD, NekDouble> fft_in (num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
365  Array<OneD, NekDouble> fft_out(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),0.0);
366 
367  if(Shuff)
368  {
369  m_transposition->Transpose(inarray,fft_in,false,LibUtilities::eXYtoZ);
370  }
371  else
372  {
373  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
374  inarray,1,fft_in,1);
375  }
376 
377  if(IsForwards)
378  {
379  for(int i = 0 ; i < num_dfts_per_proc ; i++)
380  {
381  m_FFT->FFTFwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
382  }
383  }
384  else
385  {
386  for(int i = 0 ; i < num_dfts_per_proc ; i++)
387  {
388  m_FFT->FFTBwdTrans(m_tmpIN = fft_in + i*m_homogeneousBasis->GetNumPoints(), m_tmpOUT = fft_out + i*m_homogeneousBasis->GetNumPoints());
389  }
390  }
391 
392  if(UnShuff)
393  {
394  m_transposition->Transpose(fft_out,outarray,false,LibUtilities::eZtoXY);
395  }
396  else
397  {
398  Vmath::Vcopy(num_dfts_per_proc*m_homogeneousBasis->GetNumPoints(),
399  fft_out,1,outarray,1);
400  }
401  }
402  else
403  {
404  DNekBlkMatSharedPtr blkmat;
405 
406  if(num_dofs == m_npoints) //transform phys space
407  {
408  if(IsForwards)
409  {
410  blkmat = GetHomogeneous1DBlockMatrix(eForwardsPhysSpace1D);
411  }
412  else
413  {
414  blkmat = GetHomogeneous1DBlockMatrix(eBackwardsPhysSpace1D);
415  }
416  }
417  else
418  {
419  if(IsForwards)
420  {
421  blkmat = GetHomogeneous1DBlockMatrix(eForwardsCoeffSpace1D,coeffstate);
422  }
423  else
424  {
425  blkmat = GetHomogeneous1DBlockMatrix(eBackwardsCoeffSpace1D,coeffstate);
426  }
427  }
428 
429  int nrows = blkmat->GetRows();
430  int ncols = blkmat->GetColumns();
431 
432  Array<OneD, NekDouble> sortedinarray(ncols,0.0);
433  Array<OneD, NekDouble> sortedoutarray(nrows,0.0);
434 
435  if(Shuff)
436  {
437  m_transposition->Transpose(inarray,sortedinarray,!IsForwards,LibUtilities::eXYtoZ);
438  }
439  else
440  {
441  Vmath::Vcopy(ncols,inarray,1,sortedinarray,1);
442  }
443 
444  // Create NekVectors from the given data arrays
445  NekVector<NekDouble> in (ncols,sortedinarray,eWrapper);
446  NekVector<NekDouble> out(nrows,sortedoutarray,eWrapper);
447 
448  // Perform matrix-vector multiply.
449  out = (*blkmat)*in;
450 
451  if(UnShuff)
452  {
453  m_transposition->Transpose(sortedoutarray,outarray,IsForwards,LibUtilities::eZtoXY);
454  }
455  else
456  {
457  Vmath::Vcopy(nrows,sortedinarray,1,outarray,1);
458  }
459 
460  }
461  }
462 
463  DNekBlkMatSharedPtr ExpListHomogeneous1D::GetHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate) const
464  {
465  Homo1DBlockMatrixMap::iterator matrixIter = m_homogeneous1DBlockMat->find(mattype);
466 
467  if(matrixIter == m_homogeneous1DBlockMat->end())
468  {
469  return ((*m_homogeneous1DBlockMat)[mattype] =
470  GenHomogeneous1DBlockMatrix(mattype,coeffstate));
471  }
472  else
473  {
474  return matrixIter->second;
475  }
476  }
477 
478 
479  DNekBlkMatSharedPtr ExpListHomogeneous1D::GenHomogeneous1DBlockMatrix(Homogeneous1DMatType mattype, CoeffState coeffstate) const
480  {
481  DNekMatSharedPtr loc_mat;
482  DNekBlkMatSharedPtr BlkMatrix;
483  int n_exp = 0;
484  int num_trans_per_proc = 0;
485 
486 
487  if((mattype == eForwardsCoeffSpace1D)
488  ||(mattype == eBackwardsCoeffSpace1D)) // will operate on m_coeffs
489  {
490  n_exp = m_planes[0]->GetNcoeffs();
491  }
492  else
493  {
494  n_exp = m_planes[0]->GetTotPoints(); // will operatore on m_phys
495  }
496 
497  num_trans_per_proc = n_exp/m_comm->GetColumnComm()->GetSize() + (n_exp%m_comm->GetColumnComm()->GetSize() > 0);
498 
499  Array<OneD,unsigned int> nrows(num_trans_per_proc);
500  Array<OneD,unsigned int> ncols(num_trans_per_proc);
501 
502  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
503  {
504  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
505  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
506  }
507  else
508  {
509  nrows = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumPoints());
510  ncols = Array<OneD, unsigned int>(num_trans_per_proc,m_homogeneousBasis->GetNumModes());
511  }
512 
513  MatrixStorage blkmatStorage = eDIAGONAL;
514  BlkMatrix = MemoryManager<DNekBlkMat>
515  ::AllocateSharedPtr(nrows,ncols,blkmatStorage);
516 
517  //Half Mode
518  if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
519  {
520  StdRegions::StdPointExp StdPoint(m_homogeneousBasis->GetBasisKey());
521 
522  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
523  {
524  StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
525  StdPoint.DetShapeType(),
526  StdPoint);
527 
528  loc_mat = StdPoint.GetStdMatrix(matkey);
529  }
530  else
531  {
532  StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
533  StdPoint.DetShapeType(),
534  StdPoint);
535 
536  loc_mat = StdPoint.GetStdMatrix(matkey);
537  }
538  }
539  //other cases
540  else
541  {
542  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
543 
544  if((mattype == eForwardsCoeffSpace1D)||(mattype == eForwardsPhysSpace1D))
545  {
546  StdRegions::StdMatrixKey matkey(StdRegions::eFwdTrans,
547  StdSeg.DetShapeType(),
548  StdSeg);
549 
550  loc_mat = StdSeg.GetStdMatrix(matkey);
551  }
552  else
553  {
554  StdRegions::StdMatrixKey matkey(StdRegions::eBwdTrans,
555  StdSeg.DetShapeType(),
556  StdSeg);
557 
558  loc_mat = StdSeg.GetStdMatrix(matkey);
559  }
560 
561  }
562 
563  // set up array of block matrices.
564  for(int i = 0; i < num_trans_per_proc; ++i)
565  {
566  BlkMatrix->SetBlock(i,i,loc_mat);
567  }
568 
569  return BlkMatrix;
570  }
571 
572  std::vector<LibUtilities::FieldDefinitionsSharedPtr> ExpListHomogeneous1D::v_GetFieldDefinitions()
573  {
574  std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
575 
576  // Set up Homogeneous length details.
577  Array<OneD,LibUtilities::BasisSharedPtr> HomoBasis(1,m_homogeneousBasis);
578 
579  std::vector<NekDouble> HomoLen;
580  HomoLen.push_back(m_lhom);
581 
582  std::vector<unsigned int> PlanesIDs;
583 
584  for(int i = 0; i < m_planes.num_elements(); i++)
585  {
586  PlanesIDs.push_back(m_transposition->GetPlaneID(i));
587  }
588 
589  m_planes[0]->GeneralGetFieldDefinitions(returnval, 1, HomoBasis, HomoLen, PlanesIDs);
590 
591  return returnval;
592  }
593 
594  void ExpListHomogeneous1D::v_GetFieldDefinitions(std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
595  {
596  // Set up Homogeneous length details.
597  Array<OneD,LibUtilities::BasisSharedPtr> HomoBasis(1,m_homogeneousBasis);
598 
599  std::vector<NekDouble> HomoLen;
600  HomoLen.push_back(m_lhom);
601 
602  std::vector<unsigned int> PlanesIDs;
603 
604  for(int i = 0; i < m_planes.num_elements(); i++)
605  {
606  PlanesIDs.push_back(m_transposition->GetPlaneID(i));
607  }
608 
609  // enforce NumHomoDir == 1 by direct call
610  m_planes[0]->GeneralGetFieldDefinitions(fielddef,1, HomoBasis,HomoLen,PlanesIDs);
611  }
612 
613 
614  /** This routine appends the data from the expansion list into
615  the output format where each element is given by looping
616  over its Fourier modes where as data in the expandion is
617  stored with all consecutive elements and then the Fourier
618  modes
619  */
620  void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata, Array<OneD, NekDouble> &coeffs)
621  {
622  int i,n;
623  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
624 
625  // Determine mapping from element ids to location in
626  // expansion list
627  map<int, int> ElmtID_to_ExpID;
628  for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
629  {
630  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
631  }
632 
633  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
634  {
635  int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
636  int datalen = (*m_exp)[eid]->GetNcoeffs();
637 
638  for(n = 0; n < m_planes.num_elements(); ++n)
639  {
640  fielddata.insert(fielddata.end(),&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane],&coeffs[m_coeff_offset[eid]+n*ncoeffs_per_plane]+datalen);
641  }
642  }
643  }
644 
645  void ExpListHomogeneous1D::v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector<NekDouble> &fielddata)
646  {
647  v_AppendFieldData(fielddef,fielddata,m_coeffs);
648  }
649 
650  //Extract the data in fielddata into the m_coeff list
651  void ExpListHomogeneous1D::v_ExtractDataToCoeffs(
652  LibUtilities::FieldDefinitionsSharedPtr &fielddef,
653  std::vector<NekDouble> &fielddata,
654  std::string &field,
655  Array<OneD, NekDouble> &coeffs)
656  {
657  int i,n;
658  int offset = 0;
659  int nzmodes;
660  int datalen = fielddata.size()/fielddef->m_fields.size();
661  std::vector<unsigned int> fieldDefHomoZids;
662 
663 
664  // Find data location according to field definition
665  for(i = 0; i < fielddef->m_fields.size(); ++i)
666  {
667  if(fielddef->m_fields[i] == field)
668  {
669  break;
670  }
671  offset += datalen;
672  }
673 
674  if(i == fielddef->m_fields.size())
675  {
676  cout << "Field "<< field<< "not found in data file. " << endl;
677  }
678  else
679  {
680 
681  int modes_offset = 0;
682  int planes_offset = 0;
683  Array<OneD, NekDouble> coeff_tmp;
684  std::map<int,int>::iterator it;
685 
686 
687  // Build map of plane IDs lying on this processor.
688  std::map<int,int> homoZids;
689  for (i = 0; i < m_planes.num_elements(); ++i)
690  {
691  homoZids[m_transposition->GetPlaneID(i)] = i;
692  }
693 
694  if(fielddef->m_numHomogeneousDir)
695  {
696  for(i = 0; i < fielddef->m_basis.size(); ++i)
697  {
698  if(fielddef->m_basis[i] == m_homogeneousBasis->GetBasisType())
699  {
700  nzmodes = fielddef->m_homogeneousZIDs.size();
701  break;
702  }
703  }
704  ASSERTL1(i != fielddef->m_basis.size(),"Failed to determine number of Homogeneous modes");
705 
706  fieldDefHomoZids = fielddef->m_homogeneousZIDs;
707  }
708  else // input file is 2D and so set nzmodes to 1
709  {
710  nzmodes = 1;
711  fieldDefHomoZids.push_back(0);
712  }
713 
714  // Determine mapping from element ids to location in expansion list.
715  map<int, int> ElmtID_to_ExpID;
716  for(i = 0; i < m_planes[0]->GetExpSize(); ++i)
717  {
718  ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
719  }
720 
721 
722  // calculate number of modes in the current partition
723  int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
724 
725  for(i = 0; i < fielddef->m_elementIDs.size(); ++i)
726  {
727  if(fielddef->m_uniOrder == true) // reset modes_offset to zero
728  {
729  modes_offset = 0;
730  }
731 
732  int datalen = LibUtilities::GetNumberOfCoefficients(fielddef->m_shapeType,
733  fielddef->m_numModes,
734  modes_offset);
735 
736  it = ElmtID_to_ExpID.find(fielddef->m_elementIDs[i]);
737 
738  // ensure element is on this partition for parallel case.
739  if(it == ElmtID_to_ExpID.end())
740  {
741  // increase offset for correct FieldData access
742  offset += datalen*nzmodes;
743  continue;
744  }
745 
746  int eid = it->second;
747 
748 
749  for(n = 0; n < nzmodes; ++n, offset += datalen)
750  {
751 
752  it = homoZids.find(fieldDefHomoZids[n]);
753 
754  // Check to make sure this mode number lies in this field.
755  if (it == homoZids.end())
756  {
757  continue;
758  }
759 
760  planes_offset = it->second;
761  if(datalen == (*m_exp)[eid]->GetNcoeffs())
762  {
763  Vmath::Vcopy(datalen,&fielddata[offset],1,&coeffs[m_coeff_offset[eid]+planes_offset*ncoeffs_per_plane],1);
764  }
765  else // unpack data to new order
766  {
767  (*m_exp)[eid]->ExtractDataToCoeffs(&fielddata[offset], fielddef->m_numModes,modes_offset,&coeffs[m_coeff_offset[eid] + planes_offset*ncoeffs_per_plane]);
768  }
769  }
770  modes_offset += (*m_exp)[0]->GetNumBases();
771  }
772  }
773  }
774 
775  //Extract the data in fielddata into the m_coeff list
776  void ExpListHomogeneous1D::v_ExtractCoeffsToCoeffs(
777  const boost::shared_ptr<ExpList> &fromExpList,const Array<OneD, const NekDouble> &fromCoeffs, Array<OneD, NekDouble> &toCoeffs)
778  {
779  int i;
780  int fromNcoeffs_per_plane = fromExpList->GetPlane(0)->GetNcoeffs();
781  Array<OneD, NekDouble> tocoeffs_tmp, fromcoeffs_tmp;
782 
783  for(i = 0; i < m_planes.num_elements(); ++i)
784  {
785  m_planes[i]->ExtractCoeffsToCoeffs(fromExpList->GetPlane(i),fromcoeffs_tmp = fromCoeffs + fromNcoeffs_per_plane*i, tocoeffs_tmp = toCoeffs + m_ncoeffs*i);
786  }
787  }
788 
789  void ExpListHomogeneous1D::v_WriteVtkPieceData(std::ofstream &outfile, int expansion,
790  std::string var)
791  {
792  int i;
793  int nq = (*m_exp)[expansion]->GetTotPoints();
794  int npoints_per_plane = m_planes[0]->GetTotPoints();
795 
796  // printing the fields of that zone
797  outfile << " <DataArray type=\"Float32\" Name=\""
798  << var << "\">" << endl;
799  outfile << " ";
800  for (int n = 0; n < m_planes.num_elements(); ++n)
801  {
802  const Array<OneD, NekDouble> phys = m_phys + m_phys_offset[expansion] + n*npoints_per_plane;
803  for(i = 0; i < nq; ++i)
804  {
805  outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
806  }
807  }
808  outfile << endl;
809  outfile << " </DataArray>" << endl;
810  }
811 
812  void ExpListHomogeneous1D::v_PhysInterp1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
813  {
814  int cnt,cnt1;
815  Array<OneD, NekDouble> tmparray;
816  cnt = m_planes[0]->GetTotPoints();
817  cnt1 = m_planes[0]->Get1DScaledTotPoints(scale);
818 
819  ASSERTL1(m_planes.num_elements()*cnt1 <= outarray.num_elements(),"size of outarray does not match internal estimage");
820 
821 
822  for(int i = 0; i < m_planes.num_elements(); i++)
823  {
824 
825  m_planes[i]->PhysInterp1DScaled(scale,inarray+i*cnt,
826  tmparray = outarray+i*cnt1);
827  }
828  }
829 
830 
831  void ExpListHomogeneous1D::v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
832  {
833  int cnt,cnt1;
834  Array<OneD, NekDouble> tmparray;
835  cnt = m_planes[0]->Get1DScaledTotPoints(scale);
836  cnt1 = m_planes[0]->GetTotPoints();
837 
838  ASSERTL1(m_planes.num_elements()*cnt <= inarray.num_elements(),"size of outarray does not match internal estimage");
839 
840 
841  for(int i = 0; i < m_planes.num_elements(); i++)
842  {
843  m_planes[i]->PhysGalerkinProjection1DScaled(scale,inarray+i*cnt,
844  tmparray = outarray+i*cnt1);
845  }
846 
847  }
848  void ExpListHomogeneous1D::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
849  Array<OneD, NekDouble> &out_d0,
850  Array<OneD, NekDouble> &out_d1,
851  Array<OneD, NekDouble> &out_d2)
852  {
853  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
854  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
855 
856  Array<OneD, NekDouble> temparray(nT_pts);
857  Array<OneD, NekDouble> outarray(nT_pts);
858  Array<OneD, NekDouble> tmp1;
859  Array<OneD, NekDouble> tmp2;
860  Array<OneD, NekDouble> tmp3;
861 
862  for(int i = 0; i < m_planes.num_elements(); i++)
863  {
864  m_planes[i]->PhysDeriv(inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
865  }
866 
867  if(out_d2 != NullNekDouble1DArray)
868  {
869  if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
870  m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
871  {
872  if(m_WaveSpace)
873  {
874  temparray = inarray;
875  }
876  else
877  {
878  HomogeneousFwdTrans(inarray,temparray);
879  }
880 
881  NekDouble sign = -1.0;
882  NekDouble beta;
883 
884  //Half Mode
885  if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
886  {
887  beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
888 
889  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
890  }
891  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
892  {
893  beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
894 
895  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
896  }
897 
898  //Fully complex
899  else
900  {
901  for(int i = 0; i < m_planes.num_elements(); i++)
902  {
903  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
904 
905  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
906 
907  sign = -1.0*sign;
908  }
909  }
910 
911  if(m_WaveSpace)
912  {
913  out_d2 = outarray;
914  }
915  else
916  {
917  HomogeneousBwdTrans(outarray,out_d2);
918  }
919  }
920  else
921  {
922  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Parallelisation in the homogeneous direction implemented just for Fourier basis");
923 
924  if(m_WaveSpace)
925  {
926 
927  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
928  }
929  else
930  {
931  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
932 
933  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
934 
935  for(int i = 0; i < nP_pts; i++)
936  {
937  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
938  }
939 
940  m_transposition->Transpose(outarray,out_d2,false,LibUtilities::eZtoXY);
941 
942  Vmath::Smul(nT_pts,2.0/m_lhom,out_d2,1,out_d2,1);
943  }
944  }
945  }
946  }
947 
948  void ExpListHomogeneous1D::v_PhysDeriv(Direction edir,
949  const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &out_d)
950 
951  {
952  int nT_pts = inarray.num_elements(); //number of total points = n. of Fourier points * n. of points per plane (nT_pts)
953  int nP_pts = nT_pts/m_planes.num_elements(); //number of points per plane = n of Fourier transform required (nP_pts)
954 
955  int dir= (int)edir;
956 
957  Array<OneD, NekDouble> temparray(nT_pts);
958  Array<OneD, NekDouble> outarray(nT_pts);
959  Array<OneD, NekDouble> tmp1;
960  Array<OneD, NekDouble> tmp2;
961 
962  if (dir < 2)
963  {
964  for(int i=0; i < m_planes.num_elements(); i++)
965  {
966  m_planes[i]->PhysDeriv(edir, inarray + i*nP_pts ,tmp2 = out_d + i*nP_pts);
967  }
968  }
969  else
970  {
971  if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
972  m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
973  {
974  if(m_WaveSpace)
975  {
976  temparray = inarray;
977  }
978  else
979  {
980  HomogeneousFwdTrans(inarray,temparray);
981  }
982 
983  NekDouble sign = -1.0;
984  NekDouble beta;
985 
986  //HalfMode
987  if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
988  {
989  beta = 2*M_PI*(m_transposition->GetK(0))/m_lhom;
990 
991  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
992  }
993  else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
994  {
995  beta = -2*M_PI*(m_transposition->GetK(0))/m_lhom;
996 
997  Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
998  }
999  //Fully complex
1000  else
1001  {
1002  for(int i = 0; i < m_planes.num_elements(); i++)
1003  {
1004  beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
1005 
1006  Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
1007 
1008  sign = -1.0*sign;
1009  }
1010  }
1011  if(m_WaveSpace)
1012  {
1013  out_d = outarray;
1014  }
1015  else
1016  {
1017  HomogeneousBwdTrans(outarray,out_d);
1018  }
1019  }
1020  else
1021  {
1022  ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Parallelisation in the homogeneous direction implemented just for Fourier basis");
1023 
1024  if(m_WaveSpace)
1025  {
1026  ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
1027  }
1028  else
1029  {
1030  StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
1031 
1032  m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
1033 
1034  for(int i = 0; i < nP_pts; i++)
1035  {
1036  StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
1037  }
1038 
1039  m_transposition->Transpose(outarray,out_d,false,LibUtilities::eZtoXY);
1040 
1041  Vmath::Smul(nT_pts,2.0/m_lhom,out_d,1,out_d,1);
1042  }
1043  }
1044  }
1045  }
1046 
1047  void ExpListHomogeneous1D::PhysDeriv(const Array<OneD, const NekDouble> &inarray,
1048  Array<OneD, NekDouble> &out_d0,
1049  Array<OneD, NekDouble> &out_d1,
1050  Array<OneD, NekDouble> &out_d2)
1051 
1052  {
1053  v_PhysDeriv(inarray,out_d0,out_d1,out_d2);
1054  }
1055 
1056  void ExpListHomogeneous1D::PhysDeriv(Direction edir,
1057  const Array<OneD, const NekDouble> &inarray,
1058  Array<OneD, NekDouble> &out_d)
1059  {
1060  v_PhysDeriv(edir,inarray,out_d);
1061  }
1062 
1063  LibUtilities::TranspositionSharedPtr ExpListHomogeneous1D::v_GetTransposition(void)
1064  {
1065  return m_transposition;
1066  }
1067 
1068  NekDouble ExpListHomogeneous1D::v_GetHomoLen(void)
1069  {
1070  return m_lhom;
1071  }
1072 
1073  Array<OneD, const unsigned int> ExpListHomogeneous1D::v_GetZIDs(void)
1074  {
1075  return m_transposition->GetPlanesIDs();
1076  }
1077  } //end of namespace
1078 } //end of namespace
1079 
1080 
1081 /**
1082 * $Log: v $
1083 *
1084 **/
1085 
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2