Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
IProductWRTDerivBase.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: IProductWRTDerivBase.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: IProductWRTDerivBase operator implementations
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <loki/Singleton.h>
37 #include <Collections/Operator.h>
38 #include <Collections/Collection.h>
39 #include <Collections/IProduct.h>
40 
41 using namespace std;
42 
43 namespace Nektar {
44 namespace Collections {
45 
53 
54 /**
55  * @brief Inner product WRT deriv base operator using standard matrix approach
56  */
58 {
59  public:
61 
63  {
64  }
65 
66  virtual void operator()(
67  const Array<OneD, const NekDouble> &entry0,
68  Array<OneD, NekDouble> &entry1,
69  Array<OneD, NekDouble> &entry2,
70  Array<OneD, NekDouble> &entry3,
72  {
73  int nPhys = m_stdExp->GetTotPoints();
74  int ntot = m_numElmt*nPhys;
75  int nmodes = m_stdExp->GetNcoeffs();
79 
80  in[0] = entry0; in[1] = entry1;
81  in[2] = entry2;
82 
83  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
84  entry2: entry1;
85 
86  for(int i = 0; i < m_dim; ++i)
87  {
88  tmp[i] = wsp + i*ntot;
89  }
90 
91  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
92  for(int i = 0; i < m_dim; ++i)
93  {
94  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
95  tmp[i],1);
96  for(int j = 1; j < m_coordim; ++j)
97  {
98  Vmath::Vvtvp (ntot,m_derivFac[i +j*m_dim],1,
99  in[j],1, tmp[i], 1, tmp[i],1);
100  }
101  }
102 
103  // calculate Iproduct WRT Std Deriv
104 
105  // First component
106  Vmath::Vmul(ntot,m_jac,1,tmp[0],1,tmp[0],1);
107  Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[0]->GetRows(),
108  m_numElmt,m_iProdWRTStdDBase[0]->GetColumns(),
109  1.0, m_iProdWRTStdDBase[0]->GetRawPtr(),
110  m_iProdWRTStdDBase[0]->GetRows(),
111  tmp[0].get(), nPhys, 0.0,
112  output.get(), nmodes);
113 
114  // Other components
115  for(int i = 1; i < m_dim; ++i)
116  {
117  Vmath::Vmul(ntot,m_jac,1,tmp[i],1,tmp[i],1);
118  Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[i]->GetRows(),
119  m_numElmt,m_iProdWRTStdDBase[i]->GetColumns(),
120  1.0, m_iProdWRTStdDBase[i]->GetRawPtr(),
121  m_iProdWRTStdDBase[i]->GetRows(),
122  tmp[i].get(), nPhys, 1.0,
123  output.get(), nmodes);
124  }
125  }
126 
127  virtual void operator()(
128  int dir,
129  const Array<OneD, const NekDouble> &input,
130  Array<OneD, NekDouble> &output,
132  {
133  ASSERTL0(false, "Not valid for this operator.");
134  }
135 
136  protected:
140  int m_dim;
142 
143  private:
145  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
146  CoalescedGeomDataSharedPtr pGeomData)
147  : Operator(pCollExp, pGeomData)
148  {
149  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
150  m_dim = PtsKey.size();
151  m_coordim = m_stdExp->GetCoordim();
152 
153  int nqtot = m_stdExp->GetTotPoints();
154  int nmodes = m_stdExp->GetNcoeffs();
155 
156  // set up a IProductWRTDerivBase StdMat.
157  m_iProdWRTStdDBase = Array<OneD, DNekMatSharedPtr>(m_dim);
158  for(int i = 0; i < m_dim; ++i)
159  {
160  Array<OneD, NekDouble> tmp(nqtot),tmp1(nmodes);
161  m_iProdWRTStdDBase[i] = MemoryManager<DNekMat>
162  ::AllocateSharedPtr(nmodes,nqtot);
163  for(int j = 0; j < nqtot; ++j)
164  {
165  Vmath::Zero(nqtot,tmp,1);
166  tmp[j] = 1.0;
167  m_stdExp->IProductWRTDerivBase(i,tmp,tmp1);
168  Vmath::Vcopy(nmodes, &tmp1[0],1,
169  &(m_iProdWRTStdDBase[i]->GetPtr())[0]
170  + j*nmodes, 1);
171  }
172  }
173  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
174  m_jac = pGeomData->GetJac(pCollExp);
175  m_wspSize = m_dim*nqtot*m_numElmt;
176  }
177 };
178 
179 /// Factory initialisation for the IProductWRTDerivBase_StdMat operators
180 OperatorKey IProductWRTDerivBase_StdMat::m_typeArr[] = {
183  IProductWRTDerivBase_StdMat::create,
184  "IProductWRTDerivBase_StdMat_Seg"),
187  IProductWRTDerivBase_StdMat::create,
188  "IProductWRTDerivBase_StdMat_Tri"),
191  IProductWRTDerivBase_StdMat::create,
192  "IProductWRTDerivBase_StdMat_NodalTri"),
195  IProductWRTDerivBase_StdMat::create,
196  "IProductWRTDerivBase_StdMat_Quad"),
199  IProductWRTDerivBase_StdMat::create,
200  "IProductWRTDerivBase_StdMat_Tet"),
203  IProductWRTDerivBase_StdMat::create,
204  "IProductWRTDerivBase_StdMat_NodalTet"),
207  IProductWRTDerivBase_StdMat::create,
208  "IProductWRTDerivBase_StdMat_Pyr"),
211  IProductWRTDerivBase_StdMat::create,
212  "IProductWRTDerivBase_StdMat_Prism"),
215  IProductWRTDerivBase_StdMat::create,
216  "IProductWRTDerivBase_StdMat_NodalPrism"),
219  IProductWRTDerivBase_StdMat::create,
220  "IProductWRTDerivBase_StdMat_Hex")
221 };
222 
223 
224 /**
225  * @brief Inner product WRT deriv base operator using element-wise operation
226  */
228 {
229  public:
231 
233  {
234  }
235 
236  virtual void operator()(
237  const Array<OneD, const NekDouble> &entry0,
238  Array<OneD, NekDouble> &entry1,
239  Array<OneD, NekDouble> &entry2,
240  Array<OneD, NekDouble> &entry3,
242  {
243  unsigned int nPhys = m_stdExp->GetTotPoints();
244  unsigned int ntot = m_numElmt*nPhys;
245  unsigned int nmodes = m_stdExp->GetNcoeffs();
246  unsigned int nmax = max(ntot,m_numElmt*nmodes);
248  Array<OneD, NekDouble> output, tmp1;
250 
251  in[0] = entry0; in[1] = entry1; in[2] = entry2;
252 
253  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
254  entry2: entry1;
255 
256  for(int i = 0; i < m_dim; ++i)
257  {
258  tmp[i] = wsp + i*nmax;
259  }
260 
261  // calculate dx/dxi in[0] + dy/dxi in[2] + dz/dxi in[3]
262  for(int i = 0; i < m_dim; ++i)
263  {
264  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
265  tmp[i],1);
266  for(int j = 1; j < m_coordim; ++j)
267  {
268  Vmath::Vvtvp (ntot,m_derivFac[i +j*m_dim],1,
269  in[j],1, tmp[i], 1, tmp[i],1);
270  }
271  }
272 
273  // calculate Iproduct WRT Std Deriv
274  // first component
275  Vmath::Vmul(ntot,m_jac,1,tmp[0],1,tmp[0],1);
276  for(int n = 0; n < m_numElmt; ++n)
277  {
278  m_stdExp->IProductWRTDerivBase(0,tmp[0]+n*nPhys,
279  tmp1 = output + n*nmodes);
280  }
281 
282  // other components
283  for(int i = 1; i < m_dim; ++i)
284  {
285  // multiply by Jacobian
286  Vmath::Vmul(ntot,m_jac,1,tmp[i],1,tmp[i],1);
287  for(int n = 0; n < m_numElmt; ++n)
288  {
289  m_stdExp->IProductWRTDerivBase(i,tmp[i]+n*nPhys,tmp[0]);
290  Vmath::Vadd(nmodes,tmp[0],1,output+n*nmodes,1,
291  tmp1 = output+n*nmodes,1);
292  }
293  }
294  }
295 
296  virtual void operator()(
297  int dir,
298  const Array<OneD, const NekDouble> &input,
299  Array<OneD, NekDouble> &output,
301  {
302  ASSERTL0(false, "Not valid for this operator.");
303  }
304 
305  protected:
308  int m_dim;
310 
311  private:
313  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
314  CoalescedGeomDataSharedPtr pGeomData)
315  : Operator(pCollExp, pGeomData)
316  {
317  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
318  m_dim = PtsKey.size();
319  m_coordim = m_stdExp->GetCoordim();
320 
321  int nqtot = m_stdExp->GetTotPoints();
322 
323  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
324  m_jac = pGeomData->GetJac(pCollExp);
325  m_wspSize = m_dim*nqtot*m_numElmt;
326  }
327 };
328 
329 /// Factory initialisation for the IProductWRTDerivBase_IterPerExp operators
330 OperatorKey IProductWRTDerivBase_IterPerExp::m_typeArr[] = {
333  IProductWRTDerivBase_IterPerExp::create,
334  "IProductWRTDerivBase_IterPerExp_Seg"),
337  IProductWRTDerivBase_IterPerExp::create,
338  "IProductWRTDerivBase_IterPerExp_Tri"),
341  IProductWRTDerivBase_IterPerExp::create,
342  "IProductWRTDerivBase_IterPerExp_NodalTri"),
345  IProductWRTDerivBase_IterPerExp::create,
346  "IProductWRTDerivBase_IterPerExp_Quad"),
349  IProductWRTDerivBase_IterPerExp::create,
350  "IProductWRTDerivBase_IterPerExp_Tet"),
353  IProductWRTDerivBase_IterPerExp::create,
354  "IProductWRTDerivBase_IterPerExp_NodalTet"),
357  IProductWRTDerivBase_IterPerExp::create,
358  "IProductWRTDerivBase_IterPerExp_Pyr"),
361  IProductWRTDerivBase_IterPerExp::create,
362  "IProductWRTDerivBase_IterPerExp_Prism"),
365  IProductWRTDerivBase_IterPerExp::create,
366  "IProductWRTDerivBase_IterPerExp_NodalPrism"),
369  IProductWRTDerivBase_IterPerExp::create,
370  "IProductWRTDerivBase_IterPerExp_Hex")
371 };
372 
373 
374 /**
375  * @brief Inner product WRT deriv base operator using LocalRegions
376  * implementation.
377  */
379 {
380  public:
382 
384  {
385  }
386 
387  virtual void operator()(
388  const Array<OneD, const NekDouble> &entry0,
389  Array<OneD, NekDouble> &entry1,
390  Array<OneD, NekDouble> &entry2,
391  Array<OneD, NekDouble> &entry3,
393  {
394  unsigned int nmodes = m_expList[0]->GetNcoeffs();
395  unsigned int nPhys = m_expList[0]->GetTotPoints();
396  Array<OneD, NekDouble> tmp(nmodes),tmp1;
397 
399  Array<OneD, NekDouble> output;
400  in[0] = entry0; in[1] = entry1; in[2] = entry2;
401 
402  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
403  entry2: entry1;
404 
405  for(int n = 0; n < m_numElmt; ++n)
406  {
407  m_expList[n]->IProductWRTDerivBase(0,
408  in[0] + n * nPhys,
409  tmp1 = output + n * nmodes);
410  }
411 
412  for(int i = 1; i < m_dim; ++i)
413  {
414  for(int n = 0; n < m_numElmt; ++n)
415  {
416  m_expList[n]->IProductWRTDerivBase(i,in[i]+n*nPhys,tmp);
417 
418  Vmath::Vadd(nmodes,tmp,1,output+n*nmodes,1,
419  tmp1 = output+n*nmodes,1);
420  }
421  }
422  }
423 
424  virtual void operator()(
425  int dir,
426  const Array<OneD, const NekDouble> &input,
427  Array<OneD, NekDouble> &output,
429  {
430  ASSERTL0(false, "Not valid for this operator.");
431  }
432 
433  protected:
434  int m_dim;
436  vector<StdRegions::StdExpansionSharedPtr> m_expList;
437 
438  private:
440  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
441  CoalescedGeomDataSharedPtr pGeomData)
442  : Operator(pCollExp, pGeomData)
443  {
444  m_expList = pCollExp;
445  m_dim = pCollExp[0]->GetNumBases();
446  m_coordim = pCollExp[0]->GetCoordim();
447  }
448 };
449 
450 /// Factory initialisation for the IProductWRTDerivBase_NoCollection operators
451 OperatorKey IProductWRTDerivBase_NoCollection::m_typeArr[] = {
454  IProductWRTDerivBase_NoCollection::create,
455  "IProductWRTDerivBase_NoCollection_Seg"),
458  IProductWRTDerivBase_NoCollection::create,
459  "IProductWRTDerivBase_NoCollection_Tri"),
462  IProductWRTDerivBase_NoCollection::create,
463  "IProductWRTDerivBase_NoCollection_NodalTri"),
466  IProductWRTDerivBase_NoCollection::create,
467  "IProductWRTDerivBase_NoCollection_Quad"),
470  IProductWRTDerivBase_NoCollection::create,
471  "IProductWRTDerivBase_NoCollection_Tet"),
474  IProductWRTDerivBase_NoCollection::create,
475  "IProductWRTDerivBase_NoCollection_NodalTet"),
478  IProductWRTDerivBase_NoCollection::create,
479  "IProductWRTDerivBase_NoCollection_Pyr"),
482  IProductWRTDerivBase_NoCollection::create,
483  "IProductWRTDerivBase_NoCollection_Prism"),
486  IProductWRTDerivBase_NoCollection::create,
487  "IProductWRTDerivBase_NoCollection_NodalPrism"),
490  IProductWRTDerivBase_NoCollection::create,
491  "IProductWRTDerivBase_NoCollection_Hex")
492 };
493 
494 
495 /**
496  * @brief Inner product WRT deriv base operator using sum-factorisation
497  * (Segment)
498  */
500 {
501  public:
503 
505  {
506  }
507 
508  virtual void operator()(
509  const Array<OneD, const NekDouble> &input,
510  Array<OneD, NekDouble> &output,
511  Array<OneD, NekDouble> &output1,
512  Array<OneD, NekDouble> &output2,
514  {
515 
516  Vmath::Vmul(m_numElmt*m_nquad0, m_jac, 1, input, 1, wsp, 1);
517  Vmath::Vmul(m_numElmt*m_nquad0, &m_derivFac[0][0], 1,
518  &wsp[0], 1,
519  &wsp[0], 1);
520 
521  // out = B0*in;
522  Blas::Dgemm('T', 'N', m_nmodes0, m_numElmt, m_nquad0,
523  1.0, m_derbase0.get(), m_nquad0,
524  &wsp[0], m_nquad0, 0.0,
525  &output[0], m_nmodes0);
526  }
527 
528  virtual void operator()(
529  int dir,
530  const Array<OneD, const NekDouble> &input,
531  Array<OneD, NekDouble> &output,
533  {
534  ASSERTL0(false, "Not valid for this operator.");
535  }
536 
537  protected:
538  const int m_nquad0;
539  const int m_nmodes0;
543 
544  private:
546  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
547  CoalescedGeomDataSharedPtr pGeomData)
548  : Operator (pCollExp, pGeomData),
549  m_nquad0 (m_stdExp->GetNumPoints(0)),
550  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
551  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata())
552  {
553  m_wspSize = m_numElmt*m_nquad0;
554  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
555  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
556  }
557 };
558 
559 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Seg operator
560 OperatorKey IProductWRTDerivBase_SumFac_Seg::m_type = GetOperatorFactory().
561  RegisterCreatorFunction(
563  IProductWRTDerivBase_SumFac_Seg::create,
564  "IProductWRTDerivBase_SumFac_Seg");
565 
566 
567 /**
568  * @brief Inner product WRT deriv base operator using sum-factorisation (Quad)
569  */
571 {
572  public:
574 
576  {
577  }
578 
579  virtual void operator()(
580  const Array<OneD, const NekDouble> &entry0,
581  Array<OneD, NekDouble> &entry1,
582  Array<OneD, NekDouble> &entry2,
583  Array<OneD, NekDouble> &entry3,
585  {
586  unsigned int nPhys = m_stdExp->GetTotPoints();
587  unsigned int ntot = m_numElmt*nPhys;
588  unsigned int nmodes = m_stdExp->GetNcoeffs();
589  unsigned int nmax = max(ntot,m_numElmt*nmodes);
591  Array<OneD, NekDouble> output, wsp1;
593 
594  in[0] = entry0; in[1] = entry1; in[2] = entry2;
595 
596  output = (m_coordim == 2)? entry2: entry3;
597 
598  tmp[0] = wsp; tmp[1] = wsp + nmax;
599  wsp1 = wsp + 2*nmax;
600 
601  // calculate dx/dxi in[0] + dy/dxi in[1]
602  for(int i = 0; i < 2; ++i)
603  {
604  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
605  tmp[i],1);
606  for(int j = 1; j < m_coordim; ++j)
607  {
608  Vmath::Vvtvp (ntot,m_derivFac[i +j*2],1,
609  in[j],1, tmp[i], 1, tmp[i],1);
610  }
611  }
612 
613  // Iproduct wrt derivative of base 0
614  QuadIProduct(false, m_colldir1,m_numElmt,
615  m_nquad0, m_nquad1,
616  m_nmodes0, m_nmodes1,
617  m_derbase0, m_base1,
618  m_jac, tmp[0], output, wsp1);
619 
620  // Iproduct wrt derivative of base 1
621  QuadIProduct(m_colldir0, false, m_numElmt,
622  m_nquad0, m_nquad1,
623  m_nmodes0, m_nmodes1,
624  m_base0, m_derbase1,
625  m_jac, tmp[1], tmp[0], wsp1);
626 
627  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
628  }
629 
630  virtual void operator()(
631  int dir,
632  const Array<OneD, const NekDouble> &input,
633  Array<OneD, NekDouble> &output,
635  {
636  ASSERTL0(false, "Not valid for this operator.");
637  }
638 
639  protected:
640  const int m_nquad0;
641  const int m_nquad1;
642  const int m_nmodes0;
643  const int m_nmodes1;
644  const bool m_colldir0;
645  const bool m_colldir1;
653 
654  private:
656  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
657  CoalescedGeomDataSharedPtr pGeomData)
658  : Operator(pCollExp, pGeomData),
659  m_nquad0 (m_stdExp->GetNumPoints(0)),
660  m_nquad1 (m_stdExp->GetNumPoints(1)),
661  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
662  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
663  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
664  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
665  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
666  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
667  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
668  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
669  {
670  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
671  m_coordim = m_stdExp->GetCoordim();
672 
673  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
674  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
675  m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
676  m_nmodes0*m_nmodes1));
677  }
678 };
679 
680 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Quad operator
681 OperatorKey IProductWRTDerivBase_SumFac_Quad::m_type =
684  IProductWRTDerivBase_SumFac_Quad::create,
685  "IProductWRTDerivBase_IterPerExp_Quad");
686 
687 
688 /**
689  * @brief Inner product WRT deriv base operator using sum-factorisation (Tri)
690  */
692 {
693  public:
695 
697  {
698  }
699 
700  virtual void operator()(
701  const Array<OneD, const NekDouble> &entry0,
702  Array<OneD, NekDouble> &entry1,
703  Array<OneD, NekDouble> &entry2,
704  Array<OneD, NekDouble> &entry3,
706  {
707  unsigned int nPhys = m_stdExp->GetTotPoints();
708  unsigned int ntot = m_numElmt*nPhys;
709  unsigned int nmodes = m_stdExp->GetNcoeffs();
710  unsigned int nmax = max(ntot,m_numElmt*nmodes);
712  Array<OneD, NekDouble> output, wsp1;
714 
715  in[0] = entry0; in[1] = entry1; in[2] = entry2;
716 
717  output = (m_coordim == 2)? entry2: entry3;
718 
719  tmp[0] = wsp; tmp[1] = wsp + nmax;
720  wsp1 = wsp + 2*nmax;
721 
722 
723  // calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
724  // dphi/dxi_1 dxi_1/dx),in[0])
725  // + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
726  // dphi/dxi_1 dxi_1/dy),in[1])
727  //
728  // Note dphi/dxi_0 =
729  // dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_1)
730  //
731  // dphi/dxi_1 =
732  // dphi/deta_1 deta_1/dxi_1 + dphi/deta_1 deta_1/dxi_1 =
733  // dphi/deta_0 (1+eta_0)/(1-eta_1) + dphi/deta_1
734  //
735  // and so the full inner products are
736  //
737  // (dphi/dx,in[0]) + (dphi/dy,in[1])
738  // = (dphi/deta_0, ((2/(1-eta_1) (dxi_0/dx in[0]+dxi_0/dy in[1])
739  // + (1_eta_0)/(1-eta_1) (dxi_1/dx in[0]+dxi_1/dy in[1]))
740  // + (dphi/deta_1, (dxi_1/dx in[0] + dxi_1/dy in[1]))
741 
742  for(int i = 0; i < 2; ++i)
743  {
744  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1, tmp[i],1);
745 
746  for(int j = 1; j < m_coordim; ++j)
747  {
748  Vmath::Vvtvp (ntot,m_derivFac[i +j*2],1,
749  in[j],1, tmp[i], 1, tmp[i],1);
750  }
751  }
752 
753  // Multiply by factor: 2/(1-z1)
754  for (int i = 0; i < m_numElmt; ++i)
755  {
756  // scale tmp[0] by geometric factor: 2/(1-z1)
757  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
758  tmp[0].get()+i*nPhys,1);
759 
760  // scale tmp[1] by geometric factor (1+z0)/(1-z1)
761  Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[1].get()+i*nPhys,1,
762  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
763  }
764 
765  // Iproduct wrt derivative of base 0
766  TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
767  m_nmodes0, m_nmodes1, m_derbase0, m_base1,
768  m_jac, tmp[0], output, wsp1);
769 
770  // Iproduct wrt derivative of base 1
771  TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
772  m_nmodes0, m_nmodes1, m_base0, m_derbase1,
773  m_jac, tmp[1], tmp[0], wsp1);
774 
775  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
776  }
777 
778  virtual void operator()(
779  int dir,
780  const Array<OneD, const NekDouble> &input,
781  Array<OneD, NekDouble> &output,
783  {
784  ASSERTL0(false, "Not valid for this operator.");
785  }
786 
787  protected:
788  const int m_nquad0;
789  const int m_nquad1;
790  const int m_nmodes0;
791  const int m_nmodes1;
792  const bool m_colldir0;
793  const bool m_colldir1;
804 
805  private:
807  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
808  CoalescedGeomDataSharedPtr pGeomData)
809  : Operator(pCollExp, pGeomData),
810  m_nquad0 (m_stdExp->GetNumPoints(0)),
811  m_nquad1 (m_stdExp->GetNumPoints(1)),
812  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
813  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
814  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
815  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
816  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
817  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
818  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
819  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
820  {
821  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
822  m_coordim = m_stdExp->GetCoordim();
823 
824  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
825  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
826  m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
827  m_nmodes0*m_nmodes1));
828 
829  if(m_stdExp->GetBasis(0)->GetBasisType()
831  {
832  m_sortTopVertex = true;
833  }
834  else
835  {
836  m_sortTopVertex = false;
837  }
838 
840  = m_stdExp->GetBasis(0)->GetZ();
842  = m_stdExp->GetBasis(1)->GetZ();
843 
844  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
845  // set up geometric factor: 2/(1-z1)
846  for (int i = 0; i < m_nquad0; ++i)
847  {
848  for(int j = 0; j < m_nquad1; ++j)
849  {
850  m_fac0[i+j*m_nquad0] = 2.0/(1-z1[j]);
851  }
852  }
853 
854  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
855  // set up geometric factor: (1+z0)/(1-z1)
856  for (int i = 0; i < m_nquad0; ++i)
857  {
858  for(int j = 0; j < m_nquad1; ++j)
859  {
860  m_fac1[i+j*m_nquad0] = (1+z0[i])/(1-z1[j]);
861  }
862  }
863  }
864 };
865 
866 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Tri operator
867 OperatorKey IProductWRTDerivBase_SumFac_Tri::m_type =
870  IProductWRTDerivBase_SumFac_Tri::create,
871  "IProductWRTDerivBase_IterPerExp_Tri");
872 
873 
874 /**
875  * @brief Inner product WRT deriv base operator using sum-factorisation (Hex)
876  */
878 {
879  public:
881 
883  {
884  }
885 
886  virtual void operator()(
887  const Array<OneD, const NekDouble> &entry0,
888  Array<OneD, NekDouble> &entry1,
889  Array<OneD, NekDouble> &entry2,
890  Array<OneD, NekDouble> &entry3,
892  {
893  unsigned int nPhys = m_stdExp->GetTotPoints();
894  unsigned int ntot = m_numElmt*nPhys;
895  unsigned int nmodes = m_stdExp->GetNcoeffs();
896  unsigned int nmax = max(ntot,m_numElmt*nmodes);
898  Array<OneD, NekDouble> output, wsp1;
900 
901  in[0] = entry0; in[1] = entry1;
902  in[2] = entry2;
903 
904  output = entry3;
905 
906  for(int i = 0; i < 3; ++i)
907  {
908  tmp[i] = wsp + i*nmax;
909  }
910 
911  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
912  for(int i = 0; i < 3; ++i)
913  {
914  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
915  tmp[i],1);
916  for(int j = 1; j < 3; ++j)
917  {
918  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
919  in[j],1, tmp[i], 1, tmp[i],1);
920  }
921  }
922 
923  wsp1 = wsp + 3*nmax;
924 
925  // calculate Iproduct WRT Std Deriv
926  HexIProduct(false,m_colldir1,m_colldir2, m_numElmt,
927  m_nquad0, m_nquad1, m_nquad2,
928  m_nmodes0, m_nmodes1, m_nmodes2,
929  m_derbase0, m_base1, m_base2,
930  m_jac,tmp[0],output,wsp1);
931 
932  HexIProduct(m_colldir0,false,m_colldir2, m_numElmt,
933  m_nquad0, m_nquad1, m_nquad2,
934  m_nmodes0, m_nmodes1, m_nmodes2,
935  m_base0, m_derbase1, m_base2,
936  m_jac,tmp[1],tmp[0],wsp1);
937  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
938 
939  HexIProduct(m_colldir0,m_colldir1,false, m_numElmt,
940  m_nquad0, m_nquad1, m_nquad2,
941  m_nmodes0, m_nmodes1, m_nmodes2,
942  m_base0, m_base1, m_derbase2,
943  m_jac,tmp[2],tmp[0],wsp1);
944  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
945  }
946 
947  virtual void operator()(
948  int dir,
949  const Array<OneD, const NekDouble> &input,
950  Array<OneD, NekDouble> &output,
952  {
953  ASSERTL0(false, "Not valid for this operator.");
954  }
955 
956  protected:
957  const int m_nquad0;
958  const int m_nquad1;
959  const int m_nquad2;
960  const int m_nmodes0;
961  const int m_nmodes1;
962  const int m_nmodes2;
963  const bool m_colldir0;
964  const bool m_colldir1;
965  const bool m_colldir2;
974 
975  private:
977  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
978  CoalescedGeomDataSharedPtr pGeomData)
979  : Operator (pCollExp, pGeomData),
980  m_nquad0 (m_stdExp->GetNumPoints(0)),
981  m_nquad1 (m_stdExp->GetNumPoints(1)),
982  m_nquad2 (m_stdExp->GetNumPoints(2)),
983  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
984  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
985  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
986  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
987  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
988  m_colldir2(m_stdExp->GetBasis(2)->Collocation()),
989  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
990  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
991  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
992  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
993  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
994  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
995 
996  {
997  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
998  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
999  m_nmodes0*m_nmodes1*m_nmodes2));
1000  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1001  }
1002 };
1003 
1004 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Hex operator
1005 OperatorKey IProductWRTDerivBase_SumFac_Hex::m_type = GetOperatorFactory().
1006  RegisterCreatorFunction(
1008  IProductWRTDerivBase_SumFac_Hex::create,
1009  "IProductWRTDerivBase_SumFac_Hex");
1010 
1011 
1012 /**
1013  * @brief Inner product WRT deriv base operator using sum-factorisation (Tet)
1014  */
1016 {
1017  public:
1019 
1020  virtual void operator()(
1021  const Array<OneD, const NekDouble> &entry0,
1022  Array<OneD, NekDouble> &entry1,
1023  Array<OneD, NekDouble> &entry2,
1024  Array<OneD, NekDouble> &entry3,
1025  Array<OneD, NekDouble> &wsp)
1026  {
1027  unsigned int nPhys = m_stdExp->GetTotPoints();
1028  unsigned int ntot = m_numElmt*nPhys;
1029  unsigned int nmodes = m_stdExp->GetNcoeffs();
1030  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1032  Array<OneD, NekDouble> output, wsp1;
1034 
1035  in[0] = entry0; in[1] = entry1;
1036  in[2] = entry2;
1037 
1038  output = entry3;
1039 
1040  for(int i = 0; i < 3; ++i)
1041  {
1042  tmp[i] = wsp + i*nmax;
1043  }
1044 
1045 
1046  // calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
1047  // dphi/dxi_1 dxi_1/dx +
1048  // dphi/dxi_2 dxi_2/dx),in[0])
1049  // + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
1050  // dphi/dxi_1 dxi_1/dy +
1051  // dphi/dxi_2 dxi_2/dy),in[1])
1052  // + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
1053  // dphi/dxi_1 dxi_1/dz +
1054  // dphi/dxi_2 dxi_2/dz),in[1])
1055  //
1056  // Note dphi/dxi_0 =
1057  // dphi/deta_0 4/((1-eta_1)(1-eta2))
1058  //
1059  // dphi/dxi_1 =
1060  // dphi/deta_0 2(1+eta_0)/((1-eta_1)(1-eta_2)) +
1061  // dphi/deta_1 2/(1-eta_2)
1062  //
1063  // dphi/dxi_2 =
1064  // dphi/deta_0 2(1+eta_0)/((1-eta_1)(1-eta_2)) +
1065  // dphi/deta_1 (1+eta_1)/(1-eta_2) + dphi/deta_2
1066  //
1067  // and so the full inner products are
1068  //
1069  // (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
1070  // = (dphi/deta_0, fac0 (tmp0 + fac1(tmp1 + tmp2)))
1071  // + (dphi/deta_1, fac2 (tmp1 + fac3 tmp2))
1072  // + (dphi/deta_2, tmp2)
1073  //
1074  // tmp0 = (dxi_0/dx in[0] + dxi_0/dy in[1] + dxi_0/dz in[2])
1075  // tmp1 = (dxi_1/dx in[0] + dxi_1/dy in[1] + dxi_1/dz in[2])
1076  // tmp2 = (dxi_2/dx in[0] + dxi_2/dy in[1] + dxi_2/dz in[2])
1077 
1078  // fac0 = 4/((1-eta_1)(1-eta2))
1079  // fac1 = (1+eta_0)/2
1080  // fac2 = 2/(1-eta_2)
1081  // fac3 = (1+eta_1)/2
1082 
1083  for(int i = 0; i < 3; ++i)
1084  {
1085  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,tmp[i],1);
1086  for(int j = 1; j < 3; ++j)
1087  {
1088  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
1089  in[j],1, tmp[i], 1, tmp[i],1);
1090  }
1091  }
1092 
1093  wsp1 = wsp + 3*nmax;
1094 
1095  // Sort into eta factors
1096  for (int i = 0; i < m_numElmt; ++i)
1097  {
1098  // add tmp[1] + tmp[2]
1099  Vmath::Vadd(nPhys, tmp[1].get() + i*nPhys, 1,
1100  tmp[2].get() + i*nPhys, 1,
1101  wsp1.get(), 1);
1102 
1103  // scale wsp1 by fac1 and add to tmp0
1104  Vmath::Vvtvp(nPhys,&m_fac1[0],1,wsp1.get(),1,
1105  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1106 
1107  // scale tmp[0] by fac0
1108  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
1109  tmp[0].get()+i*nPhys,1);
1110 
1111  // scale tmp[2] by fac3 and add to tmp1
1112  Vmath::Vvtvp(nPhys,&m_fac3[0],1,tmp[2].get()+i*nPhys,1,
1113  tmp[1].get()+i*nPhys,1,tmp[1].get()+i*nPhys,1);
1114 
1115  // scale tmp[1] by fac2
1116  Vmath::Vmul(nPhys,&m_fac2[0],1,tmp[1].get()+i*nPhys,1,
1117  tmp[1].get()+i*nPhys,1);
1118  }
1119 
1120 
1121  // calculate Iproduct WRT Std Deriv
1122  TetIProduct(m_sortTopEdge, m_numElmt,
1123  m_nquad0, m_nquad1, m_nquad2,
1124  m_nmodes0, m_nmodes1, m_nmodes2,
1125  m_derbase0, m_base1, m_base2,
1126  m_jac,tmp[0],output,wsp1);
1127 
1128  TetIProduct(m_sortTopEdge, m_numElmt,
1129  m_nquad0, m_nquad1, m_nquad2,
1130  m_nmodes0, m_nmodes1, m_nmodes2,
1131  m_base0, m_derbase1, m_base2,
1132  m_jac,tmp[1],tmp[0],wsp1);
1133  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1134 
1135  TetIProduct(m_sortTopEdge, m_numElmt,
1136  m_nquad0, m_nquad1, m_nquad2,
1137  m_nmodes0, m_nmodes1, m_nmodes2,
1138  m_base0, m_base1, m_derbase2,
1139  m_jac,tmp[2],tmp[0],wsp1);
1140  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1141  }
1142 
1143  virtual void operator()(
1144  int dir,
1145  const Array<OneD, const NekDouble> &input,
1146  Array<OneD, NekDouble> &output,
1148  {
1149  ASSERTL0(false, "Not valid for this operator.");
1150  }
1151 
1152  protected:
1153  const int m_nquad0;
1154  const int m_nquad1;
1155  const int m_nquad2;
1156  const int m_nmodes0;
1157  const int m_nmodes1;
1158  const int m_nmodes2;
1172 
1173  private:
1175  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1176  CoalescedGeomDataSharedPtr pGeomData)
1177  : Operator (pCollExp, pGeomData),
1178  m_nquad0 (m_stdExp->GetNumPoints(0)),
1179  m_nquad1 (m_stdExp->GetNumPoints(1)),
1180  m_nquad2 (m_stdExp->GetNumPoints(2)),
1181  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1182  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1183  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1184  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1185  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1186  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1187  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1188  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1189  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1190 
1191  {
1192  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
1193  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1194  m_nmodes0*m_nmodes1*m_nmodes2));
1195  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1196 
1197 
1199  = m_stdExp->GetBasis(0)->GetZ();
1201  = m_stdExp->GetBasis(1)->GetZ();
1203  = m_stdExp->GetBasis(2)->GetZ();
1204 
1205  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1206  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1207  m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1208  m_fac3 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1209  // calculate 2.0/((1-eta_1)(1-eta_2))
1210  for (int i = 0; i < m_nquad0; ++i)
1211  {
1212  for(int j = 0; j < m_nquad1; ++j)
1213  {
1214  for(int k = 0; k < m_nquad2; ++k)
1215  {
1216 
1217  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1218  = 4.0/((1-z1[j])*(1-z2[k]));
1219  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1220  = (1+z0[i])*0.5;
1221  m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1222  = 2.0/(1-z2[k]);
1223  m_fac3[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1224  = (1+z1[j])*0.5;
1225  }
1226  }
1227  }
1228 
1229  if(m_stdExp->GetBasis(0)->GetBasisType()
1231  {
1232  m_sortTopEdge = true;
1233  }
1234  else
1235  {
1236  m_sortTopEdge = false;
1237  }
1238 
1239  }
1240 };
1241 
1242 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Tet operator
1243 OperatorKey IProductWRTDerivBase_SumFac_Tet::m_type = GetOperatorFactory().
1244  RegisterCreatorFunction(
1246  IProductWRTDerivBase_SumFac_Tet::create,
1247  "IProductWRTDerivBase_SumFac_Tet");
1248 
1249 
1250 /**
1251  * @brief Inner product WRT deriv base operator using sum-factorisation (Prism)
1252  */
1254 {
1255  public:
1257 
1259  {
1260  }
1261 
1262  virtual void operator()(
1263  const Array<OneD, const NekDouble> &entry0,
1264  Array<OneD, NekDouble> &entry1,
1265  Array<OneD, NekDouble> &entry2,
1266  Array<OneD, NekDouble> &entry3,
1268  {
1269  unsigned int nPhys = m_stdExp->GetTotPoints();
1270  unsigned int ntot = m_numElmt*nPhys;
1271  unsigned int nmodes = m_stdExp->GetNcoeffs();
1272  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1274  Array<OneD, NekDouble> output, wsp1;
1276 
1277  in[0] = entry0; in[1] = entry1;
1278  in[2] = entry2;
1279 
1280  output = entry3;
1281 
1282  for(int i = 0; i < 3; ++i)
1283  {
1284  tmp[i] = wsp + i*nmax;
1285  }
1286 
1287  // calculate (dphi/dx,in[0]) = ((dphi/dxi_0 dxi_0/dx +
1288  // dphi/dxi_1 dxi_1/dx),in[0])
1289  // + (dphi/dy,in[1]) = ((dphi/dxi_0 dxi_0/dy +
1290  // dphi/dxi_1 dxi_1/dy),in[1])
1291  // + (dphi/dz,in[2]) = ((dphi/dxi_0 dxi_0/dz +
1292  // dphi/dxi_1 dxi_1/dz),in[2])
1293  //
1294  // Note dphi/dxi_0 =
1295  // dphi/deta_0 deta_0/dxi_0 = dphi/deta_0 2/(1-eta_2)
1296  //
1297  // dphi/dxi_2 =
1298  // dphi/deta_0 deta_0/dxi_2 + dphi/deta_2 deta_2/dxi_2 =
1299  // dphi/deta_0 (1+eta_0)/(1-eta_2) + dphi/deta_2
1300  //
1301  // and so the full inner products are
1302  //
1303  // (dphi/dx,in[0]) + (dphi/dy,in[1]) + (dphi/dz,in[2])
1304  // = (dphi/deta_0, ((2/(1-eta_2) (dxi_0/dx in[0] + dxi_0/dy in[1]
1305  // + dxi_0/dz in[2])
1306  // + (1_eta_0)/(1-eta_2) (dxi_2/dx in[0] + dxi_2/dy in[1]
1307  // + dxi_2/dz in[2] ))
1308  // + (dphi/deta_1, (dxi_1/dx in[0] + dxi_1/dy in[1]
1309  // + dxi_1/dz in[2]))
1310  // + (dphi/deta_2, (dxi_2/dx in[0] + dxi_2/dy in[1]
1311  // + dxi_2/dz in[2]))
1312 
1313  for(int i = 0; i < 3; ++i)
1314  {
1315  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
1316  tmp[i],1);
1317  for(int j = 1; j < 3; ++j)
1318  {
1319  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
1320  in[j],1, tmp[i], 1, tmp[i],1);
1321  }
1322  }
1323  wsp1 = wsp + 3*nmax;
1324 
1325  // Sort into eta factors
1326  for (int i = 0; i < m_numElmt; ++i)
1327  {
1328  // scale tmp[0] by fac0
1329  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
1330  tmp[0].get()+i*nPhys,1);
1331 
1332  // scale tmp[2] by fac1 and add to tmp0
1333  Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[2].get()+i*nPhys,1,
1334  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1335  }
1336 
1337  // calculate Iproduct WRT Std Deriv
1338  PrismIProduct(m_sortTopVertex, m_numElmt,
1339  m_nquad0, m_nquad1, m_nquad2,
1340  m_nmodes0, m_nmodes1, m_nmodes2,
1341  m_derbase0, m_base1, m_base2,
1342  m_jac,tmp[0],output,wsp1);
1343 
1344  PrismIProduct(m_sortTopVertex, m_numElmt,
1345  m_nquad0, m_nquad1, m_nquad2,
1346  m_nmodes0, m_nmodes1, m_nmodes2,
1347  m_base0, m_derbase1, m_base2,
1348  m_jac,tmp[1],tmp[0],wsp1);
1349  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1350 
1351  PrismIProduct(m_sortTopVertex, m_numElmt,
1352  m_nquad0, m_nquad1, m_nquad2,
1353  m_nmodes0, m_nmodes1, m_nmodes2,
1354  m_base0, m_base1, m_derbase2,
1355  m_jac,tmp[2],tmp[0],wsp1);
1356  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1357  }
1358 
1359  virtual void operator()(
1360  int dir,
1361  const Array<OneD, const NekDouble> &input,
1362  Array<OneD, NekDouble> &output,
1364  {
1365  ASSERTL0(false, "Not valid for this operator.");
1366  }
1367 
1368 
1369  protected:
1370  const int m_nquad0;
1371  const int m_nquad1;
1372  const int m_nquad2;
1373  const int m_nmodes0;
1374  const int m_nmodes1;
1375  const int m_nmodes2;
1387 
1388  private:
1390  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1391  CoalescedGeomDataSharedPtr pGeomData)
1392  : Operator (pCollExp, pGeomData),
1393  m_nquad0 (m_stdExp->GetNumPoints(0)),
1394  m_nquad1 (m_stdExp->GetNumPoints(1)),
1395  m_nquad2 (m_stdExp->GetNumPoints(2)),
1396  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1397  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1398  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1399  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1400  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1401  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1402  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1403  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1404  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1405 
1406  {
1407  m_jac = pGeomData->GetJacWithStdWeights(pCollExp);
1408  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1409  m_nmodes0*m_nmodes1*m_nmodes2));
1410  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1411 
1412  if(m_stdExp->GetBasis(0)->GetBasisType()
1414  {
1415  m_sortTopVertex = true;
1416  }
1417  else
1418  {
1419  m_sortTopVertex = false;
1420  }
1421 
1423  = m_stdExp->GetBasis(0)->GetZ();
1425  = m_stdExp->GetBasis(2)->GetZ();
1426 
1427  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1428  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1429 
1430  for (int i = 0; i < m_nquad0; ++i)
1431  {
1432  for(int j = 0; j < m_nquad1; ++j)
1433  {
1434  for(int k = 0; k < m_nquad2; ++k)
1435  {
1436  // set up geometric factor: 2/(1-z1)
1437  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1438  = 2.0/(1-z2[k]);
1439  // set up geometric factor: (1+z0)/(1-z1)
1440  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1441  = (1+z0[i])/(1-z2[k]);
1442 
1443  }
1444  }
1445  }
1446  }
1447 };
1448 
1449 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Prism operator
1450 OperatorKey IProductWRTDerivBase_SumFac_Prism::m_type = GetOperatorFactory().
1451  RegisterCreatorFunction(
1453  IProductWRTDerivBase_SumFac_Prism::create,
1454  "IProductWRTDerivBase_SumFac_Prism");
1455 
1456 }
1457 }
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Inner product WRT deriv base operator using sum-factorisation (Tri)
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Inner product WRT deriv base operator using sum-factorisation (Segment)
void HexIProduct(bool colldir0, bool colldir1, bool colldir2, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:179
IProductWRTDerivBase_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
boost::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:159
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
IProductWRTDerivBase_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:442
Principle Modified Functions .
Definition: BasisType.h:49
Inner product WRT deriv base operator using LocalRegions implementation.
Inner product WRT deriv base operator using standard matrix approach.
STL namespace.
IProductWRTDerivBase_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
void QuadIProduct(bool colldir0, bool colldir1, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:49
Base class for operators on a collection of elements.
Definition: Operator.h:108
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
void TetIProduct(bool sortTopEdge, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:428
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
IProductWRTDerivBase_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
void TriIProduct(bool sortTopVertex, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:135
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
Inner product WRT deriv base operator using sum-factorisation (Tet)
virtual void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp)
Perform operation.
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
IProductWRTDerivBase_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
IProductWRTDerivBase_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Inner product WRT deriv base operator using element-wise operation.
Inner product WRT deriv base operator using sum-factorisation (Quad)
IProductWRTDerivBase_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Inner product WRT deriv base operator using sum-factorisation (Hex)
double NekDouble
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:110
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.
void PrismIProduct(bool sortTopVertex, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:356
Inner product WRT deriv base operator using sum-factorisation (Prism)
#define OPERATOR_CREATE(cname)
Definition: Operator.h:44
IProductWRTDerivBase_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
boost::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1061
IProductWRTDerivBase_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
vector< StdRegions::StdExpansionSharedPtr > m_expList
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
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215