Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
BwdTrans.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: BwdTrans.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: BwdTrans operator implementations
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <loki/Singleton.h>
37 #include <Collections/Operator.h>
38 #include <Collections/Collection.h>
39 
40 using namespace std;
41 
42 namespace Nektar {
43 namespace Collections {
44 
52 
53 /**
54  * @brief Backward transform operator using standard matrix approach.
55  */
56 class BwdTrans_StdMat : public Operator
57 {
58  public:
60 
61  virtual ~BwdTrans_StdMat()
62  {
63  }
64 
65  virtual void operator()(
66  const Array<OneD, const NekDouble> &input,
67  Array<OneD, NekDouble> &output,
68  Array<OneD, NekDouble> &output1,
69  Array<OneD, NekDouble> &output2,
71  {
72  Blas::Dgemm('N', 'N', m_mat->GetRows(), m_numElmt,
73  m_mat->GetColumns(), 1.0, m_mat->GetRawPtr(),
74  m_mat->GetRows(), input.get(), m_stdExp->GetNcoeffs(),
75  0.0, output.get(), m_stdExp->GetTotPoints());
76  }
77 
78  virtual void operator()(
79  int dir,
80  const Array<OneD, const NekDouble> &input,
81  Array<OneD, NekDouble> &output,
83  {
84  ASSERTL0(false, "Not valid for this operator.");
85  }
86 
87  protected:
89 
90  private:
92  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
94  : Operator(pCollExp, pGeomData)
95  {
97  m_stdExp->DetShapeType(), *m_stdExp);
98  m_mat = m_stdExp->GetStdMatrix(key);
99  }
100 };
101 
102 /// Factory initialisation for the BwdTrans_StdMat operators
103 OperatorKey BwdTrans_StdMat::m_typeArr[] = {
106  BwdTrans_StdMat::create, "BwdTrans_StdMat_Seg"),
109  BwdTrans_StdMat::create, "BwdTrans_StdMat_Tri"),
112  BwdTrans_StdMat::create, "BwdTrans_StdMat_NodalTri"),
115  BwdTrans_StdMat::create, "BwdTrans_StdMat_Quad"),
118  BwdTrans_StdMat::create, "BwdTrans_StdMat_Tet"),
121  BwdTrans_StdMat::create, "BwdTrans_StdMat_NodalTet"),
124  BwdTrans_StdMat::create, "BwdTrans_StdMat_Pyr"),
127  BwdTrans_StdMat::create, "BwdTrans_StdMat_Prism"),
130  BwdTrans_StdMat::create, "BwdTrans_StdMat_NodalPrism"),
133  BwdTrans_StdMat::create, "BwdTrans_StdMat_Hex"),
134 };
135 
136 
137 
138 /**
139  * @brief Backward transform operator using default StdRegions operator
140  */
142 {
143  public:
145 
147  {
148  }
149 
150  virtual void operator()(
151  const Array<OneD, const NekDouble> &input,
152  Array<OneD, NekDouble> &output,
153  Array<OneD, NekDouble> &output1,
154  Array<OneD, NekDouble> &output2,
156  {
157  const int nCoeffs = m_stdExp->GetNcoeffs();
158  const int nPhys = m_stdExp->GetTotPoints();
160 
161  for (int i = 0; i < m_numElmt; ++i)
162  {
163  m_stdExp->BwdTrans(input + i*nCoeffs, tmp = output + i*nPhys);
164  }
165  }
166 
167  virtual void operator()(
168  int dir,
169  const Array<OneD, const NekDouble> &input,
170  Array<OneD, NekDouble> &output,
172  {
173  ASSERTL0(false, "Not valid for this operator.");
174  }
175 
176  private:
178  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
179  CoalescedGeomDataSharedPtr pGeomData)
180  : Operator(pCollExp, pGeomData)
181  {
182  }
183 };
184 
185 /// Factory initialisation for the BwdTrans_IterPerExp operators
186 OperatorKey BwdTrans_IterPerExp::m_typeArr[] = {
189  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Seg"),
192  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Tri"),
195  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_NodalTri"),
198  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Quad"),
201  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Tet"),
204  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_NodalTet"),
207  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Pyr"),
210  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Prism"),
213  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_NodalPrism"),
216  BwdTrans_IterPerExp::create, "BwdTrans_IterPerExp_Hex"),
217 };
218 
219 
220 /**
221  * @brief Backward transform operator using LocalRegions implementation.
222  */
224 {
225  public:
227 
229  {
230  }
231 
232  virtual void operator()(
233  const Array<OneD, const NekDouble> &input,
234  Array<OneD, NekDouble> &output,
235  Array<OneD, NekDouble> &output1,
236  Array<OneD, NekDouble> &output2,
238  {
239  const int nCoeffs = m_expList[0]->GetNcoeffs();
240  const int nPhys = m_expList[0]->GetTotPoints();
242 
243  for (int i = 0; i < m_numElmt; ++i)
244  {
245  m_expList[i]->BwdTrans(input + i*nCoeffs,
246  tmp = output + i*nPhys);
247  }
248  }
249 
250  virtual void operator()(
251  int dir,
252  const Array<OneD, const NekDouble> &input,
253  Array<OneD, NekDouble> &output,
255  {
256  ASSERTL0(false, "Not valid for this operator.");
257  }
258 
259  protected:
260  vector<StdRegions::StdExpansionSharedPtr> m_expList;
261 
262  private:
264  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
265  CoalescedGeomDataSharedPtr pGeomData)
266  : Operator(pCollExp, pGeomData)
267  {
268  m_expList = pCollExp;
269  }
270 };
271 
272 /// Factory initialisation for the BwdTrans_NoCollection operators
273 OperatorKey BwdTrans_NoCollection::m_typeArr[] = {
276  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Seg"),
279  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Tri"),
282  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_NodalTri"),
285  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Quad"),
288  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Tet"),
291  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_NodalTet"),
294  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Pyr"),
297  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Prism"),
300  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_NodalPrism"),
303  BwdTrans_NoCollection::create, "BwdTrans_NoCollection_Hex"),
304 };
305 
306 
307 /**
308  * @brief Backward transform operator using sum-factorisation (Segment)
309  */
311 {
312  public:
314 
316  {
317  }
318 
319  virtual void operator()(
320  const Array<OneD, const NekDouble> &input,
321  Array<OneD, NekDouble> &output,
322  Array<OneD, NekDouble> &output1,
323  Array<OneD, NekDouble> &output2,
325  {
326  if(m_colldir0)
327  {
328  Vmath::Vcopy(m_numElmt*m_nmodes0,input.get(),1,output.get(),1);
329  }
330  else
331  {
332  // out = B0*in;
333  Blas::Dgemm('N','N', m_nquad0, m_numElmt, m_nmodes0,
334  1.0, m_base0.get(), m_nquad0,
335  &input[0], m_nmodes0, 0.0,
336  &output[0], m_nquad0);
337  }
338  }
339 
340  virtual void operator()(
341  int dir,
342  const Array<OneD, const NekDouble> &input,
343  Array<OneD, NekDouble> &output,
345  {
346  ASSERTL0(false, "Not valid for this operator.");
347  }
348 
349 
350  protected:
351  const int m_nquad0;
352  const int m_nmodes0;
353  const bool m_colldir0;
355 
356  private:
358  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
359  CoalescedGeomDataSharedPtr pGeomData)
360  : Operator (pCollExp, pGeomData),
361  m_nquad0 (m_stdExp->GetNumPoints(0)),
362  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
363  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
364  m_base0 (m_stdExp->GetBasis(0)->GetBdata())
365  {
366  m_wspSize = 0;
367  }
368 };
369 
370 /// Factory initialisation for the BwdTrans_SumFac_Seg operator
371 OperatorKey BwdTrans_SumFac_Seg::m_type = GetOperatorFactory().
372  RegisterCreatorFunction(
374  BwdTrans_SumFac_Seg::create, "BwdTrans_SumFac_Seg");
375 
376 
377 
378 /**
379  * @brief Backward transform operator using sum-factorisation (Quad)
380  */
382 {
383  public:
385 
387  {
388  }
389 
390  virtual void operator()(
391  const Array<OneD, const NekDouble> &input,
392  Array<OneD, NekDouble> &output,
393  Array<OneD, NekDouble> &output1,
394  Array<OneD, NekDouble> &output2,
396  {
397  int i = 0;
398  if(m_colldir0 && m_colldir1)
399  {
400  Vmath::Vcopy(m_numElmt*m_nmodes0*m_nmodes1,
401  input.get(), 1, output.get(), 1);
402  }
403  else if(m_colldir0)
404  {
406  = m_stdExp->GetBasis(1)->GetBdata();
407 
408  for(i = 0; i < m_numElmt; ++i)
409  {
410  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nmodes1,
411  1.0, &input[i*m_nquad0*m_nmodes1], m_nquad0,
412  base1.get(), m_nquad1, 0.0,
413  &output[i*m_nquad0*m_nquad1], m_nquad0);
414  }
415  }
416  else if(m_colldir1)
417  {
418  Blas::Dgemm('N', 'N', m_nquad0, m_nmodes1*m_numElmt, m_nmodes0,
419  1.0, m_base0.get(), m_nquad0,
420  &input[0], m_nmodes0, 0.0,
421  &output[0], m_nquad0);
422  }
423  else
424  {
425  ASSERTL1(wsp.num_elements() == m_wspSize,
426  "Incorrect workspace size");
427 
428  // Those two calls correpsond to the operation
429  // out = B0*in*Transpose(B1);
430  Blas::Dgemm('N', 'N', m_nquad0, m_nmodes1*m_numElmt, m_nmodes0,
431  1.0, m_base0.get(), m_nquad0,
432  &input[0], m_nmodes0, 0.0,
433  &wsp[0], m_nquad0);
434 
435  for(i = 0; i < m_numElmt; ++i)
436  {
437  Blas::Dgemm('N','T', m_nquad0, m_nquad1, m_nmodes1,
438  1.0, &wsp[i*m_nquad0*m_nmodes1], m_nquad0,
439  m_base1.get(), m_nquad1, 0.0,
440  &output[i*m_nquad0*m_nquad1], m_nquad0);
441  }
442  }
443  }
444 
445  virtual void operator()(
446  int dir,
447  const Array<OneD, const NekDouble> &input,
448  Array<OneD, NekDouble> &output,
450  {
451  ASSERTL0(false, "Not valid for this operator.");
452  }
453 
454  protected:
455  const int m_nquad0;
456  const int m_nquad1;
457  const int m_nmodes0;
458  const int m_nmodes1;
459  const bool m_colldir0;
460  const bool m_colldir1;
463 
464  private:
466  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
467  CoalescedGeomDataSharedPtr pGeomData)
468  : Operator (pCollExp, pGeomData),
469  m_nquad0 (m_stdExp->GetNumPoints(0)),
470  m_nquad1 (m_stdExp->GetNumPoints(1)),
471  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
472  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
473  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
474  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
475  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
476  m_base1 (m_stdExp->GetBasis(1)->GetBdata())
477  {
478  m_wspSize = m_nquad0*m_nmodes1*m_numElmt;
479  }
480 };
481 
482 /// Factory initialisation for the BwdTrans_SumFac_Quad operator
483 OperatorKey BwdTrans_SumFac_Quad::m_type = GetOperatorFactory().
484  RegisterCreatorFunction(
486  BwdTrans_SumFac_Quad::create, "BwdTrans_SumFac_Quad");
487 
488 
489 /**
490  * @brief Backward transform operator using sum-factorisation (Tri)
491  */
493 {
494  public:
496 
498  {
499  }
500 
501  virtual void operator()(
502  const Array<OneD, const NekDouble> &input,
503  Array<OneD, NekDouble> &output,
504  Array<OneD, NekDouble> &output1,
505  Array<OneD, NekDouble> &output2,
507  {
508 
509  ASSERTL1(wsp.num_elements() == m_wspSize,
510  "Incorrect workspace size");
511 
512  int ncoeffs = m_stdExp->GetNcoeffs();
513  int i = 0;
514  int mode = 0;
515 
516  for (i = mode = 0; i < m_nmodes0; ++i)
517  {
518  Blas::Dgemm('N', 'N', m_nquad1, m_numElmt, m_nmodes1-i,
519  1.0, m_base1.get()+mode*m_nquad1, m_nquad1,
520  &input[0]+mode, ncoeffs, 0.0,
521  &wsp[i*m_nquad1*m_numElmt], m_nquad1);
522  mode += m_nmodes1-i;
523  }
524 
525  // fix for modified basis by splitting top vertex mode
526  if(m_sortTopVertex)
527  {
528  for(i = 0; i < m_numElmt; ++i)
529  {
530  Blas::Daxpy(m_nquad1, input[1+i*ncoeffs],
531  m_base1.get() + m_nquad1, 1,
532  &wsp[m_nquad1*m_numElmt] + i*m_nquad1, 1);
533  }
534 
535  }
536 
537  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1*m_numElmt, m_nmodes0,
538  1.0, m_base0.get(), m_nquad0,
539  &wsp[0], m_nquad1*m_numElmt, 0.0,
540  &output[0], m_nquad0);
541  }
542 
543  virtual void operator()(
544  int dir,
545  const Array<OneD, const NekDouble> &input,
546  Array<OneD, NekDouble> &output,
548  {
549  ASSERTL0(false, "Not valid for this operator.");
550  }
551 
552 
553  protected:
554  const int m_nquad0;
555  const int m_nquad1;
556  const int m_nmodes0;
557  const int m_nmodes1;
561 
562  private:
564  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
565  CoalescedGeomDataSharedPtr pGeomData)
566  : Operator (pCollExp, pGeomData),
567  m_nquad0 (m_stdExp->GetNumPoints(0)),
568  m_nquad1 (m_stdExp->GetNumPoints(1)),
569  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
570  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
571  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
572  m_base1 (m_stdExp->GetBasis(1)->GetBdata())
573  {
574  m_wspSize = m_nquad0 * m_nmodes1 * m_numElmt;
575  if(m_stdExp->GetBasis(0)->GetBasisType()
577  {
578  m_sortTopVertex = true;
579  }
580  else
581  {
582  m_sortTopVertex = false;
583  }
584  }
585 
586 };
587 
588 /// Factory initialisation for the BwdTrans_SumFac_Tri operator
589 OperatorKey BwdTrans_SumFac_Tri::m_type = GetOperatorFactory().
590  RegisterCreatorFunction(
592  BwdTrans_SumFac_Tri::create, "BwdTrans_SumFac_Tri");
593 
594 
595 /// Backward transform operator using sum-factorisation (Hex)
597 {
598  public:
600 
602  {
603  }
604 
605  virtual void operator()(
606  const Array<OneD, const NekDouble> &input,
607  Array<OneD, NekDouble> &output,
608  Array<OneD, NekDouble> &output1,
609  Array<OneD, NekDouble> &output2,
611  {
612 
613  if(m_colldir0 && m_colldir1 && m_colldir2)
614  {
615  Vmath::Vcopy(m_numElmt * m_nmodes0 * m_nmodes1 * m_nmodes2,
616  input.get(), 1,
617  output.get(), 1);
618  }
619  else
620  {
621  ASSERTL1(wsp.num_elements() == m_wspSize,
622  "Incorrect workspace size");
623 
624  // Assign second half of workspace for 2nd DGEMM operation.
625  int totmodes = m_nmodes0*m_nmodes1*m_nmodes2;
626 
628  = wsp + m_nmodes0*m_nmodes1*m_nquad2*m_numElmt;
629 
630  //loop over elements and do bwd trans wrt c
631  for(int n = 0; n < m_numElmt; ++n)
632  {
633  Blas::Dgemm('N', 'T', m_nquad2, m_nmodes0*m_nmodes1,
634  m_nmodes2, 1.0, m_base2.get(), m_nquad2,
635  &input[n*totmodes], m_nmodes0*m_nmodes1, 0.0,
636  &wsp[n*m_nquad2], m_nquad2*m_numElmt);
637  }
638 
639  // trans wrt b
640  Blas::Dgemm('N', 'T', m_nquad1, m_nquad2*m_numElmt*m_nmodes0,
641  m_nmodes1, 1.0, m_base1.get(), m_nquad1,
642  wsp.get(), m_nquad2*m_numElmt*m_nmodes0, 0.0,
643  wsp2.get(), m_nquad1);
644 
645  // trans wrt a
646  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1*m_nquad2*m_numElmt,
647  m_nmodes0, 1.0, m_base0.get(), m_nquad0,
648  wsp2.get(), m_nquad1*m_nquad2*m_numElmt, 0.0,
649  output.get(), m_nquad0);
650  }
651  }
652 
653  virtual void operator()(
654  int dir,
655  const Array<OneD, const NekDouble> &input,
656  Array<OneD, NekDouble> &output,
658  {
659  ASSERTL0(false, "Not valid for this operator.");
660  }
661 
662  protected:
663  const int m_nquad0;
664  const int m_nquad1;
665  const int m_nquad2;
666  const int m_nmodes0;
667  const int m_nmodes1;
668  const int m_nmodes2;
672  const bool m_colldir0;
673  const bool m_colldir1;
674  const bool m_colldir2;
675 
676  private:
678  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
679  CoalescedGeomDataSharedPtr pGeomData)
680  : Operator (pCollExp, pGeomData),
681  m_nquad0 (pCollExp[0]->GetNumPoints(0)),
682  m_nquad1 (pCollExp[0]->GetNumPoints(1)),
683  m_nquad2 (pCollExp[0]->GetNumPoints(2)),
684  m_nmodes0 (pCollExp[0]->GetBasisNumModes(0)),
685  m_nmodes1 (pCollExp[0]->GetBasisNumModes(1)),
686  m_nmodes2 (pCollExp[0]->GetBasisNumModes(2)),
687  m_base0 (pCollExp[0]->GetBasis(0)->GetBdata()),
688  m_base1 (pCollExp[0]->GetBasis(1)->GetBdata()),
689  m_base2 (pCollExp[0]->GetBasis(2)->GetBdata()),
690  m_colldir0(pCollExp[0]->GetBasis(0)->Collocation()),
691  m_colldir1(pCollExp[0]->GetBasis(1)->Collocation()),
692  m_colldir2(pCollExp[0]->GetBasis(2)->Collocation())
693  {
694  m_wspSize = m_numElmt*m_nmodes0*(m_nmodes1*m_nquad2 +
695  m_nquad1*m_nquad2);
696  }
697 };
698 
699 /// Factory initialisation for the BwdTrans_SumFac_Hex operator
700 OperatorKey BwdTrans_SumFac_Hex::m_type = GetOperatorFactory().
701  RegisterCreatorFunction(
703  BwdTrans_SumFac_Hex::create, "BwdTrans_SumFac_Hex");
704 
705 
706 /**
707  * @brief Backward transform operator using sum-factorisation (Tet)
708  */
710 {
711  public:
713 
715  {
716  }
717 
718  virtual void operator()(
719  const Array<OneD, const NekDouble> &input,
720  Array<OneD, NekDouble> &output,
721  Array<OneD, NekDouble> &output1,
722  Array<OneD, NekDouble> &output2,
724  {
725  ASSERTL1(wsp.num_elements() == m_wspSize,
726  "Incorrect workspace size");
727 
728  Array<OneD, NekDouble > tmp = wsp;
729  Array<OneD, NekDouble > tmp1 = tmp
730  + m_numElmt*m_nquad2*m_nmodes0*(2*m_nmodes1-m_nmodes0+1)/2;
731 
732  int i = 0;
733  int j = 0;
734  int mode = 0;
735  int mode1 = 0;
736  int cnt = 0;
737  int ncoeffs = m_stdExp->GetNcoeffs();
738 
739  // Perform summation over '2' direction
740  for(i = 0; i < m_nmodes0; ++i)
741  {
742  for(j = 0; j < m_nmodes1-i; ++j, ++cnt)
743  {
744  Blas::Dgemm('N', 'N', m_nquad2, m_numElmt, m_nmodes2-i-j,
745  1.0, m_base2.get()+mode*m_nquad2, m_nquad2,
746  input.get()+ mode1, ncoeffs, 0.0,
747  tmp.get() + cnt*m_nquad2*m_numElmt, m_nquad2);
748  mode += m_nmodes2-i-j;
749  mode1 += m_nmodes2-i-j;
750  }
751 
752  //increment mode in case m_nmodes1!=m_nmodes2
753  mode += (m_nmodes2-m_nmodes1)*(m_nmodes2-m_nmodes1+1)/2;
754  }
755 
756  // vertex mode - currently (1+c)/2 x (1-b)/2 x (1-a)/2
757  // component is evaluated
758  if(m_sortTopEdge)
759  {
760  for(i = 0; i < m_numElmt; ++i)
761  {
762  // top singular vertex
763  // (1+c)/2 x (1+b)/2 x (1-a)/2 component
764  Blas::Daxpy(m_nquad2, input[1+i*ncoeffs],
765  m_base2.get() + m_nquad2, 1,
766  &tmp[m_nquad2*m_numElmt] + i*m_nquad2, 1);
767 
768  // top singular vertex
769  // (1+c)/2 x (1-b)/2 x (1+a)/2 component
770  Blas::Daxpy(m_nquad2, input[1+i*ncoeffs],
771  m_base2.get() + m_nquad2, 1,
772  &tmp[m_nmodes1*m_nquad2*m_numElmt]
773  + i*m_nquad2, 1);
774  }
775  }
776 
777  // Perform summation over '1' direction
778  mode = 0;
779  for(i = 0; i < m_nmodes0; ++i)
780  {
781  Blas::Dgemm('N', 'T', m_nquad1, m_nquad2*m_numElmt, m_nmodes1-i,
782  1.0, m_base1.get()+mode*m_nquad1, m_nquad1,
783  tmp.get() + mode*m_nquad2*m_numElmt,
784  m_nquad2*m_numElmt,
785  0.0, tmp1.get() + i*m_nquad1*m_nquad2*m_numElmt,
786  m_nquad1);
787  mode += m_nmodes1-i;
788  }
789 
790  // fix for modified basis by adding additional split of
791  // top and base singular vertex modes as well as singular
792  // edge
793  if(m_sortTopEdge)
794  {
795  // this could probably be a dgemv or higher if we
796  // made a specialised m_base1[m_nuqad1] array
797  // containing multiply copies
798  for(i = 0; i < m_numElmt; ++i)
799  {
800  // sort out singular vertices and singular
801  // edge components with (1+b)/2 (1+a)/2 form
802  for(int j = 0; j < m_nquad2; ++j)
803  {
804  Blas::Daxpy(m_nquad1,
805  tmp[m_nquad2*m_numElmt + i*m_nquad2 + j],
806  m_base1.get() + m_nquad1,1,
807  &tmp1[m_nquad1*m_nquad2*m_numElmt]
808  + i*m_nquad1*m_nquad2 + j*m_nquad1, 1);
809  }
810  }
811  }
812 
813  // Perform summation over '0' direction
814  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1*m_nquad2*m_numElmt,
815  m_nmodes0, 1.0, m_base0.get(), m_nquad0,
816  tmp1.get(), m_nquad1*m_nquad2*m_numElmt, 0.0,
817  output.get(), m_nquad0);
818 
819  }
820 
821  virtual void operator()(
822  int dir,
823  const Array<OneD, const NekDouble> &input,
824  Array<OneD, NekDouble> &output,
826  {
827  ASSERTL0(false, "Not valid for this operator.");
828  }
829 
830  protected:
831  const int m_nquad0;
832  const int m_nquad1;
833  const int m_nquad2;
834  const int m_nmodes0;
835  const int m_nmodes1;
836  const int m_nmodes2;
841 
842  private:
844  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
845  CoalescedGeomDataSharedPtr pGeomData)
846  : Operator (pCollExp, pGeomData),
847  m_nquad0 (m_stdExp->GetNumPoints(0)),
848  m_nquad1 (m_stdExp->GetNumPoints(1)),
849  m_nquad2 (m_stdExp->GetNumPoints(2)),
850  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
851  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
852  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
853  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
854  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
855  m_base2 (m_stdExp->GetBasis(2)->GetBdata())
856  {
857  m_wspSize = m_numElmt
858  * (m_nquad2*m_nmodes0*(2*m_nmodes1-m_nmodes0+1)/2
859  + m_nquad2*m_nquad1*m_nmodes0);
860 
861  if(m_stdExp->GetBasis(0)->GetBasisType()
863  {
864  m_sortTopEdge = true;
865  }
866  else
867  {
868  m_sortTopEdge = false;
869  }
870  }
871 };
872 
873 /// Factory initialisation for the BwdTrans_SumFac_Tet operator
874 OperatorKey BwdTrans_SumFac_Tet::m_type = GetOperatorFactory().
875  RegisterCreatorFunction(
877  BwdTrans_SumFac_Tet::create, "BwdTrans_SumFac_Tet");
878 
879 
880 /**
881  * @brief Backward transform operator using sum-factorisation (Prism)
882  */
884 {
885  public:
887 
889  {
890  }
891 
892  virtual void operator()(
893  const Array<OneD, const NekDouble> &input,
894  Array<OneD, NekDouble> &output,
895  Array<OneD, NekDouble> &output1,
896  Array<OneD, NekDouble> &output2,
898  {
899  ASSERTL1(wsp.num_elements() == m_wspSize,
900  "Incorrect workspace size");
901 
902  // Assign second half of workspace for 2nd DGEMM operation.
903  int totmodes = m_stdExp->GetNcoeffs();
904 
906  = wsp + m_nmodes0*m_nmodes1*m_nquad2*m_numElmt;
907 
908  Vmath::Zero(m_nmodes0*m_nmodes1*m_nquad2*m_numElmt, wsp, 1);
909  int i = 0;
910  int j = 0;
911  int mode = 0;
912  int mode1 = 0;
913  int cnt = 0;
914  for (i = mode = mode1 = 0; i < m_nmodes0; ++i)
915  {
916  cnt = i * m_nquad2 * m_numElmt;
917  for (j = 0; j < m_nmodes1; ++j)
918  {
919  Blas::Dgemm('N', 'N', m_nquad2, m_numElmt, m_nmodes2-i,
920  1.0, m_base2.get()+mode*m_nquad2, m_nquad2,
921  &input[0]+mode1, totmodes, 0.0,
922  &wsp[j*m_nquad2*m_numElmt*m_nmodes0+ cnt],
923  m_nquad2);
924  mode1 += m_nmodes2-i;
925  }
926  mode += m_nmodes2-i;
927  }
928 
929  // fix for modified basis by splitting top vertex mode
930  if(m_sortTopVertex)
931  {
932  for(j = 0; j < m_nmodes1; ++j)
933  {
934  for(i = 0; i < m_numElmt; ++i)
935  {
936  Blas::Daxpy(m_nquad2,input[1+i*totmodes+j*m_nmodes2],
937  m_base2.get()+m_nquad2,1,
938  &wsp[j*m_nquad2*m_numElmt*m_nmodes0 +
939  m_nquad2*m_numElmt]+i*m_nquad2,1);
940  }
941  }
942  // Believe this could be made into a m_nmodes1
943  // dgemv if we made an array of m_numElmt copies
944  // of m_base2[m_quad2] (which are of size
945  // m_nquad2.
946  }
947 
948  // Perform summation over '1' direction
949  Blas::Dgemm('N', 'T', m_nquad1, m_nquad2*m_numElmt*m_nmodes0,
950  m_nmodes1,1.0, m_base1.get(), m_nquad1,
951  wsp.get(), m_nquad2*m_numElmt*m_nmodes0,
952  0.0, wsp2.get(), m_nquad1);
953 
954 
955  // Perform summation over '0' direction
956  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1*m_nquad2*m_numElmt,
957  m_nmodes0, 1.0, m_base0.get(), m_nquad0,
958  wsp2.get(), m_nquad1*m_nquad2*m_numElmt,
959  0.0, output.get(), m_nquad0);
960  }
961 
962  virtual void operator()(
963  int dir,
964  const Array<OneD, const NekDouble> &input,
965  Array<OneD, NekDouble> &output,
967  {
968  ASSERTL0(false, "Not valid for this operator.");
969  }
970 
971 
972  protected:
973  const int m_nquad0;
974  const int m_nquad1;
975  const int m_nquad2;
976  const int m_nmodes0;
977  const int m_nmodes1;
978  const int m_nmodes2;
983 
984  private:
986  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
987  CoalescedGeomDataSharedPtr pGeomData)
988  : Operator (pCollExp, pGeomData),
989  m_nquad0 (m_stdExp->GetNumPoints(0)),
990  m_nquad1 (m_stdExp->GetNumPoints(1)),
991  m_nquad2 (m_stdExp->GetNumPoints(2)),
992  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
993  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
994  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
995  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
996  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
997  m_base2 (m_stdExp->GetBasis(2)->GetBdata())
998  {
999  m_wspSize = m_numElmt*m_nmodes0*(m_nmodes1*m_nquad2
1000  + m_nquad1*m_nquad2);
1001 
1002  if(m_stdExp->GetBasis(0)->GetBasisType()
1004  {
1005  m_sortTopVertex = true;
1006  }
1007  else
1008  {
1009  m_sortTopVertex = false;
1010  }
1011 
1012  }
1013 };
1014 
1015 /// Factory initialisation for the BwdTrans_SumFac_Prism operator
1016 OperatorKey BwdTrans_SumFac_Prism::m_type = GetOperatorFactory().
1017  RegisterCreatorFunction(
1019  BwdTrans_SumFac_Prism::create, "BwdTrans_SumFac_Prism");
1020 
1021 }
1022 }
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:892
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:150
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Array< OneD, const NekDouble > m_base1
Definition: BwdTrans.cpp:980
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:501
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:78
boost::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:159
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:979
Backward transform operator using sum-factorisation (Tet)
Definition: BwdTrans.cpp:709
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:250
Backward transform operator using sum-factorisation (Hex)
Definition: BwdTrans.cpp:596
Principle Modified Functions .
Definition: BasisType.h:49
STL namespace.
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:543
Array< OneD, const NekDouble > m_base1
Definition: BwdTrans.cpp:670
Backward transform operator using sum-factorisation (Tri)
Definition: BwdTrans.cpp:492
Base class for operators on a collection of elements.
Definition: Operator.h:108
BwdTrans_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:563
boost::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:70
Backward transform operator using sum-factorisation (Prism)
Definition: BwdTrans.cpp:883
Backward transform operator using sum-factorisation (Segment)
Definition: BwdTrans.cpp:310
BwdTrans_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:357
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:319
Backward transform operator using LocalRegions implementation.
Definition: BwdTrans.cpp:223
Array< OneD, const NekDouble > m_base2
Definition: BwdTrans.cpp:981
BwdTrans_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:263
vector< StdRegions::StdExpansionSharedPtr > m_expList
Definition: BwdTrans.cpp:260
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:962
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:558
BwdTrans_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:465
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:354
BwdTrans_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:677
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:110
Array< OneD, const NekDouble > m_base1
Definition: BwdTrans.cpp:462
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:340
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:837
Backward transform operator using default StdRegions operator.
Definition: BwdTrans.cpp:141
Array< OneD, const NekDouble > m_base2
Definition: BwdTrans.cpp:671
#define OPERATOR_CREATE(cname)
Definition: Operator.h:44
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:669
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:232
Backward transform operator using sum-factorisation (Quad)
Definition: BwdTrans.cpp:381
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:653
Array< OneD, const NekDouble > m_base1
Definition: BwdTrans.cpp:559
BwdTrans_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:843
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:167
Array< OneD, const NekDouble > m_base2
Definition: BwdTrans.cpp:839
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:390
Array< OneD, const NekDouble > m_base1
Definition: BwdTrans.cpp:838
boost::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:65
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:605
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
#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
BwdTrans_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:985
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:821
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: BwdTrans.cpp:445
Array< OneD, const NekDouble > m_base0
Definition: BwdTrans.cpp:461
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: BwdTrans.cpp:718
Backward transform operator using standard matrix approach.
Definition: BwdTrans.cpp:56
BwdTrans_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:177
BwdTrans_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: BwdTrans.cpp:91
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215