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 <Collections/Operator.h>
38 #include <Collections/Collection.h>
39 
40 using namespace std;
41 
42 namespace Nektar {
43 namespace Collections {
44 
45 using LibUtilities::eSegment;
46 using LibUtilities::eQuadrilateral;
47 using LibUtilities::eTriangle;
48 using LibUtilities::eHexahedron;
49 using LibUtilities::eTetrahedron;
50 using LibUtilities::ePrism;
51 using LibUtilities::ePyramid;
52 
53 /**
54  * @brief Phys deriv operator using standard matrix approach
55  */
56 class PhysDeriv_StdMat : public Operator
57 {
58  public:
59  OPERATOR_CREATE(PhysDeriv_StdMat)
60 
61  virtual ~PhysDeriv_StdMat()
62  {
63  }
64 
65  virtual void operator()(
66  const Array<OneD, const NekDouble> &input,
67  Array<OneD, NekDouble> &output0,
68  Array<OneD, NekDouble> &output1,
69  Array<OneD, NekDouble> &output2,
70  Array<OneD, NekDouble> &wsp)
71  {
72  int nPhys = m_stdExp->GetTotPoints();
73  int ntot = m_numElmt*nPhys;
74  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
75  Array<OneD, Array<OneD, NekDouble> > Diff(3);
76  Array<OneD, Array<OneD, NekDouble> > out(3);
77  out[0] = output0; out[1] = output1; out[2] = output2;
78 
79  for(int i = 0; i < m_dim; ++i)
80  {
81  Diff[i] = wsp + i*ntot;
82  }
83 
84  // calculate local derivatives
85  for(int i = 0; i < m_dim; ++i)
86  {
87  Blas::Dgemm('N', 'N', m_derivMat[i]->GetRows(), m_numElmt,
88  m_derivMat[i]->GetColumns(), 1.0,
89  m_derivMat[i]->GetRawPtr(),
90  m_derivMat[i]->GetRows(), input.get(), nPhys,
91  0.0, &Diff[i][0],nPhys);
92  }
93 
94  // calculate full derivative
95  for(int i = 0; i < m_coordim; ++i)
96  {
97  Vmath::Zero(ntot,out[i],1);
98  for(int j = 0; j < m_dim; ++j)
99  {
100  Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
101  Diff[j], 1,
102  out[i], 1,
103  out[i], 1);
104  }
105  }
106  }
107 
108  virtual void operator()(
109  int dir,
110  const Array<OneD, const NekDouble> &input,
111  Array<OneD, NekDouble> &output,
112  Array<OneD, NekDouble> &wsp)
113  {
114  int nPhys = m_stdExp->GetTotPoints();
115  int ntot = m_numElmt*nPhys;
116  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
117  Array<OneD, Array<OneD, NekDouble> > Diff(3);
118 
119  for(int i = 0; i < m_dim; ++i)
120  {
121  Diff[i] = wsp + i*ntot;
122  }
123 
124  // calculate local derivatives
125  for(int i = 0; i < m_dim; ++i)
126  {
127  Blas::Dgemm('N', 'N', m_derivMat[i]->GetRows(), m_numElmt,
128  m_derivMat[i]->GetColumns(), 1.0,
129  m_derivMat[i]->GetRawPtr(),
130  m_derivMat[i]->GetRows(), input.get(), nPhys,
131  0.0, &Diff[i][0],nPhys);
132  }
133 
134  // calculate full derivative
135  Vmath::Zero(ntot,output,1);
136  for(int j = 0; j < m_dim; ++j)
137  {
138  Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
139  Diff[j], 1,
140  output, 1,
141  output, 1);
142  }
143  }
144 
145  protected:
146  Array<OneD, DNekMatSharedPtr> m_derivMat;
147  Array<TwoD, const NekDouble> m_derivFac;
148  int m_dim;
149  int m_coordim;
150 
151  private:
152  PhysDeriv_StdMat(
153  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
154  CoalescedGeomDataSharedPtr pGeomData)
155  : Operator(pCollExp, pGeomData)
156  {
157  int nqtot = 1;
158  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
159  m_dim = PtsKey.size();
160  m_coordim = pCollExp[0]->GetCoordim();
161 
162  for(int i = 0; i < m_dim; ++i)
163  {
164  nqtot *= PtsKey[i].GetNumPoints();
165  }
166  // set up a PhysDeriv StdMat.
167  m_derivMat = Array<OneD, DNekMatSharedPtr>(m_dim);
168  for(int i = 0; i < m_dim; ++i)
169  {
170  Array<OneD, NekDouble> tmp(nqtot),tmp1(nqtot);
171  m_derivMat[i] = MemoryManager<DNekMat>
172  ::AllocateSharedPtr(nqtot,nqtot);
173  for(int j = 0; j < nqtot; ++j)
174  {
175  Vmath::Zero(nqtot,tmp,1);
176  tmp[j] = 1.0;
177  m_stdExp->PhysDeriv(i,tmp,tmp1);
178  Vmath::Vcopy(nqtot, &tmp1[0], 1,
179  &(m_derivMat[i]->GetPtr())[0] + j*nqtot, 1);
180  }
181  }
182  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
183  m_wspSize = 3*nqtot*m_numElmt;
184  }
185 };
186 
187 /// Factory initialisation for the PhysDeriv_StdMat operators
188 OperatorKey PhysDeriv_StdMat::m_typeArr[] =
189 {
190  GetOperatorFactory().RegisterCreatorFunction(
191  OperatorKey(eSegment, ePhysDeriv, eStdMat, false),
192  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Seg"),
193  GetOperatorFactory().RegisterCreatorFunction(
194  OperatorKey(eTriangle, ePhysDeriv, eStdMat, false),
195  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Tri"),
196  GetOperatorFactory().RegisterCreatorFunction(
197  OperatorKey(eTriangle, ePhysDeriv, eStdMat, true),
198  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_NodalTri"),
199  GetOperatorFactory().RegisterCreatorFunction(
200  OperatorKey(eQuadrilateral, ePhysDeriv, eStdMat, false),
201  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Quad"),
202  GetOperatorFactory().RegisterCreatorFunction(
203  OperatorKey(eTetrahedron, ePhysDeriv, eStdMat, false),
204  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Tet"),
205  GetOperatorFactory().RegisterCreatorFunction(
206  OperatorKey(eTetrahedron, ePhysDeriv, eStdMat, true),
207  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_NodalTet"),
208  GetOperatorFactory().RegisterCreatorFunction(
209  OperatorKey(ePyramid, ePhysDeriv, eStdMat, false),
210  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Pyr"),
211  GetOperatorFactory().RegisterCreatorFunction(
212  OperatorKey(ePrism, ePhysDeriv, eStdMat, false),
213  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Prism"),
214  GetOperatorFactory().RegisterCreatorFunction(
215  OperatorKey(ePrism, ePhysDeriv, eStdMat, true),
216  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_NodalPrism"),
217  GetOperatorFactory().RegisterCreatorFunction(
218  OperatorKey(eHexahedron, ePhysDeriv, eStdMat, false),
219  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Hex"),
220  GetOperatorFactory().RegisterCreatorFunction(
221  OperatorKey(ePyramid, ePhysDeriv, eSumFac, false),
222  PhysDeriv_StdMat::create, "PhysDeriv_SumFac_Pyr")
223 };
224 
225 
226 /**
227  * @brief Phys deriv operator using element-wise operation
228  */
229 class PhysDeriv_IterPerExp : public Operator
230 {
231  public:
232  OPERATOR_CREATE(PhysDeriv_IterPerExp)
233 
234  virtual ~PhysDeriv_IterPerExp()
235  {
236  }
237 
238  virtual void operator()(
239  const Array<OneD, const NekDouble> &input,
240  Array<OneD, NekDouble> &output0,
241  Array<OneD, NekDouble> &output1,
242  Array<OneD, NekDouble> &output2,
243  Array<OneD, NekDouble> &wsp)
244  {
245  int nPhys = m_stdExp->GetTotPoints();
246  int ntot = m_numElmt*nPhys;
247  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
248  Array<OneD, Array<OneD, NekDouble> > Diff(3);
249  Array<OneD, Array<OneD, NekDouble> > out(3);
250  out[0] = output0; out[1] = output1; out[2] = output2;
251 
252  for(int i = 0; i < m_dim; ++i)
253  {
254  Diff[i] = wsp + i*ntot;
255  }
256 
257  // calculate local derivatives
258  for (int i = 0; i < m_numElmt; ++i)
259  {
260  m_stdExp->PhysDeriv(input + i*nPhys,
261  tmp0 = Diff[0] + i*nPhys,
262  tmp1 = Diff[1] + i*nPhys,
263  tmp2 = Diff[2] + i*nPhys);
264  }
265 
266  // calculate full derivative
267  for(int i = 0; i < m_coordim; ++i)
268  {
269  Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
270  for(int j = 1; j < m_dim; ++j)
271  {
272  Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
273  Diff[j], 1,
274  out[i], 1,
275  out[i], 1);
276  }
277  }
278  }
279 
280  virtual void operator()(
281  int dir,
282  const Array<OneD, const NekDouble> &input,
283  Array<OneD, NekDouble> &output,
284  Array<OneD, NekDouble> &wsp)
285  {
286  int nPhys = m_stdExp->GetTotPoints();
287  int ntot = m_numElmt*nPhys;
288  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
289  Array<OneD, Array<OneD, NekDouble> > Diff(3);
290 
291  for(int i = 0; i < m_dim; ++i)
292  {
293  Diff[i] = wsp + i*ntot;
294  }
295 
296  // calculate local derivatives
297  for (int i = 0; i < m_numElmt; ++i)
298  {
299  m_stdExp->PhysDeriv(input + i*nPhys,
300  tmp0 = Diff[0] + i*nPhys,
301  tmp1 = Diff[1] + i*nPhys,
302  tmp2 = Diff[2] + i*nPhys);
303  }
304 
305  // calculate full derivative
306  Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
307  for(int j = 1; j < m_dim; ++j)
308  {
309  Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
310  Diff[j], 1,
311  output, 1,
312  output, 1);
313  }
314  }
315 
316  protected:
317  Array<TwoD, const NekDouble> m_derivFac;
318  int m_dim;
319  int m_coordim;
320 
321  private:
322  PhysDeriv_IterPerExp(
323  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
324  CoalescedGeomDataSharedPtr pGeomData)
325  : Operator(pCollExp, pGeomData)
326  {
327  int nqtot = 1;
328  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
329  m_dim = PtsKey.size();
330  m_coordim = pCollExp[0]->GetCoordim();
331 
332  for(int i = 0; i < m_dim; ++i)
333  {
334  nqtot *= PtsKey[i].GetNumPoints();
335  }
336  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
337  m_wspSize = 3*nqtot*m_numElmt;
338  }
339 };
340 
341 /// Factory initialisation for the PhysDeriv_IterPerExp operators
342 OperatorKey PhysDeriv_IterPerExp::m_typeArr[] =
343 {
344  GetOperatorFactory().RegisterCreatorFunction(
345  OperatorKey(eSegment, ePhysDeriv, eIterPerExp,false),
346  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Seg"),
347  GetOperatorFactory().RegisterCreatorFunction(
348  OperatorKey(eTriangle, ePhysDeriv, eIterPerExp,false),
349  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Tri"),
350  GetOperatorFactory().RegisterCreatorFunction(
351  OperatorKey(eTriangle, ePhysDeriv, eIterPerExp,true),
352  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_NodalTri"),
353  GetOperatorFactory().RegisterCreatorFunction(
354  OperatorKey(eQuadrilateral, ePhysDeriv, eIterPerExp,false),
355  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Quad"),
356  GetOperatorFactory().RegisterCreatorFunction(
357  OperatorKey(eTetrahedron, ePhysDeriv, eIterPerExp,false),
358  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Tet"),
359  GetOperatorFactory().RegisterCreatorFunction(
360  OperatorKey(eTetrahedron, ePhysDeriv, eIterPerExp,true),
361  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_NodalTet"),
362  GetOperatorFactory().RegisterCreatorFunction(
363  OperatorKey(ePyramid, ePhysDeriv, eIterPerExp,false),
364  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Pyr"),
365  GetOperatorFactory().RegisterCreatorFunction(
366  OperatorKey(ePrism, ePhysDeriv, eIterPerExp,false),
367  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Prism"),
368  GetOperatorFactory().RegisterCreatorFunction(
369  OperatorKey(ePrism, ePhysDeriv, eIterPerExp,true),
370  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_NodalPrism"),
371  GetOperatorFactory().RegisterCreatorFunction(
372  OperatorKey(eHexahedron, ePhysDeriv, eIterPerExp,false),
373  PhysDeriv_IterPerExp::create, "PhysDeriv_IterPerExp_Hex")
374 };
375 
376 
377 /**
378  * @brief Phys deriv operator using original LocalRegions implementation.
379  */
380 class PhysDeriv_NoCollection : public Operator
381 {
382  public:
383  OPERATOR_CREATE(PhysDeriv_NoCollection)
384 
385  virtual ~PhysDeriv_NoCollection()
386  {
387  }
388 
389  virtual void operator()(
390  const Array<OneD, const NekDouble> &input,
391  Array<OneD, NekDouble> &output0,
392  Array<OneD, NekDouble> &output1,
393  Array<OneD, NekDouble> &output2,
394  Array<OneD, NekDouble> &wsp)
395  {
396  boost::ignore_unused(wsp);
397 
398  const int nPhys = m_expList[0]->GetTotPoints();
399  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
400 
401  // calculate local derivatives
402  switch (m_expList[0]->GetShapeDimension())
403  {
404  case 1:
405  {
406  for (int i = 0; i < m_numElmt; ++i)
407  {
408  m_expList[i]->PhysDeriv(input + i*nPhys,
409  tmp0 = output0 + i*nPhys);
410  }
411  break;
412  }
413  case 2:
414  {
415  for (int i = 0; i < m_numElmt; ++i)
416  {
417  m_expList[i]->PhysDeriv(input + i*nPhys,
418  tmp0 = output0 + i*nPhys,
419  tmp1 = output1 + i*nPhys);
420  }
421  break;
422  }
423  case 3:
424  {
425  for (int i = 0; i < m_numElmt; ++i)
426  {
427  m_expList[i]->PhysDeriv(input + i*nPhys,
428  tmp0 = output0 + i*nPhys,
429  tmp1 = output1 + i*nPhys,
430  tmp2 = output2 + i*nPhys);
431  }
432  break;
433  }
434  default:
435  ASSERTL0(false, "Unknown dimension.");
436  }
437  }
438 
439  virtual void operator()(
440  int dir,
441  const Array<OneD, const NekDouble> &input,
442  Array<OneD, NekDouble> &output,
443  Array<OneD, NekDouble> &wsp)
444  {
445  boost::ignore_unused(wsp);
446 
447  const int nPhys = m_expList[0]->GetTotPoints();
448  Array<OneD, NekDouble> tmp;
449 
450  // calculate local derivatives
451  for (int i = 0; i < m_numElmt; ++i)
452  {
453  m_expList[i]->PhysDeriv(dir, input + i*nPhys,
454  tmp = output + i*nPhys);
455  }
456  }
457 
458  protected:
459  vector<StdRegions::StdExpansionSharedPtr> m_expList;
460 
461  private:
462  PhysDeriv_NoCollection(
463  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
464  CoalescedGeomDataSharedPtr pGeomData)
465  : Operator(pCollExp,pGeomData)
466  {
467  m_expList = pCollExp;
468  }
469 };
470 
471 /// Factory initialisation for the PhysDeriv_NoCollection operators
472 OperatorKey PhysDeriv_NoCollection::m_typeArr[] =
473 {
474  GetOperatorFactory().RegisterCreatorFunction(
475  OperatorKey(eSegment, ePhysDeriv, eNoCollection,false),
476  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Seg"),
477  GetOperatorFactory().RegisterCreatorFunction(
478  OperatorKey(eTriangle, ePhysDeriv, eNoCollection,false),
479  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Tri"),
480  GetOperatorFactory().RegisterCreatorFunction(
481  OperatorKey(eTriangle, ePhysDeriv, eNoCollection,true),
482  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalTri"),
483  GetOperatorFactory().RegisterCreatorFunction(
484  OperatorKey(eQuadrilateral, ePhysDeriv, eNoCollection,false),
485  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Quad"),
486  GetOperatorFactory().RegisterCreatorFunction(
487  OperatorKey(eTetrahedron, ePhysDeriv, eNoCollection,false),
488  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Tet"),
489  GetOperatorFactory().RegisterCreatorFunction(
490  OperatorKey(eTetrahedron, ePhysDeriv, eNoCollection,true),
491  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalTet"),
492  GetOperatorFactory().RegisterCreatorFunction(
493  OperatorKey(ePyramid, ePhysDeriv, eNoCollection,false),
494  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Pyr"),
495  GetOperatorFactory().RegisterCreatorFunction(
496  OperatorKey(ePrism, ePhysDeriv, eNoCollection,false),
497  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Prism"),
498  GetOperatorFactory().RegisterCreatorFunction(
499  OperatorKey(ePrism, ePhysDeriv, eNoCollection,true),
500  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalPrism"),
501  GetOperatorFactory().RegisterCreatorFunction(
502  OperatorKey(eHexahedron, ePhysDeriv, eNoCollection,false),
503  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Hex")
504 };
505 
506 
507 /**
508  * @brief Phys deriv operator using sum-factorisation (Segment)
509  */
510 class PhysDeriv_SumFac_Seg : public Operator
511 {
512  public:
513  OPERATOR_CREATE(PhysDeriv_SumFac_Seg)
514 
515  virtual ~PhysDeriv_SumFac_Seg()
516  {
517  }
518 
519  virtual void operator()(
520  const Array<OneD, const NekDouble> &input,
521  Array<OneD, NekDouble> &output0,
522  Array<OneD, NekDouble> &output1,
523  Array<OneD, NekDouble> &output2,
524  Array<OneD, NekDouble> &wsp)
525  {
526  const int nqcol = m_nquad0*m_numElmt;
527 
528  ASSERTL1(wsp.num_elements() == m_wspSize,
529  "Incorrect workspace size");
530  ASSERTL1(input.num_elements() >= nqcol,
531  "Incorrect input size");
532 
533  Array<OneD, NekDouble> diff0(nqcol, wsp);
534 
535  Blas::Dgemm('N', 'N', m_nquad0, m_numElmt,
536  m_nquad0, 1.0, m_Deriv0, m_nquad0,
537  input.get(), m_nquad0, 0.0,
538  diff0.get(), m_nquad0);
539 
540  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
541 
542  if (m_coordim == 2)
543  {
544  Vmath::Vmul (nqcol, m_derivFac[1], 1, diff0, 1, output1, 1);
545  }
546  else if (m_coordim == 3)
547  {
548  Vmath::Vmul (nqcol, m_derivFac[1], 1, diff0, 1, output1, 1);
549  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output2, 1);
550  }
551  }
552 
553  virtual void operator()(
554  int dir,
555  const Array<OneD, const NekDouble> &input,
556  Array<OneD, NekDouble> &output,
557  Array<OneD, NekDouble> &wsp)
558  {
559  const int nqcol = m_nquad0*m_numElmt;
560 
561  ASSERTL1(wsp.num_elements() == m_wspSize,
562  "Incorrect workspace size");
563  ASSERTL1(input.num_elements() >= nqcol,
564  "Incorrect input size");
565 
566  Array<OneD, NekDouble> diff0(nqcol, wsp);
567 
568  Blas::Dgemm('N', 'N', m_nquad0, m_numElmt,
569  m_nquad0, 1.0, m_Deriv0, m_nquad0,
570  input.get(), m_nquad0, 0.0,
571  diff0.get(), m_nquad0);
572 
573  Vmath::Vmul(nqcol, m_derivFac[dir], 1, diff0, 1, output, 1);
574  }
575 
576  protected:
577  int m_coordim;
578  const int m_nquad0;
579  Array<TwoD, const NekDouble> m_derivFac;
580  NekDouble *m_Deriv0;
581 
582  private:
583  PhysDeriv_SumFac_Seg(
584  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
585  CoalescedGeomDataSharedPtr pGeomData)
586  : Operator (pCollExp, pGeomData),
587  m_nquad0 (m_stdExp->GetNumPoints(0))
588  {
589  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
590  m_coordim = pCollExp[0]->GetCoordim();
591 
592  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
593 
594  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
595  m_wspSize = m_nquad0*m_numElmt;
596  }
597 
598 };
599 
600 /// Factory initialisation for the PhysDeriv_SumFac_Seg operators
601 OperatorKey PhysDeriv_SumFac_Seg::m_type = GetOperatorFactory().
602  RegisterCreatorFunction(
603  OperatorKey(eSegment, ePhysDeriv, eSumFac,false),
604  PhysDeriv_SumFac_Seg::create, "PhysDeriv_SumFac_Seg");
605 
606 
607 
608 /**
609  * @brief Phys deriv operator using sum-factorisation (Quad)
610  */
611 class PhysDeriv_SumFac_Quad : public Operator
612 {
613  public:
614  OPERATOR_CREATE(PhysDeriv_SumFac_Quad)
615 
616  virtual ~PhysDeriv_SumFac_Quad()
617  {
618  }
619 
620  virtual void operator()(
621  const Array<OneD, const NekDouble> &input,
622  Array<OneD, NekDouble> &output0,
623  Array<OneD, NekDouble> &output1,
624  Array<OneD, NekDouble> &output2,
625  Array<OneD, NekDouble> &wsp)
626  {
627  const int nqtot = m_nquad0 * m_nquad1;
628  const int nqcol = nqtot*m_numElmt;
629 
630  ASSERTL1(wsp.num_elements() == m_wspSize,
631  "Incorrect workspace size");
632  ASSERTL1(input.num_elements() >= nqcol,
633  "Incorrect input size");
634 
635  Array<OneD, NekDouble> diff0(nqcol, wsp );
636  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
637 
638  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
639  m_nquad0, 1.0, m_Deriv0, m_nquad0,
640  input.get(), m_nquad0, 0.0,
641  diff0.get(), m_nquad0);
642 
643  int cnt = 0;
644  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
645  {
646  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
647  input.get() + cnt, m_nquad0,
648  m_Deriv1, m_nquad1, 0.0,
649  diff1.get() + cnt, m_nquad0);
650  }
651 
652  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
653  Vmath::Vvtvp (nqcol, m_derivFac[1], 1, diff1, 1, output0, 1,
654  output0, 1);
655  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output1, 1);
656  Vmath::Vvtvp (nqcol, m_derivFac[3], 1, diff1, 1, output1, 1,
657  output1, 1);
658 
659  if (m_coordim == 3)
660  {
661  Vmath::Vmul (nqcol, m_derivFac[4], 1, diff0, 1, output2, 1);
662  Vmath::Vvtvp (nqcol, m_derivFac[5], 1, diff1, 1, output2, 1,
663  output2, 1);
664  }
665  }
666 
667  virtual void operator()(
668  int dir,
669  const Array<OneD, const NekDouble> &input,
670  Array<OneD, NekDouble> &output,
671  Array<OneD, NekDouble> &wsp)
672  {
673  const int nqtot = m_nquad0 * m_nquad1;
674  const int nqcol = nqtot*m_numElmt;
675 
676  ASSERTL1(wsp.num_elements() == m_wspSize,
677  "Incorrect workspace size");
678  ASSERTL1(input.num_elements() >= nqcol,
679  "Incorrect input size");
680 
681  Array<OneD, NekDouble> diff0(nqcol, wsp );
682  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
683 
684  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
685  m_nquad0, 1.0, m_Deriv0, m_nquad0,
686  input.get(), m_nquad0, 0.0,
687  diff0.get(), m_nquad0);
688 
689  int cnt = 0;
690  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
691  {
692  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
693  input.get() + cnt, m_nquad0,
694  m_Deriv1, m_nquad1, 0.0,
695  diff1.get() + cnt, m_nquad0);
696  }
697 
698  Vmath::Vmul (nqcol, m_derivFac[2*dir] , 1, diff0, 1, output, 1);
699  Vmath::Vvtvp (nqcol, m_derivFac[2*dir+1], 1, diff1, 1, output, 1,
700  output, 1);
701  }
702 
703  protected:
704  int m_coordim;
705  const int m_nquad0;
706  const int m_nquad1;
707  Array<TwoD, const NekDouble> m_derivFac;
708  NekDouble *m_Deriv0;
709  NekDouble *m_Deriv1;
710 
711  private:
712  PhysDeriv_SumFac_Quad(
713  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
714  CoalescedGeomDataSharedPtr pGeomData)
715  : Operator (pCollExp, pGeomData),
716  m_nquad0 (m_stdExp->GetNumPoints(0)),
717  m_nquad1 (m_stdExp->GetNumPoints(1))
718  {
719  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
720  m_coordim = pCollExp[0]->GetCoordim();
721 
722  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
723 
724  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
725  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
726  m_wspSize = 2 * m_nquad0*m_nquad1*m_numElmt;
727  }
728 };
729 
730 /// Factory initialisation for the PhysDeriv_SumFac_Quad operators
731 OperatorKey PhysDeriv_SumFac_Quad::m_type = GetOperatorFactory().
732  RegisterCreatorFunction(
733  OperatorKey(eQuadrilateral, ePhysDeriv, eSumFac, false),
734  PhysDeriv_SumFac_Quad::create, "PhysDeriv_SumFac_Quad");
735 
736 
737 /**
738  * @brief Phys deriv operator using sum-factorisation (Tri)
739  */
740 class PhysDeriv_SumFac_Tri : public Operator
741 {
742  public:
743  OPERATOR_CREATE(PhysDeriv_SumFac_Tri)
744 
745  virtual ~PhysDeriv_SumFac_Tri()
746  {
747  }
748 
749  virtual void operator()(
750  const Array<OneD, const NekDouble> &input,
751  Array<OneD, NekDouble> &output0,
752  Array<OneD, NekDouble> &output1,
753  Array<OneD, NekDouble> &output2,
754  Array<OneD, NekDouble> &wsp)
755  {
756  const int nqtot = m_nquad0 * m_nquad1;
757  const int nqcol = nqtot*m_numElmt;
758 
759  ASSERTL1(wsp.num_elements() == m_wspSize,
760  "Incorrect workspace size");
761  ASSERTL1(input.num_elements() >= nqcol,
762  "Incorrect input size");
763 
764  Array<OneD, NekDouble> diff0(nqcol, wsp );
765  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
766 
767  // Tensor Product Derivative
768  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
769  m_nquad0, 1.0, m_Deriv0, m_nquad0,
770  input.get(), m_nquad0, 0.0,
771  diff0.get(), m_nquad0);
772 
773  int cnt = 0;
774  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
775  {
776  // scale diff0 by geometric factor: 2/(1-z1)
777  Vmath::Vmul(nqtot,&m_fac1[0],1,diff0.get()+cnt,1,
778  diff0.get()+cnt,1);
779 
780  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
781  input.get() + cnt, m_nquad0,
782  m_Deriv1, m_nquad1, 0.0,
783  diff1.get() + cnt, m_nquad0);
784 
785  // add to diff1 by diff0 scaled by: (1_z0)/(1-z1)
786  Vmath::Vvtvp(nqtot,m_fac0.get(),1,diff0.get()+cnt,1,
787  diff1.get()+cnt,1,diff1.get()+cnt,1);
788  }
789 
790 
791  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
792  Vmath::Vvtvp (nqcol, m_derivFac[1], 1, diff1, 1,
793  output0, 1, output0, 1);
794  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output1, 1);
795  Vmath::Vvtvp (nqcol, m_derivFac[3], 1, diff1, 1,
796  output1, 1, output1, 1);
797 
798  if (m_coordim == 3)
799  {
800  Vmath::Vmul (nqcol, m_derivFac[4], 1, diff0, 1,
801  output2, 1);
802  Vmath::Vvtvp (nqcol, m_derivFac[5], 1, diff1, 1,
803  output2, 1, output2, 1);
804  }
805  }
806 
807  virtual void operator()(
808  int dir,
809  const Array<OneD, const NekDouble> &input,
810  Array<OneD, NekDouble> &output,
811  Array<OneD, NekDouble> &wsp)
812  {
813  const int nqtot = m_nquad0 * m_nquad1;
814  const int nqcol = nqtot*m_numElmt;
815 
816  ASSERTL1(wsp.num_elements() == m_wspSize,
817  "Incorrect workspace size");
818  ASSERTL1(input.num_elements() >= nqcol,
819  "Incorrect input size");
820 
821  Array<OneD, NekDouble> diff0(nqcol, wsp );
822  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
823 
824  // Tensor Product Derivative
825  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
826  m_nquad0, 1.0, m_Deriv0, m_nquad0,
827  input.get(), m_nquad0, 0.0,
828  diff0.get(), m_nquad0);
829 
830  int cnt = 0;
831  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
832  {
833  // scale diff0 by geometric factor: 2/(1-z1)
834  Vmath::Vmul(nqtot,&m_fac1[0],1,diff0.get()+cnt,1,
835  diff0.get()+cnt,1);
836 
837  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
838  input.get() + cnt, m_nquad0,
839  m_Deriv1, m_nquad1, 0.0,
840  diff1.get() + cnt, m_nquad0);
841 
842  // add to diff1 by diff0 scaled by: (1_z0)/(1-z1)
843  Vmath::Vvtvp(nqtot,m_fac0.get(),1,diff0.get()+cnt,1,
844  diff1.get()+cnt,1,diff1.get()+cnt,1);
845  }
846 
847 
848  Vmath::Vmul (nqcol, m_derivFac[2*dir] , 1, diff0, 1, output, 1);
849  Vmath::Vvtvp (nqcol, m_derivFac[2*dir+1], 1, diff1, 1, output, 1,
850  output, 1);
851  }
852 
853  protected:
854  int m_coordim;
855  const int m_nquad0;
856  const int m_nquad1;
857  Array<TwoD, const NekDouble> m_derivFac;
858  NekDouble *m_Deriv0;
859  NekDouble *m_Deriv1;
860  Array<OneD, NekDouble> m_fac0;
861  Array<OneD, NekDouble> m_fac1;
862 
863  private:
864  PhysDeriv_SumFac_Tri(
865  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
866  CoalescedGeomDataSharedPtr pGeomData)
867  : Operator (pCollExp, pGeomData),
868  m_nquad0 (m_stdExp->GetNumPoints(0)),
869  m_nquad1 (m_stdExp->GetNumPoints(1))
870  {
871  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
872  m_coordim = pCollExp[0]->GetCoordim();
873 
874  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
875 
876  const Array<OneD, const NekDouble>& z0
877  = m_stdExp->GetBasis(0)->GetZ();
878  const Array<OneD, const NekDouble>& z1
879  = m_stdExp->GetBasis(1)->GetZ();
880  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
881  // set up geometric factor: 0.5*(1+z0)
882  for (int i = 0; i < m_nquad0; ++i)
883  {
884  for(int j = 0; j < m_nquad1; ++j)
885  {
886  m_fac0[i+j*m_nquad0] = 0.5*(1+z0[i]);
887  }
888  }
889 
890  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
891  // set up geometric factor: 2/(1-z1)
892  for (int i = 0; i < m_nquad0; ++i)
893  {
894  for(int j = 0; j < m_nquad1; ++j)
895  {
896  m_fac1[i+j*m_nquad0] = 2.0/(1-z1[j]);
897  }
898  }
899 
900 
901  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
902  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
903  m_wspSize = 2 * m_nquad0*m_nquad1*m_numElmt;
904  }
905 };
906 
907 /// Factory initialisation for the PhysDeriv_SumFac_Tri operators
908 OperatorKey PhysDeriv_SumFac_Tri::m_typeArr[] =
909 {
910  GetOperatorFactory().RegisterCreatorFunction(
911  OperatorKey(eTriangle, ePhysDeriv, eSumFac,false),
912  PhysDeriv_SumFac_Tri::create, "PhysDeriv_SumFac_Tri"),
913  GetOperatorFactory().RegisterCreatorFunction(
914  OperatorKey(eTriangle, ePhysDeriv, eSumFac,true),
915  PhysDeriv_SumFac_Tri::create, "PhysDeriv_SumFac_NodalTri")
916 };
917 
918 
919 /**
920  * @brief Phys deriv operator using sum-factorisation (Hex)
921  */
922 class PhysDeriv_SumFac_Hex : public Operator
923 {
924  public:
925  OPERATOR_CREATE(PhysDeriv_SumFac_Hex)
926 
927  virtual ~PhysDeriv_SumFac_Hex()
928  {
929  }
930 
931  virtual void operator()(
932  const Array<OneD, const NekDouble> &input,
933  Array<OneD, NekDouble> &output0,
934  Array<OneD, NekDouble> &output1,
935  Array<OneD, NekDouble> &output2,
936  Array<OneD, NekDouble> &wsp)
937  {
938  int nPhys = m_stdExp->GetTotPoints();
939  int ntot = m_numElmt*nPhys;
940  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
941  Array<OneD, Array<OneD, NekDouble> > Diff(3);
942  Array<OneD, Array<OneD, NekDouble> > out(3);
943  out[0] = output0; out[1] = output1; out[2] = output2;
944 
945  for(int i = 0; i < 3; ++i)
946  {
947  Diff[i] = wsp + i*ntot;
948  }
949 
950  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
951  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
952  m_nquad0,0.0,&Diff[0][0],m_nquad0);
953 
954  for(int i = 0; i < m_numElmt; ++i)
955  {
956  for (int j = 0; j < m_nquad2; ++j)
957  {
958  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
959  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
960  m_nquad0, m_Deriv1, m_nquad1, 0.0,
961  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
962  m_nquad0);
963  }
964 
965  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
966  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
967  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
968  m_nquad0*m_nquad1);
969  }
970 
971  // calculate full derivative
972  for(int i = 0; i < m_coordim; ++i)
973  {
974  Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
975  for(int j = 1; j < 3; ++j)
976  {
977  Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
978  Diff[j], 1,
979  out[i], 1,
980  out[i], 1);
981  }
982  }
983  }
984 
985  virtual void operator()(
986  int dir,
987  const Array<OneD, const NekDouble> &input,
988  Array<OneD, NekDouble> &output,
989  Array<OneD, NekDouble> &wsp)
990  {
991  int nPhys = m_stdExp->GetTotPoints();
992  int ntot = m_numElmt*nPhys;
993  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
994  Array<OneD, Array<OneD, NekDouble> > Diff(3);
995 
996  for(int i = 0; i < 3; ++i)
997  {
998  Diff[i] = wsp + i*ntot;
999  }
1000 
1001  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1002  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1003  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1004 
1005  for(int i = 0; i < m_numElmt; ++i)
1006  {
1007  for (int j = 0; j < m_nquad2; ++j)
1008  {
1009  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1010  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1011  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1012  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1013  m_nquad0);
1014  }
1015 
1016  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1017  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1018  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1019  m_nquad0*m_nquad1);
1020  }
1021 
1022  // calculate full derivative
1023  Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
1024  for(int j = 1; j < 3; ++j)
1025  {
1026  Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
1027  Diff[j], 1,
1028  output, 1,
1029  output, 1);
1030  }
1031  }
1032 
1033  protected:
1034  Array<TwoD, const NekDouble> m_derivFac;
1035  int m_coordim;
1036  const int m_nquad0;
1037  const int m_nquad1;
1038  const int m_nquad2;
1039  NekDouble *m_Deriv0;
1040  NekDouble *m_Deriv1;
1041  NekDouble *m_Deriv2;
1042 
1043  private:
1044  PhysDeriv_SumFac_Hex(
1045  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1046  CoalescedGeomDataSharedPtr pGeomData)
1047  : Operator(pCollExp, pGeomData),
1048  m_nquad0 (m_stdExp->GetNumPoints(0)),
1049  m_nquad1 (m_stdExp->GetNumPoints(1)),
1050  m_nquad2 (m_stdExp->GetNumPoints(2))
1051  {
1052  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1053 
1054  m_coordim = pCollExp[0]->GetCoordim();
1055 
1056  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1057 
1058  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1059  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1060  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1061 
1062  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1063  }
1064 };
1065 
1066 /// Factory initialisation for the PhysDeriv_SumFac_Hex operators
1067 OperatorKey PhysDeriv_SumFac_Hex::m_typeArr[] =
1068 {
1069  GetOperatorFactory().RegisterCreatorFunction(
1070  OperatorKey(eHexahedron, ePhysDeriv, eSumFac, false),
1071  PhysDeriv_SumFac_Hex::create, "PhysDeriv_SumFac_Hex")
1072 };
1073 
1074 
1075 /**
1076  * @brief Phys deriv operator using sum-factorisation (Tet)
1077  */
1078 class PhysDeriv_SumFac_Tet : public Operator
1079 {
1080  public:
1081  OPERATOR_CREATE(PhysDeriv_SumFac_Tet)
1082 
1083  virtual ~PhysDeriv_SumFac_Tet()
1084  {
1085  }
1086 
1087  virtual void operator()(
1088  const Array<OneD, const NekDouble> &input,
1089  Array<OneD, NekDouble> &output0,
1090  Array<OneD, NekDouble> &output1,
1091  Array<OneD, NekDouble> &output2,
1092  Array<OneD, NekDouble> &wsp)
1093  {
1094  int nPhys = m_stdExp->GetTotPoints();
1095  int ntot = m_numElmt*nPhys;
1096  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1097  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1098  Array<OneD, Array<OneD, NekDouble> > out(3);
1099  out[0] = output0; out[1] = output1; out[2] = output2;
1100 
1101  for(int i = 0; i < 3; ++i)
1102  {
1103  Diff[i] = wsp + i*ntot;
1104  }
1105 
1106  // dEta0
1107  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1108  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1109  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1110 
1111  // dEta2
1112  for(int i = 0; i < m_numElmt; ++i)
1113  {
1114  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1115  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1116  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1117  m_nquad0*m_nquad1);
1118  }
1119 
1120  for(int i = 0; i < m_numElmt; ++i)
1121  {
1122 
1123  // dEta1
1124  for (int j = 0; j < m_nquad2; ++j)
1125  {
1126  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1127  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1128  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1129  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1130  m_nquad0);
1131  }
1132 
1133  // dxi2 = (1 + eta_1)/(1 -eta_2)*dEta1 + dEta2
1134  Vmath::Vvtvp(nPhys, m_fac3.get(), 1,
1135  Diff[1].get() + i*nPhys, 1,
1136  Diff[2].get() + i*nPhys, 1,
1137  Diff[2].get() + i*nPhys, 1);
1138 
1139  // dxi1 = 2/(1 - eta_2) dEta1
1140  Vmath::Vmul(nPhys, m_fac2.get(), 1,
1141  Diff[1].get() + i*nPhys, 1,
1142  Diff[1].get() + i*nPhys, 1);
1143 
1144  // dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi1
1145  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1146  Diff[0].get() + i*nPhys, 1,
1147  Diff[1].get() + i*nPhys, 1,
1148  Diff[1].get() + i*nPhys, 1);
1149 
1150  // dxi2 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi2
1151  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1152  Diff[0].get() + i*nPhys, 1,
1153  Diff[2].get() + i*nPhys, 1,
1154  Diff[2].get() + i*nPhys, 1);
1155 
1156  // dxi0 = 4.0/((1-eta_1)(1-eta_2)) dEta0
1157  Vmath::Vmul(nPhys, m_fac0.get(), 1,
1158  Diff[0].get() + i*nPhys, 1,
1159  Diff[0].get() + i*nPhys, 1);
1160 
1161  }
1162 
1163  // calculate full derivative
1164  for(int i = 0; i < m_coordim; ++i)
1165  {
1166  Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
1167  for(int j = 1; j < 3; ++j)
1168  {
1169  Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
1170  Diff[j], 1, out[i], 1, out[i], 1);
1171  }
1172  }
1173  }
1174 
1175  virtual void operator()(
1176  int dir,
1177  const Array<OneD, const NekDouble> &input,
1178  Array<OneD, NekDouble> &output,
1179  Array<OneD, NekDouble> &wsp)
1180  {
1181  int nPhys = m_stdExp->GetTotPoints();
1182  int ntot = m_numElmt*nPhys;
1183  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1184  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1185 
1186  for(int i = 0; i < 3; ++i)
1187  {
1188  Diff[i] = wsp + i*ntot;
1189  }
1190 
1191  // dEta0
1192  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1193  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1194  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1195 
1196  // dEta2
1197  for(int i = 0; i < m_numElmt; ++i)
1198  {
1199  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1200  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1201  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1202  m_nquad0*m_nquad1);
1203  }
1204 
1205  for(int i = 0; i < m_numElmt; ++i)
1206  {
1207 
1208  // dEta1
1209  for (int j = 0; j < m_nquad2; ++j)
1210  {
1211  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1212  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1213  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1214  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1215  m_nquad0);
1216  }
1217 
1218  // dxi2 = (1 + eta_1)/(1 -eta_2)*dEta1 + dEta2
1219  Vmath::Vvtvp(nPhys, m_fac3.get(), 1,
1220  Diff[1].get() + i*nPhys, 1,
1221  Diff[2].get() + i*nPhys, 1,
1222  Diff[2].get() + i*nPhys, 1);
1223 
1224  // dxi1 = 2/(1 - eta_2) dEta1
1225  Vmath::Vmul(nPhys, m_fac2.get(), 1,
1226  Diff[1].get() + i*nPhys, 1,
1227  Diff[1].get() + i*nPhys, 1);
1228 
1229  // dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi1
1230  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1231  Diff[0].get() + i*nPhys, 1,
1232  Diff[1].get() + i*nPhys, 1,
1233  Diff[1].get() + i*nPhys, 1);
1234 
1235  // dxi2 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi2
1236  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1237  Diff[0].get() + i*nPhys, 1,
1238  Diff[2].get() + i*nPhys, 1,
1239  Diff[2].get() + i*nPhys, 1);
1240 
1241  // dxi0 = 4.0/((1-eta_1)(1-eta_2)) dEta0
1242  Vmath::Vmul(nPhys, m_fac0.get(), 1,
1243  Diff[0].get() + i*nPhys, 1,
1244  Diff[0].get() + i*nPhys, 1);
1245 
1246  }
1247 
1248  // calculate full derivative
1249  Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
1250  for(int j = 1; j < 3; ++j)
1251  {
1252  Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
1253  Diff[j], 1, output, 1, output, 1);
1254  }
1255  }
1256 
1257  protected:
1258  Array<TwoD, const NekDouble> m_derivFac;
1259  int m_coordim;
1260  const int m_nquad0;
1261  const int m_nquad1;
1262  const int m_nquad2;
1263  NekDouble *m_Deriv0;
1264  NekDouble *m_Deriv1;
1265  NekDouble *m_Deriv2;
1266  Array<OneD, NekDouble> m_fac0;
1267  Array<OneD, NekDouble> m_fac1;
1268  Array<OneD, NekDouble> m_fac2;
1269  Array<OneD, NekDouble> m_fac3;
1270 
1271  private:
1272  PhysDeriv_SumFac_Tet(
1273  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1274  CoalescedGeomDataSharedPtr pGeomData)
1275  : Operator(pCollExp, pGeomData),
1276  m_nquad0 (m_stdExp->GetNumPoints(0)),
1277  m_nquad1 (m_stdExp->GetNumPoints(1)),
1278  m_nquad2 (m_stdExp->GetNumPoints(2))
1279  {
1280  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1281 
1282  m_coordim = pCollExp[0]->GetCoordim();
1283 
1284  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1285 
1286  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1287  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1288  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1289 
1290  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1291 
1292  const Array<OneD, const NekDouble>& z0
1293  = m_stdExp->GetBasis(0)->GetZ();
1294  const Array<OneD, const NekDouble>& z1
1295  = m_stdExp->GetBasis(1)->GetZ();
1296  const Array<OneD, const NekDouble>& z2
1297  = m_stdExp->GetBasis(2)->GetZ();
1298 
1299  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1300  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1301  m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1302  m_fac3 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1303  // calculate 2.0/((1-eta_1)(1-eta_2))
1304  for (int i = 0; i < m_nquad0; ++i)
1305  {
1306  for(int j = 0; j < m_nquad1; ++j)
1307  {
1308  for(int k = 0; k < m_nquad2; ++k)
1309  {
1310 
1311  m_fac0[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1312  = 4.0/((1-z1[j])*(1-z2[k]));
1313  m_fac1[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1314  = 2.0*(1+z0[i])/((1-z1[j])*(1-z2[k]));
1315  m_fac2[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1316  = 2.0/(1-z2[k]);
1317  m_fac3[i + j*m_nquad0 + k*m_nquad0*m_nquad1]
1318  = (1+z1[j])/(1-z2[k]);
1319  }
1320  }
1321  }
1322 
1323  }
1324 };
1325 
1326 /// Factory initialisation for the PhysDeriv_SumFac_Tet operators
1327 OperatorKey PhysDeriv_SumFac_Tet::m_typeArr[] =
1328 {
1329  GetOperatorFactory().RegisterCreatorFunction(
1330  OperatorKey(eTetrahedron, ePhysDeriv, eSumFac, false),
1331  PhysDeriv_SumFac_Tet::create, "PhysDeriv_SumFac_Tet")
1332 };
1333 
1334 
1335 /**
1336  * @brief Phys deriv operator using sum-factorisation (Prism)
1337  */
1338 class PhysDeriv_SumFac_Prism : public Operator
1339 {
1340  public:
1341  OPERATOR_CREATE(PhysDeriv_SumFac_Prism)
1342 
1343  virtual ~PhysDeriv_SumFac_Prism()
1344  {
1345  }
1346 
1347  virtual void operator()(
1348  const Array<OneD, const NekDouble> &input,
1349  Array<OneD, NekDouble> &output0,
1350  Array<OneD, NekDouble> &output1,
1351  Array<OneD, NekDouble> &output2,
1352  Array<OneD, NekDouble> &wsp)
1353  {
1354  int nPhys = m_stdExp->GetTotPoints();
1355  int ntot = m_numElmt*nPhys;
1356  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1357  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1358  Array<OneD, Array<OneD, NekDouble> > out(3);
1359  out[0] = output0; out[1] = output1; out[2] = output2;
1360 
1361  for(int i = 0; i < 3; ++i)
1362  {
1363  Diff[i] = wsp + i*ntot;
1364  }
1365 
1366  // dEta0
1367  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1368  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1369  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1370 
1371  int cnt = 0;
1372  for(int i = 0; i < m_numElmt; ++i)
1373  {
1374 
1375  // dEta 1
1376  for (int j = 0; j < m_nquad2; ++j)
1377  {
1378  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1379  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1380  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1381  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1382  m_nquad0);
1383  }
1384 
1385  // dEta 2
1386  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1387  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1388  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1389  m_nquad0*m_nquad1);
1390 
1391  // dxi0 = 2/(1-eta_2) d Eta_0
1392  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
1393  Diff[0].get()+cnt,1);
1394 
1395  // dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
1396  Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
1397  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1398  cnt += nPhys;
1399  }
1400 
1401  // calculate full derivative
1402  for(int i = 0; i < m_coordim; ++i)
1403  {
1404  Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
1405  for(int j = 1; j < 3; ++j)
1406  {
1407  Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
1408  Diff[j], 1, out[i], 1, out[i], 1);
1409  }
1410  }
1411  }
1412 
1413  virtual void operator()(
1414  int dir,
1415  const Array<OneD, const NekDouble> &input,
1416  Array<OneD, NekDouble> &output,
1417  Array<OneD, NekDouble> &wsp)
1418  {
1419  int nPhys = m_stdExp->GetTotPoints();
1420  int ntot = m_numElmt*nPhys;
1421  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1422  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1423 
1424  for(int i = 0; i < 3; ++i)
1425  {
1426  Diff[i] = wsp + i*ntot;
1427  }
1428 
1429  // dEta0
1430  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1431  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1432  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1433 
1434  int cnt = 0;
1435  for(int i = 0; i < m_numElmt; ++i)
1436  {
1437 
1438  // dEta 1
1439  for (int j = 0; j < m_nquad2; ++j)
1440  {
1441  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1442  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1443  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1444  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1445  m_nquad0);
1446  }
1447 
1448  // dEta 2
1449  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1450  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1451  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1452  m_nquad0*m_nquad1);
1453 
1454  // dxi0 = 2/(1-eta_2) d Eta_0
1455  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
1456  Diff[0].get()+cnt,1);
1457 
1458  // dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
1459  Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
1460  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1461  cnt += nPhys;
1462  }
1463 
1464  // calculate full derivative
1465  Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
1466  for(int j = 1; j < 3; ++j)
1467  {
1468  Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
1469  Diff[j], 1, output, 1, output, 1);
1470  }
1471  }
1472 
1473  protected:
1474  Array<TwoD, const NekDouble> m_derivFac;
1475  int m_coordim;
1476  const int m_nquad0;
1477  const int m_nquad1;
1478  const int m_nquad2;
1479  NekDouble *m_Deriv0;
1480  NekDouble *m_Deriv1;
1481  NekDouble *m_Deriv2;
1482  Array<OneD, NekDouble> m_fac0;
1483  Array<OneD, NekDouble> m_fac1;
1484 
1485  private:
1486  PhysDeriv_SumFac_Prism(
1487  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1488  CoalescedGeomDataSharedPtr pGeomData)
1489  : Operator(pCollExp, pGeomData),
1490  m_nquad0 (m_stdExp->GetNumPoints(0)),
1491  m_nquad1 (m_stdExp->GetNumPoints(1)),
1492  m_nquad2 (m_stdExp->GetNumPoints(2))
1493  {
1494  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1495 
1496  m_coordim = pCollExp[0]->GetCoordim();
1497 
1498  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1499 
1500  const Array<OneD, const NekDouble>& z0
1501  = m_stdExp->GetBasis(0)->GetZ();
1502  const Array<OneD, const NekDouble>& z2
1503  = m_stdExp->GetBasis(2)->GetZ();
1504  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1505  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1506  for (int i = 0; i < m_nquad0; ++i)
1507  {
1508  for(int j = 0; j < m_nquad1; ++j)
1509  {
1510  for(int k = 0; k < m_nquad2; ++k)
1511  {
1512  m_fac0[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
1513  2.0/(1-z2[k]);
1514  m_fac1[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
1515  0.5*(1+z0[i]);
1516  }
1517  }
1518  }
1519 
1520 
1521 
1522  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1523  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1524  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1525 
1526  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1527  }
1528 };
1529 
1530 /// Factory initialisation for the PhysDeriv_SumFac_Prism operators
1531 OperatorKey PhysDeriv_SumFac_Prism::m_typeArr[] = {
1532  GetOperatorFactory().RegisterCreatorFunction(
1533  OperatorKey(ePrism, ePhysDeriv, eSumFac, false),
1534  PhysDeriv_SumFac_Prism::create, "PhysDeriv_SumFac_Prism")
1535 };
1536 
1537 
1538 /**
1539  * @brief Phys deriv operator using sum-factorisation (Pyramid)
1540  */
1541 class PhysDeriv_SumFac_Pyr : public Operator
1542 {
1543  public:
1544  OPERATOR_CREATE(PhysDeriv_SumFac_Pyr)
1545 
1546  virtual ~PhysDeriv_SumFac_Pyr()
1547  {
1548  }
1549 
1550  virtual void operator()(
1551  const Array<OneD, const NekDouble> &input,
1552  Array<OneD, NekDouble> &output0,
1553  Array<OneD, NekDouble> &output1,
1554  Array<OneD, NekDouble> &output2,
1555  Array<OneD, NekDouble> &wsp)
1556  {
1557  int nPhys = m_stdExp->GetTotPoints();
1558  int ntot = m_numElmt*nPhys;
1559  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1560  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1561  Array<OneD, Array<OneD, NekDouble> > out(3);
1562  out[0] = output0; out[1] = output1; out[2] = output2;
1563 
1564  for(int i = 0; i < 3; ++i)
1565  {
1566  Diff[i] = wsp + i*ntot;
1567  }
1568 
1569  // dEta0
1570  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1571  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1572  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1573 
1574  int cnt = 0;
1575  for(int i = 0; i < m_numElmt; ++i)
1576  {
1577 
1578  // dEta 1
1579  for (int j = 0; j < m_nquad2; ++j)
1580  {
1581  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1582  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1583  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1584  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1585  m_nquad0);
1586  }
1587 
1588  // dEta 2
1589  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1590  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1591  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1592  m_nquad0*m_nquad1);
1593 
1594  // dxi0 = 2/(1-eta_2) d Eta_0
1595  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
1596  Diff[0].get()+cnt,1);
1597 
1598  // dxi1 = 2/(1-eta_2) d Eta_1
1599  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
1600  Diff[1].get()+cnt,1);
1601 
1602  // dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
1603  Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
1604  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1605  // dxi2 += (1+eta1)/(1-eta_2) d Eta_1
1606  Vmath::Vvtvp(nPhys,&m_fac2[0],1,Diff[1].get()+cnt,1,
1607  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1608  cnt += nPhys;
1609  }
1610 
1611  // calculate full derivative
1612  for(int i = 0; i < m_coordim; ++i)
1613  {
1614  Vmath::Vmul(ntot,m_derivFac[i*3],1,Diff[0],1,out[i],1);
1615  for(int j = 1; j < 3; ++j)
1616  {
1617  Vmath::Vvtvp (ntot, m_derivFac[i*3+j], 1,
1618  Diff[j], 1, out[i], 1, out[i], 1);
1619  }
1620  }
1621  }
1622 
1623  virtual void operator()(
1624  int dir,
1625  const Array<OneD, const NekDouble> &input,
1626  Array<OneD, NekDouble> &output,
1627  Array<OneD, NekDouble> &wsp)
1628  {
1629  int nPhys = m_stdExp->GetTotPoints();
1630  int ntot = m_numElmt*nPhys;
1631  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1632  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1633 
1634  for(int i = 0; i < 3; ++i)
1635  {
1636  Diff[i] = wsp + i*ntot;
1637  }
1638 
1639  // dEta0
1640  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1641  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1642  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1643 
1644  int cnt = 0;
1645  for(int i = 0; i < m_numElmt; ++i)
1646  {
1647  // dEta 1
1648  for (int j = 0; j < m_nquad2; ++j)
1649  {
1650  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1651  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1652  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1653  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1654  m_nquad0);
1655  }
1656 
1657  // dEta 2
1658  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1659  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1660  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1661  m_nquad0*m_nquad1);
1662 
1663  // dxi0 = 2/(1-eta_2) d Eta_0
1664  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[0].get()+cnt,1,
1665  Diff[0].get()+cnt,1);
1666 
1667  // dxi1 = 2/(1-eta_2) d Eta_1
1668  Vmath::Vmul(nPhys,&m_fac0[0],1,Diff[1].get()+cnt,1,
1669  Diff[1].get()+cnt,1);
1670 
1671  // dxi2 = (1+eta0)/(1-eta_2) d Eta_0 + d/dEta2;
1672  Vmath::Vvtvp(nPhys,&m_fac1[0],1,Diff[0].get()+cnt,1,
1673  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1674  // dxi2 = (1+eta1)/(1-eta_2) d Eta_1 + d/dEta2;
1675  Vmath::Vvtvp(nPhys,&m_fac2[0],1,Diff[1].get()+cnt,1,
1676  Diff[2].get()+cnt,1,Diff[2].get()+cnt,1);
1677  cnt += nPhys;
1678  }
1679 
1680  // calculate full derivative
1681  Vmath::Vmul(ntot,m_derivFac[dir*3],1,Diff[0],1,output,1);
1682  for(int j = 1; j < 3; ++j)
1683  {
1684  Vmath::Vvtvp (ntot, m_derivFac[dir*3+j], 1,
1685  Diff[j], 1, output, 1, output, 1);
1686  }
1687  }
1688 
1689  protected:
1690  Array<TwoD, const NekDouble> m_derivFac;
1691  int m_coordim;
1692  const int m_nquad0;
1693  const int m_nquad1;
1694  const int m_nquad2;
1695  NekDouble *m_Deriv0;
1696  NekDouble *m_Deriv1;
1697  NekDouble *m_Deriv2;
1698  Array<OneD, NekDouble> m_fac0;
1699  Array<OneD, NekDouble> m_fac1;
1700  Array<OneD, NekDouble> m_fac2;
1701 
1702  private:
1703  PhysDeriv_SumFac_Pyr(
1704  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1705  CoalescedGeomDataSharedPtr pGeomData)
1706  : Operator(pCollExp, pGeomData),
1707  m_nquad0 (m_stdExp->GetNumPoints(0)),
1708  m_nquad1 (m_stdExp->GetNumPoints(1)),
1709  m_nquad2 (m_stdExp->GetNumPoints(2))
1710  {
1711  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1712 
1713  m_coordim = pCollExp[0]->GetCoordim();
1714 
1715  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1716 
1717  const Array<OneD, const NekDouble>& z0
1718  = m_stdExp->GetBasis(0)->GetZ();
1719  const Array<OneD, const NekDouble>& z1
1720  = m_stdExp->GetBasis(1)->GetZ();
1721  const Array<OneD, const NekDouble>& z2
1722  = m_stdExp->GetBasis(2)->GetZ();
1723  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1724  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1725  m_fac2 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1726 
1727  int nq0_nq1 = m_nquad0*m_nquad1;
1728  for (int i = 0; i < m_nquad0; ++i)
1729  {
1730  for(int j = 0; j < m_nquad1; ++j)
1731  {
1732  int ifac = i+j*m_nquad0;
1733  for(int k = 0; k < m_nquad2; ++k)
1734  {
1735  m_fac0[ifac + k*nq0_nq1] =
1736  2.0/(1-z2[k]);
1737  m_fac1[ifac + k*nq0_nq1] =
1738  0.5*(1+z0[i]);
1739  m_fac2[ifac + k*nq0_nq1] =
1740  0.5*(1+z1[j]);
1741  }
1742  }
1743  }
1744 
1745  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1746  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1747  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1748 
1749  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1750  }
1751 };
1752 
1753 /// Factory initialisation for the PhysDeriv_SumFac_Pyr operators
1754 OperatorKey PhysDeriv_SumFac_Pyr::m_typeArr[] = {
1755  GetOperatorFactory().RegisterCreatorFunction(
1756  OperatorKey(ePyramid, ePhysDeriv, eSumFac, false),
1757  PhysDeriv_SumFac_Pyr::create, "PhysDeriv_SumFac_Pyr")
1758 };
1759 
1760 }
1761 }
Nektar::Collections::PhysDeriv_SumFac_Prism
Phys deriv operator using sum-factorisation (Prism)
Definition: PhysDeriv.cpp:1338
Nektar::Collections::CoalescedGeomDataSharedPtr
std::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
Definition: CoalescedGeomData.h:74
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::Collections::PhysDeriv_StdMat::m_dim
int m_dim
Definition: PhysDeriv.cpp:148
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:246
Nektar::Collections::PhysDeriv_SumFac_Quad::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:706
Nektar::Collections::PhysDeriv_SumFac_Quad::PhysDeriv_SumFac_Quad
PhysDeriv_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:712
Nektar::Collections::PhysDeriv_SumFac_Quad::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:704
Nektar::Collections::PhysDeriv_SumFac_Hex::PhysDeriv_SumFac_Hex
PhysDeriv_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1044
Nektar::Collections::PhysDeriv_SumFac_Tet::PhysDeriv_SumFac_Tet
PhysDeriv_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1272
Nektar::Collections::PhysDeriv_SumFac_Prism::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:1347
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1696
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1690
Nektar::Collections::PhysDeriv_IterPerExp::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:317
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: PhysDeriv.cpp:1699
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1693
Nektar::Collections::PhysDeriv_SumFac_Seg::PhysDeriv_SumFac_Seg
PhysDeriv_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:583
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1691
Nektar::Collections::PhysDeriv_SumFac_Seg::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:553
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:83
Nektar::Collections::PhysDeriv_SumFac_Tri::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:749
Nektar::Collections::OperatorKey
std::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:161
Nektar::Collections::PhysDeriv_IterPerExp::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:280
Nektar::Collections::PhysDeriv_SumFac_Seg
Phys deriv operator using sum-factorisation (Segment)
Definition: PhysDeriv.cpp:510
Nektar::Collections::PhysDeriv_NoCollection::m_expList
vector< StdRegions::StdExpansionSharedPtr > m_expList
Definition: PhysDeriv.cpp:459
Operator.h
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: PhysDeriv.cpp:1268
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:445
Nektar::Collections::PhysDeriv_SumFac_Hex::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1034
Nektar::Collections::PhysDeriv_SumFac_Prism::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1481
Nektar::Collections::PhysDeriv_SumFac_Quad::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:707
Nektar::Collections::PhysDeriv_StdMat::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:149
std
STL namespace.
Nektar::Collections::PhysDeriv_SumFac_Hex::m_nquad2
const int m_nquad2
Definition: PhysDeriv.cpp:1038
Nektar::Collections::PhysDeriv_IterPerExp
Phys deriv operator using element-wise operation.
Definition: PhysDeriv.cpp:229
Nektar::Collections::PhysDeriv_SumFac_Tri::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:807
Collection.h
Nektar::Collections::PhysDeriv_SumFac_Tet::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1261
Nektar::Collections::PhysDeriv_SumFac_Prism::PhysDeriv_SumFac_Prism
PhysDeriv_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1486
Nektar::Collections::PhysDeriv_SumFac_Tri::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:860
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 A[m x n], B[n x k], C[m x k].
Definition: Blas.hpp:213
Nektar::Collections::PhysDeriv_SumFac_Tri
Phys deriv operator using sum-factorisation (Tri)
Definition: PhysDeriv.cpp:740
Nektar::Collections::Operator
Base class for operators on a collection of elements.
Definition: Operator.h:109
Nektar::Collections::PhysDeriv_SumFac_Hex::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1035
Nektar::Collections::PhysDeriv_SumFac_Prism::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1477
Nektar::Collections::PhysDeriv_SumFac_Seg::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:578
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1692
Nektar::Collections::PhysDeriv_SumFac_Quad::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:620
Nektar::Collections::PhysDeriv_SumFac_Prism::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:1482
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1041
Nektar::Collections::PhysDeriv_SumFac_Tri::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:857
Nektar::Collections::PhysDeriv_StdMat
Phys deriv operator using standard matrix approach.
Definition: PhysDeriv.cpp:56
Nektar::Collections::PhysDeriv_SumFac_Prism::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1479
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1697
Nektar::Collections::PhysDeriv_SumFac_Seg::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:580
Nektar::Collections::PhysDeriv_SumFac_Quad::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:705
Nektar::Collections::PhysDeriv_SumFac_Hex::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:931
Nektar::Collections::PhysDeriv_SumFac_Tet::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:1087
Nektar::Collections::PhysDeriv_StdMat::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:108
Nektar::Collections::PhysDeriv_IterPerExp::m_dim
int m_dim
Definition: PhysDeriv.cpp:318
Nektar::Collections::PhysDeriv_SumFac_Tet::m_nquad2
const int m_nquad2
Definition: PhysDeriv.cpp:1262
Nektar::Collections::PhysDeriv_SumFac_Tet::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1260
Nektar::Collections::PhysDeriv_SumFac_Seg::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:579
Nektar::Collections::PhysDeriv_SumFac_Hex::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:985
Nektar::Collections::PhysDeriv_SumFac_Prism::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:1413
Nektar::Collections::PhysDeriv_SumFac_Prism::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: PhysDeriv.cpp:1483
Nektar::Collections::PhysDeriv_SumFac_Tet::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:1175
Nektar::Collections::PhysDeriv_SumFac_Prism::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1474
Nektar::Collections::PhysDeriv_SumFac_Prism::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1475
Nektar::Collections::PhysDeriv_SumFac_Quad::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:667
Nektar::Collections::PhysDeriv_SumFac_Tri::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:858
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::Collections::GetOperatorFactory
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:108
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1039
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1040
Nektar::Collections::PhysDeriv_SumFac_Pyr
Phys deriv operator using sum-factorisation (Pyramid)
Definition: PhysDeriv.cpp:1541
Nektar::Collections::PhysDeriv_SumFac_Tet
Phys deriv operator using sum-factorisation (Tet)
Definition: PhysDeriv.cpp:1078
Nektar::Collections::PhysDeriv_SumFac_Hex::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1037
Nektar::Collections::PhysDeriv_SumFac_Prism::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1480
Nektar::Collections::PhysDeriv_SumFac_Seg::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:519
Nektar::Collections::PhysDeriv_SumFac_Hex::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1036
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: PhysDeriv.cpp:1267
Nektar::Collections::PhysDeriv_SumFac_Quad::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:708
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:1266
Nektar::Collections::PhysDeriv_SumFac_Tet::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1259
Nektar::Collections::PhysDeriv_NoCollection::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:439
Nektar::Collections::PhysDeriv_SumFac_Tri::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:854
OPERATOR_CREATE
#define OPERATOR_CREATE(cname)
Definition: Operator.h:45
Nektar::Collections::PhysDeriv_SumFac_Tet::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1258
Nektar::Collections::PhysDeriv_SumFac_Pyr::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:1550
Nektar::Collections::PhysDeriv_IterPerExp::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:319
Nektar::Collections::PhysDeriv_SumFac_Tri::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:859
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac3
Array< OneD, NekDouble > m_fac3
Definition: PhysDeriv.cpp:1269
Nektar::Collections::PhysDeriv_StdMat::m_derivMat
Array< OneD, DNekMatSharedPtr > m_derivMat
Definition: PhysDeriv.cpp:146
Nektar::Collections::PhysDeriv_NoCollection
Phys deriv operator using original LocalRegions implementation.
Definition: PhysDeriv.cpp:380
Nektar::Collections::PhysDeriv_SumFac_Tri::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: PhysDeriv.cpp:861
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: PhysDeriv.cpp:1700
Nektar::Collections::PhysDeriv_SumFac_Pyr::PhysDeriv_SumFac_Pyr
PhysDeriv_SumFac_Pyr(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1703
Nektar::Collections::PhysDeriv_SumFac_Tri::PhysDeriv_SumFac_Tri
PhysDeriv_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:864
Nektar::Collections::PhysDeriv_SumFac_Hex
Phys deriv operator using sum-factorisation (Hex)
Definition: PhysDeriv.cpp:922
Nektar::Collections::PhysDeriv_SumFac_Prism::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1476
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_nquad2
const int m_nquad2
Definition: PhysDeriv.cpp:1694
Nektar::Collections::PhysDeriv_StdMat::PhysDeriv_StdMat
PhysDeriv_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:152
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:199
Nektar::Collections::PhysDeriv_SumFac_Tet::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1263
Nektar::Collections::PhysDeriv_SumFac_Tet::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1265
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1695
Nektar::Collections::PhysDeriv_IterPerExp::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:238
Nektar::Collections::PhysDeriv_SumFac_Tet::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1264
Nektar::Collections::PhysDeriv_SumFac_Pyr::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:1698
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:376
Nektar::Collections::PhysDeriv_SumFac_Tri::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:855
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:250
Nektar::Collections::PhysDeriv_SumFac_Quad
Phys deriv operator using sum-factorisation (Quad)
Definition: PhysDeriv.cpp:611
Nektar::Collections::PhysDeriv_StdMat::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:147
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1064
Nektar::Collections::PhysDeriv_SumFac_Tri::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:856
Nektar::Collections::PhysDeriv_NoCollection::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:389
Nektar::Collections::PhysDeriv_StdMat::operator()
virtual void operator()(const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp)
Perform operation.
Definition: PhysDeriv.cpp:65
Nektar::Collections::PhysDeriv_SumFac_Prism::m_nquad2
const int m_nquad2
Definition: PhysDeriv.cpp:1478
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:186
Nektar::Collections::PhysDeriv_SumFac_Seg::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:577
Nektar::Collections::PhysDeriv_SumFac_Pyr::operator()
virtual void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: PhysDeriv.cpp:1623
Nektar::Collections::PhysDeriv_SumFac_Quad::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:709
Nektar::Collections::PhysDeriv_NoCollection::PhysDeriv_NoCollection
PhysDeriv_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:462
Nektar::Collections::PhysDeriv_IterPerExp::PhysDeriv_IterPerExp
PhysDeriv_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:322
Generated on Wed Jan 20 2021 03:33:39 for Nektar++ by   doxygen 1.8.13