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