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