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