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