Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: IProductWRTDerivBase operator implementations
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
37 #include <MatrixFreeOps/Operator.hpp>
38 
39 #include <Collections/Operator.h>
40 #include <Collections/MatrixFreeBase.h>
41 #include <Collections/Collection.h>
42 #include <Collections/IProduct.h>
43 
44 using namespace std;
45 
46 namespace Nektar {
47 namespace Collections {
48 
49 using LibUtilities::eSegment;
50 using LibUtilities::eQuadrilateral;
51 using LibUtilities::eTriangle;
52 using LibUtilities::eHexahedron;
53 using LibUtilities::eTetrahedron;
54 using LibUtilities::ePrism;
55 using LibUtilities::ePyramid;
56 
57 /**
58  * @brief Inner product WRT deriv base operator using standard matrix approach
59  */
60 class IProductWRTDerivBase_StdMat : public Operator
61 {
62  public:
63  OPERATOR_CREATE(IProductWRTDerivBase_StdMat)
64 
65  ~IProductWRTDerivBase_StdMat() final
66  {
67  }
68 
69  void operator()(
70  const Array<OneD, const NekDouble> &entry0,
71  Array<OneD, NekDouble> &entry1,
72  Array<OneD, NekDouble> &entry2,
73  Array<OneD, NekDouble> &entry3,
74  Array<OneD, NekDouble> &wsp) final
75  {
76 
77  int nPhys = m_stdExp->GetTotPoints();
78  int ntot = m_numElmt*nPhys;
79  int nmodes = m_stdExp->GetNcoeffs();
80  Array<OneD, Array<OneD, const NekDouble> > in(3);
81  Array<OneD, NekDouble> output;
82  Array<OneD, Array<OneD, NekDouble> > tmp(3);
83 
84  in[0] = entry0; in[1] = entry1;
85  in[2] = entry2;
86 
87  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
88  entry2: entry1;
89 
90  for(int i = 0; i < m_dim; ++i)
91  {
92  tmp[i] = wsp + i*ntot;
93  }
94 
95  // calculate Iproduct WRT Std Deriv
96 
97  // First component
98  if(m_isDeformed)
99  {
100  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
101  for(int i = 0; i < m_dim; ++i)
102  {
103  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
104  tmp[i],1);
105  for(int j = 1; j < m_coordim; ++j)
106  {
107  Vmath::Vvtvp (ntot,m_derivFac[i +j*m_dim],1,
108  in[j],1, tmp[i], 1, tmp[i],1);
109  }
110  }
111 
112  Vmath::Vmul(ntot,m_jac,1,tmp[0],1,tmp[0],1);
113  }
114  else
115  {
116  Array<OneD,NekDouble> t;
117  for(int e = 0; e < m_numElmt; ++e)
118  {
119  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
120  for(int i = 0; i < m_dim; ++i)
121  {
122  Vmath::Smul (m_nqe,m_derivFac[i][e],
123  in[0] + e*m_nqe,1,
124  t = tmp[i] + e*m_nqe,1);
125  for(int j = 1; j < m_coordim; ++j)
126  {
127  Vmath::Svtvp (m_nqe,m_derivFac[i +j*m_dim][e],
128  in[j] + e*m_nqe,1,
129  tmp[i] + e*m_nqe, 1,
130  t = tmp[i] + e*m_nqe,1);
131  }
132  }
133 
134  Vmath::Smul(m_nqe,m_jac[e],tmp[0]+e*m_nqe,1,t=tmp[0]+e*m_nqe,1);
135  }
136  }
137 
138  Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[0]->GetRows(),
139  m_numElmt,m_iProdWRTStdDBase[0]->GetColumns(),
140  1.0, m_iProdWRTStdDBase[0]->GetRawPtr(),
141  m_iProdWRTStdDBase[0]->GetRows(),
142  tmp[0].get(), nPhys, 0.0,
143  output.get(), nmodes);
144 
145  // Other components
146  for(int i = 1; i < m_dim; ++i)
147  {
148  if(m_isDeformed)
149  {
150  Vmath::Vmul(ntot,m_jac,1,tmp[i],1,tmp[i],1);
151  }
152  else
153  {
154  Array<OneD,NekDouble> t;
155  for(int e = 0; e < m_numElmt; ++e)
156  {
157  Vmath::Smul(m_nqe,m_jac[e],tmp[i]+e*m_nqe,1,
158  t=tmp[i]+e*m_nqe,1);
159  }
160  }
161  Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[i]->GetRows(),
162  m_numElmt,m_iProdWRTStdDBase[i]->GetColumns(),
163  1.0, m_iProdWRTStdDBase[i]->GetRawPtr(),
164  m_iProdWRTStdDBase[i]->GetRows(),
165  tmp[i].get(), nPhys, 1.0,
166  output.get(), nmodes);
167  }
168  }
169 
170  void operator()(int dir,
171  const Array<OneD, const NekDouble> &input,
172  Array<OneD, NekDouble> &output,
173  Array<OneD, NekDouble> &wsp) final
174  {
175  boost::ignore_unused(dir, input, output, wsp);
176  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
177  }
178 
179  virtual void CheckFactors(StdRegions::FactorMap factors,
180  int coll_phys_offset)
181  {
182  boost::ignore_unused(factors, coll_phys_offset);
183  ASSERTL0(false, "Not valid for this operator.");
184  }
185 
186  protected:
187  Array<OneD, DNekMatSharedPtr> m_iProdWRTStdDBase;
188  Array<TwoD, const NekDouble> m_derivFac;
189  Array<OneD, const NekDouble> m_jac;
190  int m_dim;
191  int m_coordim;
192 
193  private:
194  IProductWRTDerivBase_StdMat(
195  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
196  CoalescedGeomDataSharedPtr pGeomData,
197  StdRegions::FactorMap factors)
198  : Operator(pCollExp, pGeomData, factors)
199  {
200  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
201  m_dim = PtsKey.size();
202  m_coordim = pCollExp[0]->GetCoordim();
203 
204  m_nqe = m_stdExp->GetTotPoints();
205  int nmodes = m_stdExp->GetNcoeffs();
206 
207  // set up a IProductWRTDerivBase StdMat.
208  m_iProdWRTStdDBase = Array<OneD, DNekMatSharedPtr>(m_dim);
209  for(int i = 0; i < m_dim; ++i)
210  {
211  Array<OneD, NekDouble> tmp(m_nqe),tmp1(nmodes);
212  m_iProdWRTStdDBase[i] = MemoryManager<DNekMat>
213  ::AllocateSharedPtr(nmodes,m_nqe);
214  for(int j = 0; j < m_nqe; ++j)
215  {
216  Vmath::Zero(m_nqe,tmp,1);
217  tmp[j] = 1.0;
218  m_stdExp->IProductWRTDerivBase(i,tmp,tmp1);
219  Vmath::Vcopy(nmodes, &tmp1[0],1,
220  &(m_iProdWRTStdDBase[i]->GetPtr())[0]
221  + j*nmodes, 1);
222  }
223  }
224  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
225  m_jac = pGeomData->GetJac(pCollExp);
226  m_wspSize = m_dim*m_nqe*m_numElmt;
227  }
228 };
229 
230 /// Factory initialisation for the IProductWRTDerivBase_StdMat operators
231 OperatorKey IProductWRTDerivBase_StdMat::m_typeArr[] = {
232  GetOperatorFactory().RegisterCreatorFunction(
233  OperatorKey(eSegment, eIProductWRTDerivBase, eStdMat, false),
234  IProductWRTDerivBase_StdMat::create,
235  "IProductWRTDerivBase_StdMat_Seg"),
236  GetOperatorFactory().RegisterCreatorFunction(
237  OperatorKey(eTriangle, eIProductWRTDerivBase, eStdMat, false),
238  IProductWRTDerivBase_StdMat::create,
239  "IProductWRTDerivBase_StdMat_Tri"),
240  GetOperatorFactory().RegisterCreatorFunction(
241  OperatorKey(eTriangle, eIProductWRTDerivBase, eStdMat, true),
242  IProductWRTDerivBase_StdMat::create,
243  "IProductWRTDerivBase_StdMat_NodalTri"),
244  GetOperatorFactory().RegisterCreatorFunction(
245  OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eStdMat, false),
246  IProductWRTDerivBase_StdMat::create,
247  "IProductWRTDerivBase_StdMat_Quad"),
248  GetOperatorFactory().RegisterCreatorFunction(
249  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eStdMat, false),
250  IProductWRTDerivBase_StdMat::create,
251  "IProductWRTDerivBase_StdMat_Tet"),
252  GetOperatorFactory().RegisterCreatorFunction(
253  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eStdMat, true),
254  IProductWRTDerivBase_StdMat::create,
255  "IProductWRTDerivBase_StdMat_NodalTet"),
256  GetOperatorFactory().RegisterCreatorFunction(
257  OperatorKey(ePyramid, eIProductWRTDerivBase, eStdMat, false),
258  IProductWRTDerivBase_StdMat::create,
259  "IProductWRTDerivBase_StdMat_Pyr"),
260  GetOperatorFactory().RegisterCreatorFunction(
261  OperatorKey(ePrism, eIProductWRTDerivBase, eStdMat, false),
262  IProductWRTDerivBase_StdMat::create,
263  "IProductWRTDerivBase_StdMat_Prism"),
264  GetOperatorFactory().RegisterCreatorFunction(
265  OperatorKey(ePrism, eIProductWRTDerivBase, eStdMat, true),
266  IProductWRTDerivBase_StdMat::create,
267  "IProductWRTDerivBase_StdMat_NodalPrism"),
268  GetOperatorFactory().RegisterCreatorFunction(
269  OperatorKey(eHexahedron, eIProductWRTDerivBase, eStdMat, false),
270  IProductWRTDerivBase_StdMat::create,
271  "IProductWRTDerivBase_StdMat_Hex"),
272  GetOperatorFactory().RegisterCreatorFunction(
273  OperatorKey(ePyramid, eIProductWRTDerivBase, eSumFac, false),
274  IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_SumFac_Pyr")
275 };
276 
277 
278 /**
279  * @brief Inner product operator using operator using matrix free operators.
280  */
281 class IProductWRTDerivBase_MatrixFree : public Operator, MatrixFreeMultiInOneOut
282 {
283  public:
284  OPERATOR_CREATE(IProductWRTDerivBase_MatrixFree)
285 
286  ~IProductWRTDerivBase_MatrixFree() final
287  {
288  }
289 
290  void operator()(
291  const Array<OneD, const NekDouble> &entry0,
292  Array<OneD, NekDouble> &entry1,
293  Array<OneD, NekDouble> &entry2,
294  Array<OneD, NekDouble> &entry3,
295  Array<OneD, NekDouble> &wsp) final
296  {
297  boost::ignore_unused(wsp);
298 
299  Array<OneD, NekDouble> output;
300 
301  if (m_isPadded)
302  {
303  // copy into padded vector
304  switch(m_coordim)
305  {
306  case 1:
307  output = entry1;
308  Vmath::Vcopy(m_nIn, entry0, 1, m_input[0], 1);
309  break;
310  case 2:
311  output = entry2;
312  Vmath::Vcopy(m_nIn, entry0, 1, m_input[0], 1);
313  Vmath::Vcopy(m_nIn, entry1, 1, m_input[1], 1);
314  break;
315  case 3:
316  output = entry3;
317  Vmath::Vcopy(m_nIn, entry0, 1, m_input[0], 1);
318  Vmath::Vcopy(m_nIn, entry1, 1, m_input[1], 1);
319  Vmath::Vcopy(m_nIn, entry2, 1, m_input[2], 1);
320  break;
321  default:
322  NEKERROR(ErrorUtil::efatal, "m_coordim value not valid");
323  break;
324  }
325 
326  // call op
327  (*m_oper)(m_input, m_output);
328  // copy out of padded vector
329  Vmath::Vcopy(m_nOut, m_output, 1, output, 1);
330  }
331  else
332  {
333  Array<OneD, Array<OneD, NekDouble>> input{m_coordim};
334 
335  // copy into padded vector
336  switch(m_coordim)
337  {
338  case 1:
339  output = entry1;
340  input[0] = entry0;
341  break;
342  case 2:
343  output = entry2;
344  input[0] = entry0;
345  input[1] = entry1;
346  break;
347  case 3:
348  output = entry3;
349  input[0] = entry0;
350  input[1] = entry1;
351  input[2] = entry2;
352  break;
353  default:
354  NEKERROR(ErrorUtil::efatal, "coordim not valid");
355  break;
356  }
357 
358  (*m_oper)(input, output);
359  }
360  }
361 
362  void operator()(int dir,
363  const Array<OneD, const NekDouble> &input,
364  Array<OneD, NekDouble> &output,
365  Array<OneD, NekDouble> &wsp) final
366  {
367  boost::ignore_unused(dir, input, output, wsp);
368  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
369  }
370 
371  virtual void CheckFactors(StdRegions::FactorMap factors,
372  int coll_phys_offset)
373  {
374  boost::ignore_unused(factors, coll_phys_offset);
375  ASSERTL0(false, "Not valid for this operator.");
376  }
377 
378  private:
379  std::shared_ptr<MatrixFree::IProductWRTDerivBase> m_oper;
380 
381  IProductWRTDerivBase_MatrixFree(
382  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
383  CoalescedGeomDataSharedPtr pGeomData,
384  StdRegions::FactorMap factors)
385  : Operator(pCollExp, pGeomData, factors),
386  MatrixFreeMultiInOneOut(pCollExp[0]->GetCoordim(),
387  pCollExp[0]->GetStdExp()->GetTotPoints(),
388  pCollExp[0]->GetStdExp()->GetNcoeffs(),
389  pCollExp.size())
390  {
391  // Check if deformed
392  const auto dim = pCollExp[0]->GetStdExp()->GetShapeDimension();
393 
394  // Basis vector
395  std::vector<LibUtilities::BasisSharedPtr> basis(dim);
396  for (unsigned int i = 0; i < dim; ++i)
397  {
398  basis[i] = pCollExp[0]->GetBasis(i);
399  }
400 
401  // Get shape type
402  auto shapeType = pCollExp[0]->GetStdExp()->DetShapeType();
403 
404  // Generate operator string and create operator.
405  std::string op_string = "IProductWRTDerivBase";
406  op_string += MatrixFree::GetOpstring(shapeType, m_isDeformed);
407  auto oper = MatrixFree::GetOperatorFactory().
408  CreateInstance(op_string, basis, m_nElmtPad);
409 
410  // Set Jacobian
411  oper->SetJac(pGeomData->GetJacInterLeave(pCollExp,m_nElmtPad));
412 
413  // Set derivative factors
414  oper->SetDF(pGeomData->GetDerivFactorsInterLeave
415  (pCollExp,m_nElmtPad));
416 
417  m_oper = std::dynamic_pointer_cast<MatrixFree::
418  IProductWRTDerivBase>(oper);
419  ASSERTL0(m_oper, "Failed to cast pointer.");
420 
421  }
422 };
423 
424 /// Factory initialisation for the IProductWRTDerivBase_MatrixFree operators
425 OperatorKey IProductWRTDerivBase_MatrixFree::m_typeArr[] = {
426  GetOperatorFactory().RegisterCreatorFunction(
427  OperatorKey(eSegment, eIProductWRTDerivBase, eMatrixFree, false),
428  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Seg"),
429  GetOperatorFactory().RegisterCreatorFunction(
430  OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eMatrixFree, false),
431  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Quad"),
432  GetOperatorFactory().RegisterCreatorFunction(
433  OperatorKey(eTriangle, eIProductWRTDerivBase, eMatrixFree, false),
434  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Tri"),
435  GetOperatorFactory().RegisterCreatorFunction(
436  OperatorKey(eHexahedron, eIProductWRTDerivBase, eMatrixFree, false),
437  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Hex"),
438  GetOperatorFactory().RegisterCreatorFunction(
439  OperatorKey(ePrism, eIProductWRTDerivBase, eMatrixFree, false),
440  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Prism"),
441  GetOperatorFactory().RegisterCreatorFunction(
442  OperatorKey(ePyramid, eIProductWRTDerivBase, eMatrixFree, false),
443  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Pyr"),
444  GetOperatorFactory().RegisterCreatorFunction(
445  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eMatrixFree, false),
446  IProductWRTDerivBase_MatrixFree::create, "IProductWRTDerivBase_MatrixFree_Tet")
447 };
448 
449 /**
450  * @brief Inner product WRT deriv base operator using element-wise operation
451  */
452 class IProductWRTDerivBase_IterPerExp : public Operator
453 {
454  public:
455  OPERATOR_CREATE(IProductWRTDerivBase_IterPerExp)
456 
457  ~IProductWRTDerivBase_IterPerExp() final
458  {
459  }
460 
461  void operator()(
462  const Array<OneD, const NekDouble> &entry0,
463  Array<OneD, NekDouble> &entry1,
464  Array<OneD, NekDouble> &entry2,
465  Array<OneD, NekDouble> &entry3,
466  Array<OneD, NekDouble> &wsp) final
467  {
468 
469  unsigned int nPhys = m_stdExp->GetTotPoints();
470  unsigned int ntot = m_numElmt*nPhys;
471  unsigned int nmodes = m_stdExp->GetNcoeffs();
472  unsigned int nmax = max(ntot,m_numElmt*nmodes);
473  Array<OneD, Array<OneD, const NekDouble> > in(3);
474  Array<OneD, NekDouble> output, tmp1;
475  Array<OneD, Array<OneD, NekDouble> > tmp(3);
476 
477  in[0] = entry0; in[1] = entry1; in[2] = entry2;
478 
479  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
480  entry2: entry1;
481 
482  for(int i = 0; i < m_dim; ++i)
483  {
484  tmp[i] = wsp + i*nmax;
485  }
486 
487 
488  // calculate Iproduct WRT Std Deriv
489  // first component
490  if(m_isDeformed)
491  {
492  // calculate dx/dxi in[0] + dy/dxi in[2] + dz/dxi in[3]
493  for(int i = 0; i < m_dim; ++i)
494  {
495  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
496  tmp[i],1);
497  for(int j = 1; j < m_coordim; ++j)
498  {
499  Vmath::Vvtvp (ntot,m_derivFac[i +j*m_dim],1,
500  in[j],1, tmp[i], 1, tmp[i],1);
501  }
502  }
503 
504  Vmath::Vmul(ntot,m_jac,1,tmp[0],1,tmp[0],1);
505  }
506  else
507  {
508  Array<OneD,NekDouble> t;
509  for(int e = 0; e < m_numElmt; ++e)
510  {
511  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
512  for(int i = 0; i < m_dim; ++i)
513  {
514  Vmath::Smul (m_nqe,m_derivFac[i][e],
515  in[0] + e*m_nqe,1,
516  t = tmp[i] + e*m_nqe,1);
517  for(int j = 1; j < m_coordim; ++j)
518  {
519  Vmath::Svtvp (m_nqe,m_derivFac[i +j*m_dim][e],
520  in[j] + e*m_nqe,1,
521  tmp[i] + e*m_nqe, 1,
522  t = tmp[i] + e*m_nqe,1);
523  }
524  }
525 
526  Vmath::Smul(m_nqe,m_jac[e],tmp[0]+e*m_nqe,1,t=tmp[0]+e*m_nqe,1);
527  }
528  }
529 
530  for(int n = 0; n < m_numElmt; ++n)
531  {
532  m_stdExp->IProductWRTDerivBase(0,tmp[0]+n*nPhys,
533  tmp1 = output + n*nmodes);
534  }
535 
536  // other components
537  for(int i = 1; i < m_dim; ++i)
538  {
539  // multiply by Jacobian
540  if(m_isDeformed)
541  {
542  Vmath::Vmul(ntot,m_jac,1,tmp[i],1,tmp[i],1);
543  }
544  else
545  {
546  Array<OneD,NekDouble> t;
547  for(int e = 0; e < m_numElmt; ++e)
548  {
549  Vmath::Smul(m_nqe,m_jac[e],tmp[i]+e*m_nqe,1,
550  t=tmp[i]+e*m_nqe,1);
551  }
552  }
553 
554  for(int n = 0; n < m_numElmt; ++n)
555  {
556  m_stdExp->IProductWRTDerivBase(i,tmp[i]+n*nPhys,tmp[0]);
557  Vmath::Vadd(nmodes,tmp[0],1,output+n*nmodes,1,
558  tmp1 = output+n*nmodes,1);
559  }
560  }
561  }
562 
563  void operator()(int dir,
564  const Array<OneD, const NekDouble> &input,
565  Array<OneD, NekDouble> &output,
566  Array<OneD, NekDouble> &wsp) final
567  {
568  boost::ignore_unused(dir, input, output, wsp);
569  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
570  }
571 
572  virtual void CheckFactors(StdRegions::FactorMap factors,
573  int coll_phys_offset)
574  {
575  boost::ignore_unused(factors, coll_phys_offset);
576  ASSERTL0(false, "Not valid for this operator.");
577  }
578 
579  protected:
580  Array<TwoD, const NekDouble> m_derivFac;
581  Array<OneD, const NekDouble> m_jac;
582  int m_dim;
583  int m_coordim;
584 
585  private:
586  IProductWRTDerivBase_IterPerExp(
587  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
588  CoalescedGeomDataSharedPtr pGeomData,
589  StdRegions::FactorMap factors)
590  : Operator(pCollExp, pGeomData, factors)
591  {
592  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
593  m_dim = PtsKey.size();
594  m_coordim = pCollExp[0]->GetCoordim();
595 
596  m_nqe = m_stdExp->GetTotPoints();
597 
598  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
599  m_jac = pGeomData->GetJac(pCollExp);
600  m_wspSize = m_dim*m_nqe*m_numElmt;
601  }
602 };
603 
604 /// Factory initialisation for the IProductWRTDerivBase_IterPerExp operators
605 OperatorKey IProductWRTDerivBase_IterPerExp::m_typeArr[] = {
606  GetOperatorFactory().RegisterCreatorFunction(
607  OperatorKey(eSegment, eIProductWRTDerivBase, eIterPerExp,false),
608  IProductWRTDerivBase_IterPerExp::create,
609  "IProductWRTDerivBase_IterPerExp_Seg"),
610  GetOperatorFactory().RegisterCreatorFunction(
611  OperatorKey(eTriangle, eIProductWRTDerivBase, eIterPerExp,false),
612  IProductWRTDerivBase_IterPerExp::create,
613  "IProductWRTDerivBase_IterPerExp_Tri"),
614  GetOperatorFactory().RegisterCreatorFunction(
615  OperatorKey(eTriangle, eIProductWRTDerivBase, eIterPerExp,true),
616  IProductWRTDerivBase_IterPerExp::create,
617  "IProductWRTDerivBase_IterPerExp_NodalTri"),
618  GetOperatorFactory().RegisterCreatorFunction(
619  OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eIterPerExp,false),
620  IProductWRTDerivBase_IterPerExp::create,
621  "IProductWRTDerivBase_IterPerExp_Quad"),
622  GetOperatorFactory().RegisterCreatorFunction(
623  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eIterPerExp,false),
624  IProductWRTDerivBase_IterPerExp::create,
625  "IProductWRTDerivBase_IterPerExp_Tet"),
626  GetOperatorFactory().RegisterCreatorFunction(
627  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eIterPerExp,true),
628  IProductWRTDerivBase_IterPerExp::create,
629  "IProductWRTDerivBase_IterPerExp_NodalTet"),
630  GetOperatorFactory().RegisterCreatorFunction(
631  OperatorKey(ePyramid, eIProductWRTDerivBase, eIterPerExp,false),
632  IProductWRTDerivBase_IterPerExp::create,
633  "IProductWRTDerivBase_IterPerExp_Pyr"),
634  GetOperatorFactory().RegisterCreatorFunction(
635  OperatorKey(ePrism, eIProductWRTDerivBase, eIterPerExp,false),
636  IProductWRTDerivBase_IterPerExp::create,
637  "IProductWRTDerivBase_IterPerExp_Prism"),
638  GetOperatorFactory().RegisterCreatorFunction(
639  OperatorKey(ePrism, eIProductWRTDerivBase, eIterPerExp,true),
640  IProductWRTDerivBase_IterPerExp::create,
641  "IProductWRTDerivBase_IterPerExp_NodalPrism"),
642  GetOperatorFactory().RegisterCreatorFunction(
643  OperatorKey(eHexahedron, eIProductWRTDerivBase, eIterPerExp,false),
644  IProductWRTDerivBase_IterPerExp::create,
645  "IProductWRTDerivBase_IterPerExp_Hex")
646 };
647 
648 
649 /**
650  * @brief Inner product WRT deriv base operator using LocalRegions
651  * implementation.
652  */
653 class IProductWRTDerivBase_NoCollection : public Operator
654 {
655  public:
656  OPERATOR_CREATE(IProductWRTDerivBase_NoCollection)
657 
658  ~IProductWRTDerivBase_NoCollection() final
659  {
660  }
661 
662  void operator()(
663  const Array<OneD, const NekDouble> &entry0,
664  Array<OneD, NekDouble> &entry1,
665  Array<OneD, NekDouble> &entry2,
666  Array<OneD, NekDouble> &entry3,
667  Array<OneD, NekDouble> &wsp) final
668  {
669  boost::ignore_unused(wsp);
670 
671  unsigned int nmodes = m_expList[0]->GetNcoeffs();
672  unsigned int nPhys = m_expList[0]->GetTotPoints();
673  Array<OneD, NekDouble> tmp(nmodes),tmp1;
674 
675  Array<OneD, Array<OneD, const NekDouble> > in(3);
676  Array<OneD, NekDouble> output;
677  in[0] = entry0; in[1] = entry1; in[2] = entry2;
678 
679  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
680  entry2: entry1;
681 
682  for(int n = 0; n < m_numElmt; ++n)
683  {
684  m_expList[n]->IProductWRTDerivBase(0,
685  in[0] + n * nPhys,
686  tmp1 = output + n * nmodes);
687  }
688 
689  for(int i = 1; i < m_dim; ++i)
690  {
691  for(int n = 0; n < m_numElmt; ++n)
692  {
693  m_expList[n]->IProductWRTDerivBase(i,in[i]+n*nPhys,tmp);
694 
695  Vmath::Vadd(nmodes,tmp,1,output+n*nmodes,1,
696  tmp1 = output+n*nmodes,1);
697  }
698  }
699  }
700 
701  void operator()(int dir,
702  const Array<OneD, const NekDouble> &input,
703  Array<OneD, NekDouble> &output,
704  Array<OneD, NekDouble> &wsp) final
705  {
706  boost::ignore_unused(dir, input, output, wsp);
707  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
708  }
709 
710  virtual void CheckFactors(StdRegions::FactorMap factors,
711  int coll_phys_offset)
712  {
713  boost::ignore_unused(factors, coll_phys_offset);
714  ASSERTL0(false, "Not valid for this operator.");
715  }
716 
717  protected:
718  int m_dim;
719  int m_coordim;
720  vector<StdRegions::StdExpansionSharedPtr> m_expList;
721 
722  private:
723  IProductWRTDerivBase_NoCollection(
724  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
725  CoalescedGeomDataSharedPtr pGeomData,
726  StdRegions::FactorMap factors)
727  : Operator(pCollExp, pGeomData, factors)
728  {
729  m_expList = pCollExp;
730  m_dim = pCollExp[0]->GetNumBases();
731  m_coordim = pCollExp[0]->GetCoordim();
732  }
733 };
734 
735 /// Factory initialisation for the IProductWRTDerivBase_NoCollection operators
736 OperatorKey IProductWRTDerivBase_NoCollection::m_typeArr[] = {
737  GetOperatorFactory().RegisterCreatorFunction(
738  OperatorKey(eSegment, eIProductWRTDerivBase, eNoCollection,false),
739  IProductWRTDerivBase_NoCollection::create,
740  "IProductWRTDerivBase_NoCollection_Seg"),
741  GetOperatorFactory().RegisterCreatorFunction(
742  OperatorKey(eTriangle, eIProductWRTDerivBase, eNoCollection,false),
743  IProductWRTDerivBase_NoCollection::create,
744  "IProductWRTDerivBase_NoCollection_Tri"),
745  GetOperatorFactory().RegisterCreatorFunction(
746  OperatorKey(eTriangle, eIProductWRTDerivBase, eNoCollection,true),
747  IProductWRTDerivBase_NoCollection::create,
748  "IProductWRTDerivBase_NoCollection_NodalTri"),
749  GetOperatorFactory().RegisterCreatorFunction(
750  OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eNoCollection,false),
751  IProductWRTDerivBase_NoCollection::create,
752  "IProductWRTDerivBase_NoCollection_Quad"),
753  GetOperatorFactory().RegisterCreatorFunction(
754  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eNoCollection,false),
755  IProductWRTDerivBase_NoCollection::create,
756  "IProductWRTDerivBase_NoCollection_Tet"),
757  GetOperatorFactory().RegisterCreatorFunction(
758  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eNoCollection,true),
759  IProductWRTDerivBase_NoCollection::create,
760  "IProductWRTDerivBase_NoCollection_NodalTet"),
761  GetOperatorFactory().RegisterCreatorFunction(
762  OperatorKey(ePyramid, eIProductWRTDerivBase, eNoCollection,false),
763  IProductWRTDerivBase_NoCollection::create,
764  "IProductWRTDerivBase_NoCollection_Pyr"),
765  GetOperatorFactory().RegisterCreatorFunction(
766  OperatorKey(ePrism, eIProductWRTDerivBase, eNoCollection,false),
767  IProductWRTDerivBase_NoCollection::create,
768  "IProductWRTDerivBase_NoCollection_Prism"),
769  GetOperatorFactory().RegisterCreatorFunction(
770  OperatorKey(ePrism, eIProductWRTDerivBase, eNoCollection,true),
771  IProductWRTDerivBase_NoCollection::create,
772  "IProductWRTDerivBase_NoCollection_NodalPrism"),
773  GetOperatorFactory().RegisterCreatorFunction(
774  OperatorKey(eHexahedron, eIProductWRTDerivBase, eNoCollection,false),
775  IProductWRTDerivBase_NoCollection::create,
776  "IProductWRTDerivBase_NoCollection_Hex")
777 };
778 
779 
780 /**
781  * @brief Inner product WRT deriv base operator using sum-factorisation
782  * (Segment)
783  */
784 class IProductWRTDerivBase_SumFac_Seg : public Operator
785 {
786  public:
787  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Seg)
788 
789  ~IProductWRTDerivBase_SumFac_Seg() final
790  {
791  }
792 
793  void operator()(
794  const Array<OneD, const NekDouble> &entry0,
795  Array<OneD, NekDouble> &entry1,
796  Array<OneD, NekDouble> &entry2,
797  Array<OneD, NekDouble> &entry3,
798  Array<OneD, NekDouble> &wsp) final
799  {
800  Array<OneD, Array<OneD, const NekDouble> > in(3);
801  Array<OneD, NekDouble> output;
802 
803  output = (m_coordim == 3)? entry3: (m_coordim == 2)?
804  entry2: entry1;
805 
806  in[0] = entry0; in[1] = entry1;
807  in[2] = entry2;
808 
809  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
810  if(m_isDeformed)
811  {
812  Vmath::Vmul (m_wspSize,m_derivFac[0], 1, in[0], 1, wsp, 1);
813  for(int i = 1; i < m_coordim; ++i)
814  {
815  Vmath::Vvtvp (m_wspSize,m_derivFac[i],1,in[i],1,wsp,1,wsp,1);
816  }
817  }
818  else
819  {
820  Array<OneD,NekDouble> t;
821  for(int e = 0; e < m_numElmt; ++e)
822  {
823  Vmath::Smul(m_nquad0, m_derivFac[0][e],
824  in[0] + e*m_nquad0, 1,
825  t = wsp + e*m_nquad0, 1);
826  }
827 
828  for(int i = 1; i < m_coordim; ++i)
829  {
830  for(int e = 0; e < m_numElmt; ++e)
831  {
832  Vmath::Svtvp (m_nquad0, m_derivFac[i][e],
833  in[i] + e*m_nquad0,1,
834  wsp + e*m_nquad0,1,
835  t = wsp + e*m_nquad0,1);
836  }
837  }
838  }
839 
840  Vmath::Vmul(m_wspSize, m_jacWStdW, 1, wsp, 1, wsp, 1);
841 
842  // out = B0*in;
843  Blas::Dgemm('T', 'N', m_nmodes0, m_numElmt, m_nquad0,
844  1.0, m_derbase0.get(), m_nquad0,
845  &wsp[0], m_nquad0, 0.0,
846  &output[0], m_nmodes0);
847  }
848 
849  void operator()(int dir,
850  const Array<OneD, const NekDouble> &input,
851  Array<OneD, NekDouble> &output,
852  Array<OneD, NekDouble> &wsp) final
853  {
854  boost::ignore_unused(dir, input, output, wsp);
855  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
856  }
857 
858  virtual void CheckFactors(StdRegions::FactorMap factors,
859  int coll_phys_offset)
860  {
861  boost::ignore_unused(factors, coll_phys_offset);
862  ASSERTL0(false, "Not valid for this operator.");
863  }
864 
865  protected:
866  const int m_nquad0;
867  const int m_nmodes0;
868  Array<OneD, const NekDouble> m_jacWStdW;
869  Array<OneD, const NekDouble> m_derbase0;
870  Array<TwoD, const NekDouble> m_derivFac;
871  int m_coordim;
872 
873  private:
874  IProductWRTDerivBase_SumFac_Seg(
875  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
876  CoalescedGeomDataSharedPtr pGeomData,
877  StdRegions::FactorMap factors)
878  : Operator(pCollExp, pGeomData, factors),
879  m_nquad0 (m_stdExp->GetNumPoints(0)),
880  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
881  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata())
882  {
883  m_coordim = pCollExp[0]->GetCoordim();
884  m_wspSize = m_numElmt*m_nquad0;
885  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
886  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
887  }
888 };
889 
890 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Seg operator
891 OperatorKey IProductWRTDerivBase_SumFac_Seg::m_type = GetOperatorFactory().
892  RegisterCreatorFunction(
893  OperatorKey(eSegment, eIProductWRTDerivBase, eSumFac, false),
894  IProductWRTDerivBase_SumFac_Seg::create,
895  "IProductWRTDerivBase_SumFac_Seg");
896 
897 
898 /**
899  * @brief Inner product WRT deriv base operator using sum-factorisation (Quad)
900  */
901 class IProductWRTDerivBase_SumFac_Quad : public Operator
902 {
903 public:
904  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Quad)
905 
906  ~IProductWRTDerivBase_SumFac_Quad() final
907  {
908  }
909 
910  void operator()(
911  const Array<OneD, const NekDouble> &entry0,
912  Array<OneD, NekDouble> &entry1,
913  Array<OneD, NekDouble> &entry2,
914  Array<OneD, NekDouble> &entry3,
915  Array<OneD, NekDouble> &wsp) final
916  {
917 
918  unsigned int nPhys = m_stdExp->GetTotPoints();
919  unsigned int ntot = m_numElmt*nPhys;
920  unsigned int nmodes = m_stdExp->GetNcoeffs();
921  unsigned int nmax = max(ntot,m_numElmt*nmodes);
922  Array<OneD, Array<OneD, const NekDouble> > in(3);
923  Array<OneD, NekDouble> output, wsp1;
924  Array<OneD, Array<OneD, NekDouble> > tmp(2);
925 
926  in[0] = entry0; in[1] = entry1; in[2] = entry2;
927 
928  output = (m_coordim == 2)? entry2: entry3;
929 
930  tmp[0] = wsp; tmp[1] = wsp + nmax;
931  wsp1 = wsp + 2*nmax;
932 
933  if(m_isDeformed)
934  {
935  // calculate dx/dxi in[0] + dy/dxi in[1]
936  for(int i = 0; i < 2; ++i)
937  {
938  Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
939  for(int j = 1; j < m_coordim; ++j)
940  {
941  Vmath::Vvtvp(ntot, m_derivFac[i + 2*j], 1, in[j], 1,
942  tmp[i], 1, tmp[i], 1);
943  }
944  }
945  }
946  else
947  {
948  Array<OneD,NekDouble> t;
949  for(int e = 0; e < m_numElmt; ++e)
950  {
951  // calculate dx/dxi in[0] + dy/dxi in[1]
952  for(int i = 0; i < 2; ++i)
953  {
954  Vmath::Smul (m_nqe,m_derivFac[i][e],
955  in[0] + e*m_nqe,1,
956  t = tmp[i] + e*m_nqe,1);
957  for(int j = 1; j < m_coordim; ++j)
958  {
959  Vmath::Svtvp (m_nqe,m_derivFac[i+2*j][e],
960  in[j] + e*m_nqe,1,
961  tmp[i] + e*m_nqe, 1,
962  t = tmp[i] + e*m_nqe,1);
963  }
964  }
965  }
966  }
967  // Iproduct wrt derivative of base 1
968  QuadIProduct(false, m_colldir1, m_numElmt,
969  m_nquad0, m_nquad1,
970  m_nmodes0, m_nmodes1,
971  m_derbase0, m_base1,
972  m_jacWStdW, tmp[0], output, wsp1);
973 
974  // Iproduct wrt derivative of base 1
975  QuadIProduct(m_colldir0, false, m_numElmt,
976  m_nquad0, m_nquad1,
977  m_nmodes0, m_nmodes1,
978  m_base0, m_derbase1,
979  m_jacWStdW, tmp[1], tmp[0], wsp1);
980 
981  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
982  }
983 
984  void operator()(int dir,
985  const Array<OneD, const NekDouble> &input,
986  Array<OneD, NekDouble> &output,
987  Array<OneD, NekDouble> &wsp) final
988  {
989  boost::ignore_unused(dir, input, output, wsp);
990  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
991  }
992 
993  virtual void CheckFactors(StdRegions::FactorMap factors,
994  int coll_phys_offset)
995  {
996  boost::ignore_unused(factors, coll_phys_offset);
997  ASSERTL0(false, "Not valid for this operator.");
998  }
999 
1000 protected:
1001  const int m_nquad0;
1002  const int m_nquad1;
1003  const int m_nmodes0;
1004  const int m_nmodes1;
1005  const bool m_colldir0;
1006  const bool m_colldir1;
1007  int m_coordim;
1008  Array<TwoD, const NekDouble> m_derivFac;
1009  Array<OneD, const NekDouble> m_jacWStdW;
1010  Array<OneD, const NekDouble> m_base0;
1011  Array<OneD, const NekDouble> m_base1;
1012  Array<OneD, const NekDouble> m_derbase0;
1013  Array<OneD, const NekDouble> m_derbase1;
1014 
1015 private:
1016  IProductWRTDerivBase_SumFac_Quad(
1017  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1018  CoalescedGeomDataSharedPtr pGeomData,
1019  StdRegions::FactorMap factors)
1020  : Operator(pCollExp, pGeomData, factors),
1021  m_nquad0 (m_stdExp->GetNumPoints(0)),
1022  m_nquad1 (m_stdExp->GetNumPoints(1)),
1023  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1024  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1025  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
1026  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
1027  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1028  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1029  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1030  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
1031  {
1032  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1033  m_coordim = pCollExp[0]->GetCoordim();
1034 
1035  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1036  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1037  m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
1038  m_nmodes0*m_nmodes1));
1039  }
1040 };
1041 
1042 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Quad operator
1043 OperatorKey IProductWRTDerivBase_SumFac_Quad::m_type =
1044  GetOperatorFactory().RegisterCreatorFunction(
1045  OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eSumFac, false),
1046  IProductWRTDerivBase_SumFac_Quad::create,
1047  "IProductWRTDerivBase_IterPerExp_Quad");
1048 
1049 
1050 /**
1051  * @brief Inner product WRT deriv base operator using sum-factorisation (Tri)
1052  */
1053 class IProductWRTDerivBase_SumFac_Tri : public Operator
1054 {
1055 public:
1056  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Tri)
1057 
1058  ~IProductWRTDerivBase_SumFac_Tri() final
1059  {
1060  }
1061 
1062  /**
1063  * This method calculates:
1064  *
1065  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) \f]
1066  *
1067  * which can be represented in terms of local cartesian
1068  * derivaties as:
1069  *
1070  * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1071  * d\phi/d\xi_1\, d\xi_1/dx),in[0]) + \f]
1072  *
1073  * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1074  * d\phi/d\xi_1\, d\xi_1/dy),in[1]) + \f]
1075  *
1076  * where we note that
1077  *
1078  * \f[ d\phi/d\xi_0 = d\phi/d\eta_0\, d\eta_0/d\xi_0 =
1079  * d\phi/d\eta_0 2/(1-\eta_1) \f]
1080  *
1081  * \f[ d\phi/d\xi_1 = d\phi/d\eta_1\, d\eta_1/d\xi_1 +
1082  * d\phi/d\eta_1\, d\eta_1/d\xi_1 = d\phi/d\eta_0 (1+\eta_0)/(1-\eta_1)
1083  * + d\phi/d\eta_1 \f]
1084  *
1085  * and so the full inner products are
1086  *
1087  * \f[ (d\phi/dx,in[0]) + (dphi/dy,in[1]) =
1088  * (d\phi/d\eta_0, ((2/(1-\eta_1) (d\xi_0/dx in[0] + d\xi_0/dy in[1])
1089  * + (1-\eta_0)/(1-\eta_1) (d\xi_1/dx in[0]+d\xi_1/dy in[1]))
1090  * + (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1])) \f]
1091  *
1092  */
1093  void operator()(
1094  const Array<OneD, const NekDouble> &entry0,
1095  Array<OneD, NekDouble> &entry1,
1096  Array<OneD, NekDouble> &entry2,
1097  Array<OneD, NekDouble> &entry3,
1098  Array<OneD, NekDouble> &wsp) final
1099  {
1100 
1101  unsigned int nPhys = m_stdExp->GetTotPoints();
1102  unsigned int ntot = m_numElmt*nPhys;
1103  unsigned int nmodes = m_stdExp->GetNcoeffs();
1104  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1105  Array<OneD, Array<OneD, const NekDouble> > in(3);
1106  Array<OneD, NekDouble> output, wsp1;
1107  Array<OneD, Array<OneD, NekDouble> > tmp(2);
1108 
1109  in[0] = entry0; in[1] = entry1; in[2] = entry2;
1110 
1111  output = (m_coordim == 2)? entry2: entry3;
1112 
1113  tmp[0] = wsp; tmp[1] = wsp + nmax;
1114  wsp1 = wsp + 2*nmax;
1115 
1116  if(m_isDeformed)
1117  {
1118  for (int i = 0; i < 2; ++i)
1119  {
1120  Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1121  for(int j = 1; j < m_coordim; ++j)
1122  {
1123  Vmath::Vvtvp(ntot, m_derivFac[i+2*j], 1, in[j], 1,
1124  tmp[i], 1, tmp[i], 1);
1125  }
1126  }
1127  }
1128  else
1129  {
1130  Array<OneD,NekDouble> t;
1131  for(int e = 0; e < m_numElmt; ++e)
1132  {
1133  // calculate dx/dxi in[0] + dy/dxi in[1]
1134  for(int i = 0; i < 2; ++i)
1135  {
1136  Vmath::Smul (m_nqe,m_derivFac[i][e],
1137  in[0] + e*m_nqe,1,
1138  t = tmp[i] + e*m_nqe,1);
1139  for(int j = 1; j < m_coordim; ++j)
1140  {
1141  Vmath::Svtvp (m_nqe,m_derivFac[i+2*j][e],
1142  in[j] + e*m_nqe,1,
1143  tmp[i] + e*m_nqe, 1,
1144  t = tmp[i] + e*m_nqe,1);
1145  }
1146  }
1147  }
1148  }
1149 
1150  // Multiply by factor: 2/(1-z1)
1151  for (int i = 0; i < m_numElmt; ++i)
1152  {
1153  // scale tmp[0] by geometric factor: 2/(1-z1)
1154  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
1155  tmp[0].get()+i*nPhys,1);
1156 
1157  // scale tmp[1] by geometric factor (1+z0)/(1-z1)
1158  Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[1].get()+i*nPhys,1,
1159  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1160  }
1161 
1162  // Iproduct wrt derivative of base 0
1163  TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
1164  m_nmodes0, m_nmodes1, m_derbase0, m_base1,
1165  m_jacWStdW, tmp[0], output, wsp1);
1166 
1167  // Iproduct wrt derivative of base 1
1168  TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1,
1169  m_nmodes0, m_nmodes1, m_base0, m_derbase1,
1170  m_jacWStdW, tmp[1], tmp[0], wsp1);
1171 
1172  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1173  }
1174 
1175  void operator()(int dir,
1176  const Array<OneD, const NekDouble> &input,
1177  Array<OneD, NekDouble> &output,
1178  Array<OneD, NekDouble> &wsp) final
1179  {
1180  boost::ignore_unused(dir, input, output, wsp);
1181  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1182  }
1183 
1184  virtual void CheckFactors(StdRegions::FactorMap factors,
1185  int coll_phys_offset)
1186  {
1187  boost::ignore_unused(factors, coll_phys_offset);
1188  ASSERTL0(false, "Not valid for this operator.");
1189  }
1190 
1191  protected:
1192  const int m_nquad0;
1193  const int m_nquad1;
1194  const int m_nmodes0;
1195  const int m_nmodes1;
1196  const bool m_colldir0;
1197  const bool m_colldir1;
1198  int m_coordim;
1199  Array<TwoD, const NekDouble> m_derivFac;
1200  Array<OneD, const NekDouble> m_jacWStdW;
1201  Array<OneD, const NekDouble> m_base0;
1202  Array<OneD, const NekDouble> m_base1;
1203  Array<OneD, const NekDouble> m_derbase0;
1204  Array<OneD, const NekDouble> m_derbase1;
1205  Array<OneD, NekDouble> m_fac0;
1206  Array<OneD, NekDouble> m_fac1;
1207  bool m_sortTopVertex;
1208 
1209  private:
1210  IProductWRTDerivBase_SumFac_Tri(
1211  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1212  CoalescedGeomDataSharedPtr pGeomData,
1213  StdRegions::FactorMap factors)
1214  : Operator(pCollExp, pGeomData, factors),
1215  m_nquad0 (m_stdExp->GetNumPoints(0)),
1216  m_nquad1 (m_stdExp->GetNumPoints(1)),
1217  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1218  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1219  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
1220  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
1221  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1222  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1223  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1224  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
1225  {
1226  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1227  m_coordim = pCollExp[0]->GetCoordim();
1228 
1229  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1230  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1231  m_wspSize = 4 * m_numElmt * (max(m_nquad0*m_nquad1,
1232  m_nmodes0*m_nmodes1));
1233 
1234  if(m_stdExp->GetBasis(0)->GetBasisType()
1235  == LibUtilities::eModified_A)
1236  {
1237  m_sortTopVertex = true;
1238  }
1239  else
1240  {
1241  m_sortTopVertex = false;
1242  }
1243 
1244  const Array<OneD, const NekDouble>& z0
1245  = m_stdExp->GetBasis(0)->GetZ();
1246  const Array<OneD, const NekDouble>& z1
1247  = m_stdExp->GetBasis(1)->GetZ();
1248 
1249  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
1250  // set up geometric factor: 2/(1-z1)
1251  for (int i = 0; i < m_nquad0; ++i)
1252  {
1253  for(int j = 0; j < m_nquad1; ++j)
1254  {
1255  m_fac0[i+j*m_nquad0] = 2.0/(1-z1[j]);
1256  }
1257  }
1258 
1259  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
1260  // set up geometric factor: (1+z0)/(1-z1)
1261  for (int i = 0; i < m_nquad0; ++i)
1262  {
1263  for(int j = 0; j < m_nquad1; ++j)
1264  {
1265  m_fac1[i+j*m_nquad0] = (1+z0[i])/(1-z1[j]);
1266  }
1267  }
1268  }
1269 };
1270 
1271 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Tri operator
1272 OperatorKey IProductWRTDerivBase_SumFac_Tri::m_type =
1273  GetOperatorFactory().RegisterCreatorFunction(
1274  OperatorKey(eTriangle, eIProductWRTDerivBase, eSumFac,false),
1275  IProductWRTDerivBase_SumFac_Tri::create,
1276  "IProductWRTDerivBase_IterPerExp_Tri");
1277 
1278 
1279 /**
1280  * @brief Inner product WRT deriv base operator using sum-factorisation (Hex)
1281  */
1282 class IProductWRTDerivBase_SumFac_Hex : public Operator
1283 {
1284  public:
1285  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Hex)
1286 
1287  ~IProductWRTDerivBase_SumFac_Hex() final
1288  {
1289  }
1290 
1291  void operator()(
1292  const Array<OneD, const NekDouble> &entry0,
1293  Array<OneD, NekDouble> &entry1,
1294  Array<OneD, NekDouble> &entry2,
1295  Array<OneD, NekDouble> &entry3,
1296  Array<OneD, NekDouble> &wsp) final
1297  {
1298 
1299  unsigned int nPhys = m_stdExp->GetTotPoints();
1300  unsigned int ntot = m_numElmt*nPhys;
1301  unsigned int nmodes = m_stdExp->GetNcoeffs();
1302  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1303  Array<OneD, Array<OneD, const NekDouble> > in(3);
1304  Array<OneD, NekDouble> output, wsp1;
1305  Array<OneD, Array<OneD, NekDouble> > tmp(3);
1306 
1307  in[0] = entry0; in[1] = entry1;
1308  in[2] = entry2;
1309 
1310  output = entry3;
1311 
1312  for(int i = 0; i < 3; ++i)
1313  {
1314  tmp[i] = wsp + i*nmax;
1315  }
1316 
1317  if(m_isDeformed)
1318  {
1319  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1320  for(int i = 0; i < 3; ++i)
1321  {
1322  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
1323  tmp[i],1);
1324  for(int j = 1; j < 3; ++j)
1325  {
1326  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
1327  in[j],1, tmp[i], 1, tmp[i],1);
1328  }
1329  }
1330  }
1331  else
1332  {
1333  Array<OneD,NekDouble> t;
1334  for(int e = 0; e < m_numElmt; ++e)
1335  {
1336  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1337  for(int i = 0; i < 3; ++i)
1338  {
1339  Vmath::Smul (m_nqe,m_derivFac[i][e],
1340  in[0] + e*m_nqe,1,
1341  t = tmp[i] + e*m_nqe,1);
1342  for(int j = 1; j < 3; ++j)
1343  {
1344  Vmath::Svtvp (m_nqe,m_derivFac[i+3*j][e],
1345  in[j] + e*m_nqe,1,
1346  tmp[i] + e*m_nqe, 1,
1347  t = tmp[i] + e*m_nqe,1);
1348  }
1349  }
1350  }
1351  }
1352 
1353  wsp1 = wsp + 3*nmax;
1354 
1355  // calculate Iproduct WRT Std Deriv
1356  HexIProduct(false,m_colldir1,m_colldir2, m_numElmt,
1357  m_nquad0, m_nquad1, m_nquad2,
1358  m_nmodes0, m_nmodes1, m_nmodes2,
1359  m_derbase0, m_base1, m_base2,
1360  m_jacWStdW,tmp[0],output,wsp1);
1361 
1362  HexIProduct(m_colldir0,false,m_colldir2, m_numElmt,
1363  m_nquad0, m_nquad1, m_nquad2,
1364  m_nmodes0, m_nmodes1, m_nmodes2,
1365  m_base0, m_derbase1, m_base2,
1366  m_jacWStdW,tmp[1],tmp[0],wsp1);
1367  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1368 
1369  HexIProduct(m_colldir0,m_colldir1,false, m_numElmt,
1370  m_nquad0, m_nquad1, m_nquad2,
1371  m_nmodes0, m_nmodes1, m_nmodes2,
1372  m_base0, m_base1, m_derbase2,
1373  m_jacWStdW,tmp[2],tmp[0],wsp1);
1374  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1375  }
1376 
1377  void operator()(int dir,
1378  const Array<OneD, const NekDouble> &input,
1379  Array<OneD, NekDouble> &output,
1380  Array<OneD, NekDouble> &wsp) final
1381  {
1382  boost::ignore_unused(dir, input, output, wsp);
1383  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1384  }
1385 
1386  virtual void CheckFactors(StdRegions::FactorMap factors,
1387  int coll_phys_offset)
1388  {
1389  boost::ignore_unused(factors, coll_phys_offset);
1390  ASSERTL0(false, "Not valid for this operator.");
1391  }
1392 
1393  protected:
1394  const int m_nquad0;
1395  const int m_nquad1;
1396  const int m_nquad2;
1397  const int m_nmodes0;
1398  const int m_nmodes1;
1399  const int m_nmodes2;
1400  const bool m_colldir0;
1401  const bool m_colldir1;
1402  const bool m_colldir2;
1403  Array<OneD, const NekDouble> m_jacWStdW;
1404  Array<OneD, const NekDouble> m_base0;
1405  Array<OneD, const NekDouble> m_base1;
1406  Array<OneD, const NekDouble> m_base2;
1407  Array<OneD, const NekDouble> m_derbase0;
1408  Array<OneD, const NekDouble> m_derbase1;
1409  Array<OneD, const NekDouble> m_derbase2;
1410  Array<TwoD, const NekDouble> m_derivFac;
1411 
1412  private:
1413  IProductWRTDerivBase_SumFac_Hex(
1414  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1415  CoalescedGeomDataSharedPtr pGeomData,
1416  StdRegions::FactorMap factors)
1417  : Operator(pCollExp, pGeomData, factors),
1418  m_nquad0 (m_stdExp->GetNumPoints(0)),
1419  m_nquad1 (m_stdExp->GetNumPoints(1)),
1420  m_nquad2 (m_stdExp->GetNumPoints(2)),
1421  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1422  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1423  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1424  m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
1425  m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
1426  m_colldir2(m_stdExp->GetBasis(2)->Collocation()),
1427  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1428  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1429  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1430  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1431  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1432  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1433 
1434  {
1435  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1436  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1437  m_nmodes0*m_nmodes1*m_nmodes2));
1438  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1439  }
1440 };
1441 
1442 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Hex operator
1443 OperatorKey IProductWRTDerivBase_SumFac_Hex::m_type = GetOperatorFactory().
1444  RegisterCreatorFunction(
1445  OperatorKey(eHexahedron, eIProductWRTDerivBase, eSumFac, false),
1446  IProductWRTDerivBase_SumFac_Hex::create,
1447  "IProductWRTDerivBase_SumFac_Hex");
1448 
1449 
1450 /**
1451  * @brief Inner product WRT deriv base operator using sum-factorisation (Tet)
1452  */
1453 class IProductWRTDerivBase_SumFac_Tet : public Operator
1454 {
1455  public:
1456  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Tet)
1457 
1458 
1459  /**
1460  * This method calculates:
1461  *
1462  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
1463  *
1464  * which can be represented in terms of local cartesian
1465  * derivaties as:
1466  *
1467  * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1468  * d\phi/d\xi_1\, d\xi_1/dx +
1469  * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
1470  *
1471  * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1472  * d\phi/d\xi_1\, d\xi_1/dy +
1473  * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
1474  *
1475  * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
1476  * d\phi/d\xi_1\, d\xi_1/dz +
1477  * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
1478  *
1479  * where we note that
1480  *
1481  * \f[ d\phi/d\xi_0 = d\phi/d\eta_0 4/((1-\eta_1)(1-\eta_2)) \f]
1482  *
1483  * \f[ d\phi/d\xi_1 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
1484  * + d\phi/d\eta_1 2/(1-\eta_2) \f]
1485  *
1486  * \f[ d\phi/d\xi_2 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
1487  * + d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
1488  *
1489  * and so the full inner products are
1490  *
1491  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
1492  *
1493  * \f[ (d\phi/d\eta_0, fac0 (tmp0 + fac1(tmp1 + tmp2)))
1494  * + (d\phi/d\eta_1, fac2 (tmp1 + fac3 tmp2))
1495  * + (d\phi/d\eta_2, tmp2) \f]
1496  *
1497  * where
1498  *
1499  * \f[ \begin{array}{lcl}
1500  * tmp0 &=& (d\xi_0/dx in[0] + d\xi_0/dy in[1] + d\xi_0/dz in[2]) \\
1501  * tmp1 &=& (d\xi_1/dx in[0] + d\xi_1/dy in[1] + d\xi_1/dz in[2]) \\
1502  * tmp2 &=& (d\xi_2/dx in[0] + d\xi_2/dy in[1] + d\xi_2/dz in[2])
1503  * \end{array} \f]
1504  *
1505  * \f[ \begin{array}{lcl}
1506  * fac0 &= & 4/((1-\eta_1)(1-\eta_2)) \\
1507  * fac1 &= & (1+\eta_0)/2 \\
1508  * fac2 &= & 2/(1-\eta_2) \\
1509  * fac3 &= & (1+\eta_1)/2 \end{array} \f]
1510  *
1511  */
1512  void operator()(
1513  const Array<OneD, const NekDouble> &entry0,
1514  Array<OneD, NekDouble> &entry1,
1515  Array<OneD, NekDouble> &entry2,
1516  Array<OneD, NekDouble> &entry3,
1517  Array<OneD, NekDouble> &wsp) final
1518  {
1519 
1520  unsigned int nPhys = m_stdExp->GetTotPoints();
1521  unsigned int ntot = m_numElmt*nPhys;
1522  unsigned int nmodes = m_stdExp->GetNcoeffs();
1523  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1524  Array<OneD, Array<OneD, const NekDouble> > in(3);
1525  Array<OneD, NekDouble> output, wsp1;
1526  Array<OneD, Array<OneD, NekDouble> > tmp(3);
1527 
1528  in[0] = entry0; in[1] = entry1;
1529  in[2] = entry2;
1530 
1531  output = entry3;
1532 
1533  for(int i = 0; i < 3; ++i)
1534  {
1535  tmp[i] = wsp + i*nmax;
1536  }
1537 
1538  if(m_isDeformed)
1539  {
1540  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1541  for(int i = 0; i < 3; ++i)
1542  {
1543  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
1544  tmp[i],1);
1545  for(int j = 1; j < 3; ++j)
1546  {
1547  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
1548  in[j],1, tmp[i], 1, tmp[i],1);
1549  }
1550  }
1551  }
1552  else
1553  {
1554  Array<OneD,NekDouble> t;
1555  for(int e = 0; e < m_numElmt; ++e)
1556  {
1557  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1558  for(int i = 0; i < 3; ++i)
1559  {
1560  Vmath::Smul (m_nqe,m_derivFac[i][e],
1561  in[0] + e*m_nqe,1,
1562  t = tmp[i] + e*m_nqe,1);
1563  for(int j = 1; j < 3; ++j)
1564  {
1565  Vmath::Svtvp (m_nqe,m_derivFac[i+3*j][e],
1566  in[j] + e*m_nqe,1,
1567  tmp[i] + e*m_nqe, 1,
1568  t = tmp[i] + e*m_nqe,1);
1569  }
1570  }
1571  }
1572  }
1573 
1574  wsp1 = wsp + 3*nmax;
1575 
1576  // Sort into eta factors
1577  for (int i = 0; i < m_numElmt; ++i)
1578  {
1579  // add tmp[1] + tmp[2]
1580  Vmath::Vadd(nPhys, tmp[1].get() + i*nPhys, 1,
1581  tmp[2].get() + i*nPhys, 1,
1582  wsp1.get(), 1);
1583 
1584  // scale wsp1 by fac1 and add to tmp0
1585  Vmath::Vvtvp(nPhys,&m_fac1[0],1,wsp1.get(),1,
1586  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1587 
1588  // scale tmp[0] by fac0
1589  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
1590  tmp[0].get()+i*nPhys,1);
1591 
1592  // scale tmp[2] by fac3 and add to tmp1
1593  Vmath::Vvtvp(nPhys,&m_fac3[0],1,tmp[2].get()+i*nPhys,1,
1594  tmp[1].get()+i*nPhys,1,tmp[1].get()+i*nPhys,1);
1595 
1596  // scale tmp[1] by fac2
1597  Vmath::Vmul(nPhys,&m_fac2[0],1,tmp[1].get()+i*nPhys,1,
1598  tmp[1].get()+i*nPhys,1);
1599  }
1600 
1601 
1602  // calculate Iproduct WRT Std Deriv
1603  TetIProduct(m_sortTopEdge, m_numElmt,
1604  m_nquad0, m_nquad1, m_nquad2,
1605  m_nmodes0, m_nmodes1, m_nmodes2,
1606  m_derbase0, m_base1, m_base2,
1607  m_jacWStdW,tmp[0],output,wsp1);
1608 
1609  TetIProduct(m_sortTopEdge, m_numElmt,
1610  m_nquad0, m_nquad1, m_nquad2,
1611  m_nmodes0, m_nmodes1, m_nmodes2,
1612  m_base0, m_derbase1, m_base2,
1613  m_jacWStdW,tmp[1],tmp[0],wsp1);
1614  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1615 
1616  TetIProduct(m_sortTopEdge, m_numElmt,
1617  m_nquad0, m_nquad1, m_nquad2,
1618  m_nmodes0, m_nmodes1, m_nmodes2,
1619  m_base0, m_base1, m_derbase2,
1620  m_jacWStdW,tmp[2],tmp[0],wsp1);
1621  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1622  }
1623 
1624  void operator()(int dir,
1625  const Array<OneD, const NekDouble> &input,
1626  Array<OneD, NekDouble> &output,
1627  Array<OneD, NekDouble> &wsp) final
1628  {
1629  boost::ignore_unused(dir, input, output, wsp);
1630  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1631  }
1632 
1633  virtual void CheckFactors(StdRegions::FactorMap factors,
1634  int coll_phys_offset)
1635  {
1636  boost::ignore_unused(factors, coll_phys_offset);
1637  ASSERTL0(false, "Not valid for this operator.");
1638  }
1639 
1640  protected:
1641  const int m_nquad0;
1642  const int m_nquad1;
1643  const int m_nquad2;
1644  const int m_nmodes0;
1645  const int m_nmodes1;
1646  const int m_nmodes2;
1647  Array<OneD, const NekDouble> m_jacWStdW;
1648  Array<OneD, const NekDouble> m_base0;
1649  Array<OneD, const NekDouble> m_base1;
1650  Array<OneD, const NekDouble> m_base2;
1651  Array<OneD, const NekDouble> m_derbase0;
1652  Array<OneD, const NekDouble> m_derbase1;
1653  Array<OneD, const NekDouble> m_derbase2;
1654  Array<TwoD, const NekDouble> m_derivFac;
1655  Array<OneD, NekDouble> m_fac0;
1656  Array<OneD, NekDouble> m_fac1;
1657  Array<OneD, NekDouble> m_fac2;
1658  Array<OneD, NekDouble> m_fac3;
1659  bool m_sortTopEdge;
1660 
1661  private:
1662  IProductWRTDerivBase_SumFac_Tet(
1663  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1664  CoalescedGeomDataSharedPtr pGeomData,
1665  StdRegions::FactorMap factors)
1666  : Operator(pCollExp, pGeomData, factors),
1667  m_nquad0 (m_stdExp->GetNumPoints(0)),
1668  m_nquad1 (m_stdExp->GetNumPoints(1)),
1669  m_nquad2 (m_stdExp->GetNumPoints(2)),
1670  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1671  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1672  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1673  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1674  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1675  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1676  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1677  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1678  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1679 
1680  {
1681  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1682  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1683  m_nmodes0*m_nmodes1*m_nmodes2));
1684  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1685 
1686 
1687  const Array<OneD, const NekDouble>& z0
1688  = m_stdExp->GetBasis(0)->GetZ();
1689  const Array<OneD, const NekDouble>& z1
1690  = m_stdExp->GetBasis(1)->GetZ();
1691  const Array<OneD, const NekDouble>& z2
1692  = m_stdExp->GetBasis(2)->GetZ();
1693 
1694  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1695  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1696  m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1697  m_fac3 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1698  // calculate 2.0/((1-eta_1)(1-eta_2))
1699  for (int i = 0; i < m_nquad0; ++i)
1700  {
1701  for(int j = 0; j < m_nquad1; ++j)
1702  {
1703  for(int k = 0; k < m_nquad2; ++k)
1704  {
1705 
1706  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1707  = 4.0/((1-z1[j])*(1-z2[k]));
1708  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1709  = (1+z0[i])*0.5;
1710  m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1711  = 2.0/(1-z2[k]);
1712  m_fac3[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1713  = (1+z1[j])*0.5;
1714  }
1715  }
1716  }
1717 
1718  if(m_stdExp->GetBasis(0)->GetBasisType()
1719  == LibUtilities::eModified_A)
1720  {
1721  m_sortTopEdge = true;
1722  }
1723  else
1724  {
1725  m_sortTopEdge = false;
1726  }
1727 
1728  }
1729 };
1730 
1731 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Tet operator
1732 OperatorKey IProductWRTDerivBase_SumFac_Tet::m_type = GetOperatorFactory().
1733  RegisterCreatorFunction(
1734  OperatorKey(eTetrahedron, eIProductWRTDerivBase, eSumFac,false),
1735  IProductWRTDerivBase_SumFac_Tet::create,
1736  "IProductWRTDerivBase_SumFac_Tet");
1737 
1738 
1739 /**
1740  * @brief Inner product WRT deriv base operator using sum-factorisation (Prism)
1741  */
1742 class IProductWRTDerivBase_SumFac_Prism : public Operator
1743 {
1744  public:
1745  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Prism)
1746 
1747  ~IProductWRTDerivBase_SumFac_Prism() final
1748  {
1749  }
1750 
1751  /**
1752  * This method calculates:
1753  *
1754  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
1755  *
1756  * which can be represented in terms of local cartesian
1757  * derivaties as:
1758  *
1759  * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1760  * d\phi/d\xi_1\, d\xi_1/dx +
1761  * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
1762  *
1763  * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1764  * d\phi/d\xi_1\, d\xi_1/dy +
1765  * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
1766  *
1767  * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
1768  * d\phi/d\xi_1\, d\xi_1/dz +
1769  * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
1770  *
1771  * where we note that
1772  *
1773  * \f[ d\phi/d\xi_0 =
1774  * d\phi/d\eta_0 d\eta_0/d\xi_0 = d\phi/d\eta_0 2/(1-\eta_2) \f]
1775  *
1776  * \f[ d\phi/d\xi_2 =
1777  * d\phi/d\eta_0 d\eta_0/d\xi_2 + d\phi/d\eta_2 d\eta_2/d\xi_2 =
1778  * d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) + d\phi/d\eta_2 \f]
1779  *
1780  *
1781  * and so the full inner products are
1782  *
1783  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
1784  *
1785  * \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] + d\xi_0/dy in[1]
1786  * + d\xi_0/dz in[2])
1787  * + (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
1788  * + d\xi_2/dz in[2] )) + \f]
1789  *
1790  * \f[ (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1]
1791  * + d\xi_1/dz in[2])) + \f]
1792  *
1793  * \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1]
1794  * + d\xi_2/dz in[2])) \f]
1795  *
1796  */
1797  void operator()(
1798  const Array<OneD, const NekDouble> &entry0,
1799  Array<OneD, NekDouble> &entry1,
1800  Array<OneD, NekDouble> &entry2,
1801  Array<OneD, NekDouble> &entry3,
1802  Array<OneD, NekDouble> &wsp) final
1803  {
1804 
1805  unsigned int nPhys = m_stdExp->GetTotPoints();
1806  unsigned int ntot = m_numElmt*nPhys;
1807  unsigned int nmodes = m_stdExp->GetNcoeffs();
1808  unsigned int nmax = max(ntot,m_numElmt*nmodes);
1809  Array<OneD, Array<OneD, const NekDouble> > in(3);
1810  Array<OneD, NekDouble> output, wsp1;
1811  Array<OneD, Array<OneD, NekDouble> > tmp(3);
1812 
1813  in[0] = entry0; in[1] = entry1;
1814  in[2] = entry2;
1815 
1816  output = entry3;
1817 
1818  for(int i = 0; i < 3; ++i)
1819  {
1820  tmp[i] = wsp + i*nmax;
1821  }
1822 
1823  if(m_isDeformed)
1824  {
1825  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1826  for(int i = 0; i < 3; ++i)
1827  {
1828  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
1829  tmp[i],1);
1830  for(int j = 1; j < 3; ++j)
1831  {
1832  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
1833  in[j],1, tmp[i], 1, tmp[i],1);
1834  }
1835  }
1836  }
1837  else
1838  {
1839  Array<OneD,NekDouble> t;
1840  for(int e = 0; e < m_numElmt; ++e)
1841  {
1842  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1843  for(int i = 0; i < 3; ++i)
1844  {
1845  Vmath::Smul (m_nqe,m_derivFac[i][e],
1846  in[0] + e*m_nqe,1,
1847  t = tmp[i] + e*m_nqe,1);
1848  for(int j = 1; j < 3; ++j)
1849  {
1850  Vmath::Svtvp (m_nqe,m_derivFac[i+3*j][e],
1851  in[j] + e*m_nqe,1,
1852  tmp[i] + e*m_nqe, 1,
1853  t = tmp[i] + e*m_nqe,1);
1854  }
1855  }
1856  }
1857  }
1858 
1859  wsp1 = wsp + 3*nmax;
1860 
1861  // Sort into eta factors
1862  for (int i = 0; i < m_numElmt; ++i)
1863  {
1864  // scale tmp[0] by fac0
1865  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
1866  tmp[0].get()+i*nPhys,1);
1867 
1868  // scale tmp[2] by fac1 and add to tmp0
1869  Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[2].get()+i*nPhys,1,
1870  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
1871  }
1872 
1873  // calculate Iproduct WRT Std Deriv
1874  PrismIProduct(m_sortTopVertex, m_numElmt,
1875  m_nquad0, m_nquad1, m_nquad2,
1876  m_nmodes0, m_nmodes1, m_nmodes2,
1877  m_derbase0, m_base1, m_base2,
1878  m_jacWStdW,tmp[0],output,wsp1);
1879 
1880  PrismIProduct(m_sortTopVertex, m_numElmt,
1881  m_nquad0, m_nquad1, m_nquad2,
1882  m_nmodes0, m_nmodes1, m_nmodes2,
1883  m_base0, m_derbase1, m_base2,
1884  m_jacWStdW,tmp[1],tmp[0],wsp1);
1885  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1886 
1887  PrismIProduct(m_sortTopVertex, m_numElmt,
1888  m_nquad0, m_nquad1, m_nquad2,
1889  m_nmodes0, m_nmodes1, m_nmodes2,
1890  m_base0, m_base1, m_derbase2,
1891  m_jacWStdW,tmp[2],tmp[0],wsp1);
1892  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
1893  }
1894 
1895  void operator()(int dir,
1896  const Array<OneD, const NekDouble> &input,
1897  Array<OneD, NekDouble> &output,
1898  Array<OneD, NekDouble> &wsp) final
1899  {
1900  boost::ignore_unused(dir, input, output, wsp);
1901  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1902  }
1903 
1904  virtual void CheckFactors(StdRegions::FactorMap factors,
1905  int coll_phys_offset)
1906  {
1907  boost::ignore_unused(factors, coll_phys_offset);
1908  ASSERTL0(false, "Not valid for this operator.");
1909  }
1910 
1911 
1912  protected:
1913  const int m_nquad0;
1914  const int m_nquad1;
1915  const int m_nquad2;
1916  const int m_nmodes0;
1917  const int m_nmodes1;
1918  const int m_nmodes2;
1919  Array<OneD, const NekDouble> m_jacWStdW;
1920  Array<OneD, const NekDouble> m_base0;
1921  Array<OneD, const NekDouble> m_base1;
1922  Array<OneD, const NekDouble> m_base2;
1923  Array<OneD, const NekDouble> m_derbase0;
1924  Array<OneD, const NekDouble> m_derbase1;
1925  Array<OneD, const NekDouble> m_derbase2;
1926  Array<TwoD, const NekDouble> m_derivFac;
1927  Array<OneD, NekDouble> m_fac0;
1928  Array<OneD, NekDouble> m_fac1;
1929  bool m_sortTopVertex;
1930 
1931  private:
1932  IProductWRTDerivBase_SumFac_Prism(
1933  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1934  CoalescedGeomDataSharedPtr pGeomData,
1935  StdRegions::FactorMap factors)
1936  : Operator(pCollExp, pGeomData, factors),
1937  m_nquad0 (m_stdExp->GetNumPoints(0)),
1938  m_nquad1 (m_stdExp->GetNumPoints(1)),
1939  m_nquad2 (m_stdExp->GetNumPoints(2)),
1940  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
1941  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
1942  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
1943  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
1944  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
1945  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
1946  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1947  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1948  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1949 
1950  {
1951  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1952  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
1953  m_nmodes0*m_nmodes1*m_nmodes2));
1954  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1955 
1956  if(m_stdExp->GetBasis(0)->GetBasisType()
1957  == LibUtilities::eModified_A)
1958  {
1959  m_sortTopVertex = true;
1960  }
1961  else
1962  {
1963  m_sortTopVertex = false;
1964  }
1965 
1966  const Array<OneD, const NekDouble>& z0
1967  = m_stdExp->GetBasis(0)->GetZ();
1968  const Array<OneD, const NekDouble>& z2
1969  = m_stdExp->GetBasis(2)->GetZ();
1970 
1971  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1972  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1973 
1974  for (int i = 0; i < m_nquad0; ++i)
1975  {
1976  for(int j = 0; j < m_nquad1; ++j)
1977  {
1978  for(int k = 0; k < m_nquad2; ++k)
1979  {
1980  // set up geometric factor: 2/(1-z1)
1981  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1982  = 2.0/(1-z2[k]);
1983  // set up geometric factor: (1+z0)/(1-z1)
1984  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1985  = (1+z0[i])/(1-z2[k]);
1986 
1987  }
1988  }
1989  }
1990  }
1991 };
1992 
1993 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Prism operator
1994 OperatorKey IProductWRTDerivBase_SumFac_Prism::m_type = GetOperatorFactory().
1995  RegisterCreatorFunction(
1996  OperatorKey(ePrism, eIProductWRTDerivBase, eSumFac, false),
1997  IProductWRTDerivBase_SumFac_Prism::create,
1998  "IProductWRTDerivBase_SumFac_Prism");
1999 
2000 
2001 
2002 /**
2003  * @brief Inner product WRT deriv base operator using sum-factorisation (Pyr)
2004  */
2005 class IProductWRTDerivBase_SumFac_Pyr : public Operator
2006 {
2007  public:
2008  OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Pyr)
2009 
2010  ~IProductWRTDerivBase_SumFac_Pyr() final
2011  {
2012  }
2013 
2014 
2015  /**
2016  * This method calculates:
2017  *
2018  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
2019  *
2020  * which can be represented in terms of local cartesian
2021  * derivaties as:
2022  *
2023  * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
2024  * d\phi/d\xi_1\, d\xi_1/dx +
2025  * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
2026  *
2027  * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
2028  * d\phi/d\xi_1\, d\xi_1/dy +
2029  * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
2030  *
2031  * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
2032  * d\phi/d\xi_1\, d\xi_1/dz +
2033  * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
2034  *
2035  * where we note that
2036  *
2037  * \f[ d\phi/d\xi_0 =
2038  * d\phi/d\eta_0\, d\eta_0/d\xi_0 =
2039  * d\phi/d\eta_0\, 2/(1-\eta_2). \f]
2040  *
2041  * \f[ d\phi/d\xi_1 =
2042  * d\phi/d\eta_1\, d\eta_1/d\xi_1 =
2043  * d\phi/d\eta_1\, 2/(1-\eta_2) \f]
2044  *
2045  * \f[ d\phi/d\xi_2 =
2046  * d\phi/d\eta_0\, d\eta_0/d\xi_2 +
2047  * d\phi/d\eta_1\, d\eta_1/d\xi_2 +
2048  * d\phi/d\eta_2\, d\eta_2/d\xi_2 =
2049  * d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) +
2050  * d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
2051  *
2052  * and so the full inner products are
2053  *
2054  * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
2055  *
2056  * \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] +
2057  * d\xi_0/dy in[1] +
2058  * (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
2059  * + d\xi_2/dz in[2] )) + \f]
2060  * \f[ (d\phi/d\eta_1, ((2/(1-\eta_2) (d\xi_1/dx in[0] +
2061  * d\xi_0/dy in[1] + d\xi_0/dz in[2]) +
2062  * (1-\eta_1)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
2063  * d\xi_2/dz in[2] )) \f]
2064  *
2065  * \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
2066  * d\xi_2/dz in[2])) \f]
2067  */
2068  void operator()(
2069  const Array<OneD, const NekDouble> &entry0,
2070  Array<OneD, NekDouble> &entry1,
2071  Array<OneD, NekDouble> &entry2,
2072  Array<OneD, NekDouble> &entry3,
2073  Array<OneD, NekDouble> &wsp) final
2074  {
2075  unsigned int nPhys = m_stdExp->GetTotPoints();
2076  unsigned int ntot = m_numElmt*nPhys;
2077  unsigned int nmodes = m_stdExp->GetNcoeffs();
2078  unsigned int nmax = max(ntot,m_numElmt*nmodes);
2079  Array<OneD, Array<OneD, const NekDouble> > in(3);
2080  Array<OneD, NekDouble> output, wsp1;
2081  Array<OneD, Array<OneD, NekDouble> > tmp(3);
2082 
2083  in[0] = entry0; in[1] = entry1;
2084  in[2] = entry2;
2085 
2086  output = entry3;
2087 
2088  for(int i = 0; i < 3; ++i)
2089  {
2090  tmp[i] = wsp + i*nmax;
2091  }
2092 
2093  if(m_isDeformed)
2094  {
2095  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
2096  for(int i = 0; i < 3; ++i)
2097  {
2098  Vmath::Vmul (ntot,m_derivFac[i],1, in[0],1,
2099  tmp[i],1);
2100  for(int j = 1; j < 3; ++j)
2101  {
2102  Vmath::Vvtvp (ntot,m_derivFac[i+3*j],1,
2103  in[j],1, tmp[i], 1, tmp[i],1);
2104  }
2105  }
2106  }
2107  else
2108  {
2109  Array<OneD,NekDouble> t;
2110  for(int e = 0; e < m_numElmt; ++e)
2111  {
2112  // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
2113  for(int i = 0; i < 3; ++i)
2114  {
2115  Vmath::Smul (m_nqe,m_derivFac[i][e],
2116  in[0] + e*m_nqe,1,
2117  t = tmp[i] + e*m_nqe,1);
2118  for(int j = 1; j < 3; ++j)
2119  {
2120  Vmath::Svtvp (m_nqe,m_derivFac[i+3*j][e],
2121  in[j] + e*m_nqe,1,
2122  tmp[i] + e*m_nqe, 1,
2123  t = tmp[i] + e*m_nqe,1);
2124  }
2125  }
2126  }
2127  }
2128 
2129  wsp1 = wsp + 3*nmax;
2130 
2131  // Sort into eta factors
2132  for (int i = 0; i < m_numElmt; ++i)
2133  {
2134  // scale tmp[0] by fac0
2135  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[0].get()+i*nPhys,1,
2136  tmp[0].get()+i*nPhys,1);
2137 
2138  // scale tmp[2] by fac1 and add to tmp0
2139  Vmath::Vvtvp(nPhys,&m_fac1[0],1,tmp[2].get()+i*nPhys,1,
2140  tmp[0].get()+i*nPhys,1,tmp[0].get()+i*nPhys,1);
2141 
2142  // scale tmp[1] by fac0
2143  Vmath::Vmul(nPhys,&m_fac0[0],1,tmp[1].get()+i*nPhys,1,
2144  tmp[1].get()+i*nPhys,1);
2145 
2146  // scale tmp[2] by fac2 and add to tmp1
2147  Vmath::Vvtvp(nPhys,&m_fac2[0],1,tmp[2].get()+i*nPhys,1,
2148  tmp[1].get()+i*nPhys,1,tmp[1].get()+i*nPhys,1);
2149  }
2150 
2151  // calculate Iproduct WRT Std Deriv
2152  PyrIProduct(m_sortTopVertex, m_numElmt,
2153  m_nquad0, m_nquad1, m_nquad2,
2154  m_nmodes0, m_nmodes1, m_nmodes2,
2155  m_derbase0, m_base1, m_base2,
2156  m_jacWStdW,tmp[0],output,wsp1);
2157 
2158  PyrIProduct(m_sortTopVertex, m_numElmt,
2159  m_nquad0, m_nquad1, m_nquad2,
2160  m_nmodes0, m_nmodes1, m_nmodes2,
2161  m_base0, m_derbase1, m_base2,
2162  m_jacWStdW,tmp[1],tmp[0],wsp1);
2163  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
2164 
2165  PyrIProduct(m_sortTopVertex, m_numElmt,
2166  m_nquad0, m_nquad1, m_nquad2,
2167  m_nmodes0, m_nmodes1, m_nmodes2,
2168  m_base0, m_base1, m_derbase2,
2169  m_jacWStdW,tmp[2],tmp[0],wsp1);
2170  Vmath::Vadd(m_numElmt*nmodes,tmp[0],1,output,1,output,1);
2171  }
2172 
2173  void operator()(int dir,
2174  const Array<OneD, const NekDouble> &input,
2175  Array<OneD, NekDouble> &output,
2176  Array<OneD, NekDouble> &wsp) final
2177  {
2178  boost::ignore_unused(dir, input, output, wsp);
2179  NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
2180  }
2181 
2182  virtual void CheckFactors(StdRegions::FactorMap factors,
2183  int coll_phys_offset)
2184  {
2185  boost::ignore_unused(factors, coll_phys_offset);
2186  ASSERTL0(false, "Not valid for this operator.");
2187  }
2188 
2189 
2190  protected:
2191  const int m_nquad0;
2192  const int m_nquad1;
2193  const int m_nquad2;
2194  const int m_nmodes0;
2195  const int m_nmodes1;
2196  const int m_nmodes2;
2197  Array<OneD, const NekDouble> m_jacWStdW;
2198  Array<OneD, const NekDouble> m_base0;
2199  Array<OneD, const NekDouble> m_base1;
2200  Array<OneD, const NekDouble> m_base2;
2201  Array<OneD, const NekDouble> m_derbase0;
2202  Array<OneD, const NekDouble> m_derbase1;
2203  Array<OneD, const NekDouble> m_derbase2;
2204  Array<TwoD, const NekDouble> m_derivFac;
2205  Array<OneD, NekDouble> m_fac0;
2206  Array<OneD, NekDouble> m_fac1;
2207  Array<OneD, NekDouble> m_fac2;
2208  bool m_sortTopVertex;
2209 
2210  private:
2211  IProductWRTDerivBase_SumFac_Pyr(
2212  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
2213  CoalescedGeomDataSharedPtr pGeomData,
2214  StdRegions::FactorMap factors)
2215  : Operator(pCollExp, pGeomData, factors),
2216  m_nquad0 (m_stdExp->GetNumPoints(0)),
2217  m_nquad1 (m_stdExp->GetNumPoints(1)),
2218  m_nquad2 (m_stdExp->GetNumPoints(2)),
2219  m_nmodes0 (m_stdExp->GetBasisNumModes(0)),
2220  m_nmodes1 (m_stdExp->GetBasisNumModes(1)),
2221  m_nmodes2 (m_stdExp->GetBasisNumModes(2)),
2222  m_base0 (m_stdExp->GetBasis(0)->GetBdata()),
2223  m_base1 (m_stdExp->GetBasis(1)->GetBdata()),
2224  m_base2 (m_stdExp->GetBasis(2)->GetBdata()),
2225  m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
2226  m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
2227  m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
2228 
2229  {
2230  m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
2231  m_wspSize = 6 * m_numElmt * (max(m_nquad0*m_nquad1*m_nquad2,
2232  m_nmodes0*m_nmodes1*m_nmodes2));
2233  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
2234 
2235  if(m_stdExp->GetBasis(0)->GetBasisType()
2236  == LibUtilities::eModified_A)
2237  {
2238  m_sortTopVertex = true;
2239  }
2240  else
2241  {
2242  m_sortTopVertex = false;
2243  }
2244 
2245  const Array<OneD, const NekDouble>& z0
2246  = m_stdExp->GetBasis(0)->GetZ();
2247  const Array<OneD, const NekDouble>& z1
2248  = m_stdExp->GetBasis(1)->GetZ();
2249  const Array<OneD, const NekDouble>& z2
2250  = m_stdExp->GetBasis(2)->GetZ();
2251 
2252  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
2253  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
2254  m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
2255 
2256  for (int i = 0; i < m_nquad0; ++i)
2257  {
2258  for(int j = 0; j < m_nquad1; ++j)
2259  {
2260  for(int k = 0; k < m_nquad2; ++k)
2261  {
2262  // set up geometric factor: 2/(1-z2)
2263  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
2264  = 2.0/(1-z2[k]);
2265  // set up geometric factor: (1+z0)/(1-z2)
2266  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
2267  = (1+z0[i])/(1-z2[k]);
2268  // set up geometric factor: (1+z1)/(1-z2)
2269  m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
2270  = (1+z1[j])/(1-z2[k]);
2271 
2272  }
2273  }
2274  }
2275  }
2276 };
2277 
2278 /// Factory initialisation for the IProductWRTDerivBase_SumFac_Pyr operator
2279 OperatorKey IProductWRTDerivBase_SumFac_Pyr::m_type = GetOperatorFactory().
2280  RegisterCreatorFunction(
2281  OperatorKey(ePyramid, eIProductWRTDerivBase, eSumFac, false),
2282  IProductWRTDerivBase_SumFac_Pyr::create,
2283  "IProductWRTDerivBase_SumFac_Pyr");
2284 
2285 
2286 }
2287 }
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
OPERATOR_CREATE
#define OPERATOR_CREATE(cname)
Definition: Operator.h:45
Nektar::Collections::IProductWRTDerivBase_IterPerExp
Inner product WRT deriv base operator using element-wise operation.
Definition: IProductWRTDerivBase.cpp:453
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:580
Nektar::Collections::IProductWRTDerivBase_IterPerExp::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:461
Nektar::Collections::IProductWRTDerivBase_IterPerExp::IProductWRTDerivBase_IterPerExp
IProductWRTDerivBase_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:586
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:582
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:583
Nektar::Collections::IProductWRTDerivBase_IterPerExp::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:563
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_jac
Array< OneD, const NekDouble > m_jac
Definition: IProductWRTDerivBase.cpp:581
Nektar::Collections::IProductWRTDerivBase_IterPerExp::~IProductWRTDerivBase_IterPerExp
~IProductWRTDerivBase_IterPerExp() final
Definition: IProductWRTDerivBase.cpp:457
Nektar::Collections::IProductWRTDerivBase_IterPerExp::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:572
Nektar::Collections::IProductWRTDerivBase_MatrixFree
Inner product operator using operator using matrix free operators.
Definition: IProductWRTDerivBase.cpp:282
Nektar::Collections::IProductWRTDerivBase_MatrixFree::~IProductWRTDerivBase_MatrixFree
~IProductWRTDerivBase_MatrixFree() final
Definition: IProductWRTDerivBase.cpp:286
Nektar::Collections::IProductWRTDerivBase_MatrixFree::IProductWRTDerivBase_MatrixFree
IProductWRTDerivBase_MatrixFree(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:381
Nektar::Collections::IProductWRTDerivBase_MatrixFree::m_oper
std::shared_ptr< MatrixFree::IProductWRTDerivBase > m_oper
Definition: IProductWRTDerivBase.cpp:379
Nektar::Collections::IProductWRTDerivBase_MatrixFree::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:362
Nektar::Collections::IProductWRTDerivBase_MatrixFree::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:371
Nektar::Collections::IProductWRTDerivBase_MatrixFree::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:290
Nektar::Collections::IProductWRTDerivBase_NoCollection
Inner product WRT deriv base operator using LocalRegions implementation.
Definition: IProductWRTDerivBase.cpp:654
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_expList
vector< StdRegions::StdExpansionSharedPtr > m_expList
Definition: IProductWRTDerivBase.cpp:720
Nektar::Collections::IProductWRTDerivBase_NoCollection::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:701
Nektar::Collections::IProductWRTDerivBase_NoCollection::IProductWRTDerivBase_NoCollection
IProductWRTDerivBase_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:723
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:719
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:718
Nektar::Collections::IProductWRTDerivBase_NoCollection::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:662
Nektar::Collections::IProductWRTDerivBase_NoCollection::~IProductWRTDerivBase_NoCollection
~IProductWRTDerivBase_NoCollection() final
Definition: IProductWRTDerivBase.cpp:658
Nektar::Collections::IProductWRTDerivBase_NoCollection::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:710
Nektar::Collections::IProductWRTDerivBase_StdMat
Inner product WRT deriv base operator using standard matrix approach.
Definition: IProductWRTDerivBase.cpp:61
Nektar::Collections::IProductWRTDerivBase_StdMat::m_jac
Array< OneD, const NekDouble > m_jac
Definition: IProductWRTDerivBase.cpp:189
Nektar::Collections::IProductWRTDerivBase_StdMat::~IProductWRTDerivBase_StdMat
~IProductWRTDerivBase_StdMat() final
Definition: IProductWRTDerivBase.cpp:65
Nektar::Collections::IProductWRTDerivBase_StdMat::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:188
Nektar::Collections::IProductWRTDerivBase_StdMat::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:69
Nektar::Collections::IProductWRTDerivBase_StdMat::m_iProdWRTStdDBase
Array< OneD, DNekMatSharedPtr > m_iProdWRTStdDBase
Definition: IProductWRTDerivBase.cpp:187
Nektar::Collections::IProductWRTDerivBase_StdMat::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:170
Nektar::Collections::IProductWRTDerivBase_StdMat::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:191
Nektar::Collections::IProductWRTDerivBase_StdMat::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:190
Nektar::Collections::IProductWRTDerivBase_StdMat::IProductWRTDerivBase_StdMat
IProductWRTDerivBase_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:194
Nektar::Collections::IProductWRTDerivBase_StdMat::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:179
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex
Inner product WRT deriv base operator using sum-factorisation (Hex)
Definition: IProductWRTDerivBase.cpp:1283
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:1291
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:1400
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1408
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1396
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1377
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:1401
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1403
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1406
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::IProductWRTDerivBase_SumFac_Hex
IProductWRTDerivBase_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1413
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1395
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1410
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1409
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:1386
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir2
const bool m_colldir2
Definition: IProductWRTDerivBase.cpp:1402
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1397
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1407
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1394
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1405
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1404
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::~IProductWRTDerivBase_SumFac_Hex
~IProductWRTDerivBase_SumFac_Hex() final
Definition: IProductWRTDerivBase.cpp:1287
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1399
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1398
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism
Inner product WRT deriv base operator using sum-factorisation (Prism)
Definition: IProductWRTDerivBase.cpp:1743
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1917
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1922
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1925
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1923
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:1904
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1914
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1919
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::~IProductWRTDerivBase_SumFac_Prism
~IProductWRTDerivBase_SumFac_Prism() final
Definition: IProductWRTDerivBase.cpp:1747
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1916
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::IProductWRTDerivBase_SumFac_Prism
IProductWRTDerivBase_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1932
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1921
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1915
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1797
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1895
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1927
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1924
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1918
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1928
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1920
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1926
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:1929
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1913
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr
Inner product WRT deriv base operator using sum-factorisation (Pyr)
Definition: IProductWRTDerivBase.cpp:2006
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:2193
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:2173
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:2204
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:2200
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:2199
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:2192
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:2208
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:2068
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:2205
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:2202
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:2182
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::~IProductWRTDerivBase_SumFac_Pyr
~IProductWRTDerivBase_SumFac_Pyr() final
Definition: IProductWRTDerivBase.cpp:2010
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:2196
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:2206
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:2197
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:2201
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:2194
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:2198
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:2195
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::IProductWRTDerivBase_SumFac_Pyr
IProductWRTDerivBase_SumFac_Pyr(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:2211
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: IProductWRTDerivBase.cpp:2207
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:2191
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:2203
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad
Inner product WRT deriv base operator using sum-factorisation (Quad)
Definition: IProductWRTDerivBase.cpp:902
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1003
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:984
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1004
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:910
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1008
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1010
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:1007
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1013
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1012
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:1006
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:993
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:1005
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1001
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1009
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::IProductWRTDerivBase_SumFac_Quad
IProductWRTDerivBase_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1016
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1002
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1011
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::~IProductWRTDerivBase_SumFac_Quad
~IProductWRTDerivBase_SumFac_Quad() final
Definition: IProductWRTDerivBase.cpp:906
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg
Inner product WRT deriv base operator using sum-factorisation (Segment)
Definition: IProductWRTDerivBase.cpp:785
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::~IProductWRTDerivBase_SumFac_Seg
~IProductWRTDerivBase_SumFac_Seg() final
Definition: IProductWRTDerivBase.cpp:789
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:849
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:858
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:868
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:871
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:869
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:867
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:866
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::IProductWRTDerivBase_SumFac_Seg
IProductWRTDerivBase_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:874
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Perform operation.
Definition: IProductWRTDerivBase.cpp:793
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:870
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet
Inner product WRT deriv base operator using sum-factorisation (Tet)
Definition: IProductWRTDerivBase.cpp:1454
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1624
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1641
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: IProductWRTDerivBase.cpp:1657
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1656
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1650
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::IProductWRTDerivBase_SumFac_Tet
IProductWRTDerivBase_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1662
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1644
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_sortTopEdge
bool m_sortTopEdge
Definition: IProductWRTDerivBase.cpp:1659
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1648
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:1633
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1643
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac3
Array< OneD, NekDouble > m_fac3
Definition: IProductWRTDerivBase.cpp:1658
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1655
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1647
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1654
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1652
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1646
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1649
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1651
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1645
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1642
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1653
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri
Inner product WRT deriv base operator using sum-factorisation (Tri)
Definition: IProductWRTDerivBase.cpp:1054
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1205
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1193
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1201
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1202
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:1198
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::~IProductWRTDerivBase_SumFac_Tri
~IProductWRTDerivBase_SumFac_Tri() final
Definition: IProductWRTDerivBase.cpp:1058
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:1196
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1200
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1195
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1206
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1175
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:1197
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1199
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1194
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1203
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::operator()
void operator()(const Array< OneD, const NekDouble > &entry0, Array< OneD, NekDouble > &entry1, Array< OneD, NekDouble > &entry2, Array< OneD, NekDouble > &entry3, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1093
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1192
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::IProductWRTDerivBase_SumFac_Tri
IProductWRTDerivBase_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1210
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::CheckFactors
virtual void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset)
Check the validity of the supplied factor map.
Definition: IProductWRTDerivBase.cpp:1184
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:1207
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1204
Nektar::Collections::Operator
Base class for operators on a collection of elements.
Definition: Operator.h:115
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:200
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:84
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:394
Nektar::Collections::TetIProduct
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:522
Nektar::Collections::PyrIProduct
void PyrIProduct(bool sortTopVert, 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:427
Nektar::Collections::TriIProduct
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:134
Nektar::Collections::QuadIProduct
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:48
Nektar::Collections::OperatorKey
std::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:181
Nektar::Collections::CoalescedGeomDataSharedPtr
std::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
Definition: CoalescedGeomData.h:94
Nektar::Collections::GetOperatorFactory
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:121
Nektar::Collections::PrismIProduct
void PrismIProduct(bool sortTopVert, 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:355
Nektar::Collections::HexIProduct
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:178
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::StdRegions::FactorMap
ConstFactorMap FactorMap
Definition: StdRegions.hpp:318
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
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:192
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:565
Vmath::Vvtvp
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:513
Vmath::Vadd
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:322
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:225
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:436
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1199
Generated on Wed Nov 24 2021 10:33:27 for Nektar++ by doxygen 1.9.1