Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: PhysDeriv operator implementations
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <loki/Singleton.h>
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 = m_stdExp->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, false),
207  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Tet"),
208  GetOperatorFactory().RegisterCreatorFunction(
209  OperatorKey(eTetrahedron, ePhysDeriv, eStdMat, true),
210  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_NodalTet"),
211  GetOperatorFactory().RegisterCreatorFunction(
212  OperatorKey(ePyramid, ePhysDeriv, eStdMat, false),
213  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Pyr"),
214  GetOperatorFactory().RegisterCreatorFunction(
215  OperatorKey(ePrism, ePhysDeriv, eStdMat, false),
216  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Prism"),
217  GetOperatorFactory().RegisterCreatorFunction(
218  OperatorKey(ePrism, ePhysDeriv, eStdMat, true),
219  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_NodalPrism"),
220  GetOperatorFactory().RegisterCreatorFunction(
221  OperatorKey(eHexahedron, ePhysDeriv, eStdMat, false),
222  PhysDeriv_StdMat::create, "PhysDeriv_StdMat_Hex")
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 = m_stdExp->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  const int nPhys = m_expList[0]->GetTotPoints();
397  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
398 
399  // calculate local derivatives
400  switch (m_expList[0]->GetShapeDimension())
401  {
402  case 1:
403  {
404  for (int i = 0; i < m_numElmt; ++i)
405  {
406  m_expList[i]->PhysDeriv(input + i*nPhys,
407  tmp0 = output0 + i*nPhys);
408  }
409  break;
410  }
411  case 2:
412  {
413  for (int i = 0; i < m_numElmt; ++i)
414  {
415  m_expList[i]->PhysDeriv(input + i*nPhys,
416  tmp0 = output0 + i*nPhys,
417  tmp1 = output1 + i*nPhys);
418  }
419  break;
420  }
421  case 3:
422  {
423  for (int i = 0; i < m_numElmt; ++i)
424  {
425  m_expList[i]->PhysDeriv(input + i*nPhys,
426  tmp0 = output0 + i*nPhys,
427  tmp1 = output1 + i*nPhys,
428  tmp2 = output2 + i*nPhys);
429  }
430  break;
431  }
432  default:
433  ASSERTL0(false, "Unknown dimension.");
434  }
435  }
436 
437  virtual void operator()(
438  int dir,
439  const Array<OneD, const NekDouble> &input,
440  Array<OneD, NekDouble> &output,
441  Array<OneD, NekDouble> &wsp)
442  {
443  const int nPhys = m_expList[0]->GetTotPoints();
444  Array<OneD, NekDouble> tmp;
445 
446  // calculate local derivatives
447  for (int i = 0; i < m_numElmt; ++i)
448  {
449  m_expList[i]->PhysDeriv(dir, input + i*nPhys,
450  tmp = output + i*nPhys);
451  }
452  }
453 
454  protected:
455  vector<StdRegions::StdExpansionSharedPtr> m_expList;
456 
457  private:
458  PhysDeriv_NoCollection(
459  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
460  CoalescedGeomDataSharedPtr pGeomData)
461  : Operator(pCollExp,pGeomData)
462  {
463  m_expList = pCollExp;
464  }
465 };
466 
467 /// Factory initialisation for the PhysDeriv_NoCollection operators
468 OperatorKey PhysDeriv_NoCollection::m_typeArr[] =
469 {
470  GetOperatorFactory().RegisterCreatorFunction(
471  OperatorKey(eSegment, ePhysDeriv, eNoCollection,false),
472  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Seg"),
473  GetOperatorFactory().RegisterCreatorFunction(
474  OperatorKey(eTriangle, ePhysDeriv, eNoCollection,false),
475  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Tri"),
476  GetOperatorFactory().RegisterCreatorFunction(
477  OperatorKey(eTriangle, ePhysDeriv, eNoCollection,true),
478  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalTri"),
479  GetOperatorFactory().RegisterCreatorFunction(
480  OperatorKey(eQuadrilateral, ePhysDeriv, eNoCollection,false),
481  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Quad"),
482  GetOperatorFactory().RegisterCreatorFunction(
483  OperatorKey(eTetrahedron, ePhysDeriv, eNoCollection,false),
484  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Tet"),
485  GetOperatorFactory().RegisterCreatorFunction(
486  OperatorKey(eTetrahedron, ePhysDeriv, eNoCollection,true),
487  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalTet"),
488  GetOperatorFactory().RegisterCreatorFunction(
489  OperatorKey(ePyramid, ePhysDeriv, eNoCollection,false),
490  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Pyr"),
491  GetOperatorFactory().RegisterCreatorFunction(
492  OperatorKey(ePrism, ePhysDeriv, eNoCollection,false),
493  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Prism"),
494  GetOperatorFactory().RegisterCreatorFunction(
495  OperatorKey(ePrism, ePhysDeriv, eNoCollection,true),
496  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_NodalPrism"),
497  GetOperatorFactory().RegisterCreatorFunction(
498  OperatorKey(eHexahedron, ePhysDeriv, eNoCollection,false),
499  PhysDeriv_NoCollection::create, "PhysDeriv_NoCollection_Hex")
500 };
501 
502 
503 /**
504  * @brief Phys deriv operator using sum-factorisation (Segment)
505  */
506 class PhysDeriv_SumFac_Seg : public Operator
507 {
508  public:
509  OPERATOR_CREATE(PhysDeriv_SumFac_Seg)
510 
511  virtual ~PhysDeriv_SumFac_Seg()
512  {
513  }
514 
515  virtual void operator()(
516  const Array<OneD, const NekDouble> &input,
517  Array<OneD, NekDouble> &output0,
518  Array<OneD, NekDouble> &output1,
519  Array<OneD, NekDouble> &output2,
520  Array<OneD, NekDouble> &wsp)
521  {
522  const int nqcol = m_nquad0*m_numElmt;
523 
524  ASSERTL1(wsp.num_elements() == m_wspSize,
525  "Incorrect workspace size");
526  ASSERTL1(input.num_elements() >= nqcol,
527  "Incorrect input size");
528 
529  Array<OneD, NekDouble> diff0(nqcol, wsp);
530 
531  Blas::Dgemm('N', 'N', m_nquad0, m_numElmt,
532  m_nquad0, 1.0, m_Deriv0, m_nquad0,
533  input.get(), m_nquad0, 0.0,
534  diff0.get(), m_nquad0);
535 
536  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
537 
538  if (m_coordim == 2)
539  {
540  Vmath::Vmul (nqcol, m_derivFac[1], 1, diff0, 1, output1, 1);
541  }
542  else if (m_coordim == 3)
543  {
544  Vmath::Vmul (nqcol, m_derivFac[1], 1, diff0, 1, output1, 1);
545  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output2, 1);
546  }
547  }
548 
549  virtual void operator()(
550  int dir,
551  const Array<OneD, const NekDouble> &input,
552  Array<OneD, NekDouble> &output,
553  Array<OneD, NekDouble> &wsp)
554  {
555  const int nqcol = m_nquad0*m_numElmt;
556 
557  ASSERTL1(wsp.num_elements() == m_wspSize,
558  "Incorrect workspace size");
559  ASSERTL1(input.num_elements() >= nqcol,
560  "Incorrect input size");
561 
562  Array<OneD, NekDouble> diff0(nqcol, wsp);
563 
564  Blas::Dgemm('N', 'N', m_nquad0, m_numElmt,
565  m_nquad0, 1.0, m_Deriv0, m_nquad0,
566  input.get(), m_nquad0, 0.0,
567  diff0.get(), m_nquad0);
568 
569  Vmath::Vmul(nqcol, m_derivFac[dir], 1, diff0, 1, output, 1);
570  }
571 
572  protected:
573  int m_coordim;
574  const int m_nquad0;
575  Array<TwoD, const NekDouble> m_derivFac;
576  NekDouble *m_Deriv0;
577 
578  private:
579  PhysDeriv_SumFac_Seg(
580  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
581  CoalescedGeomDataSharedPtr pGeomData)
582  : Operator (pCollExp, pGeomData),
583  m_nquad0 (m_stdExp->GetNumPoints(0))
584  {
585  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
586  m_coordim = m_stdExp->GetCoordim();
587 
588  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
589 
590  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
591  m_wspSize = m_nquad0*m_numElmt;
592  }
593 
594 };
595 
596 /// Factory initialisation for the PhysDeriv_SumFac_Seg operators
597 OperatorKey PhysDeriv_SumFac_Seg::m_type = GetOperatorFactory().
598  RegisterCreatorFunction(
599  OperatorKey(eSegment, ePhysDeriv, eSumFac,false),
600  PhysDeriv_SumFac_Seg::create, "PhysDeriv_SumFac_Seg");
601 
602 
603 
604 /**
605  * @brief Phys deriv operator using sum-factorisation (Quad)
606  */
607 class PhysDeriv_SumFac_Quad : public Operator
608 {
609  public:
610  OPERATOR_CREATE(PhysDeriv_SumFac_Quad)
611 
612  virtual ~PhysDeriv_SumFac_Quad()
613  {
614  }
615 
616  virtual void operator()(
617  const Array<OneD, const NekDouble> &input,
618  Array<OneD, NekDouble> &output0,
619  Array<OneD, NekDouble> &output1,
620  Array<OneD, NekDouble> &output2,
621  Array<OneD, NekDouble> &wsp)
622  {
623  const int nqtot = m_nquad0 * m_nquad1;
624  const int nqcol = nqtot*m_numElmt;
625 
626  ASSERTL1(wsp.num_elements() == m_wspSize,
627  "Incorrect workspace size");
628  ASSERTL1(input.num_elements() >= nqcol,
629  "Incorrect input size");
630 
631  Array<OneD, NekDouble> diff0(nqcol, wsp );
632  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
633 
634  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
635  m_nquad0, 1.0, m_Deriv0, m_nquad0,
636  input.get(), m_nquad0, 0.0,
637  diff0.get(), m_nquad0);
638 
639  int cnt = 0;
640  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
641  {
642  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
643  input.get() + cnt, m_nquad0,
644  m_Deriv1, m_nquad1, 0.0,
645  diff1.get() + cnt, m_nquad0);
646  }
647 
648  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
649  Vmath::Vvtvp (nqcol, m_derivFac[1], 1, diff1, 1, output0, 1,
650  output0, 1);
651  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output1, 1);
652  Vmath::Vvtvp (nqcol, m_derivFac[3], 1, diff1, 1, output1, 1,
653  output1, 1);
654 
655  if (m_coordim == 3)
656  {
657  Vmath::Vmul (nqcol, m_derivFac[4], 1, diff0, 1, output2, 1);
658  Vmath::Vvtvp (nqcol, m_derivFac[5], 1, diff1, 1, output2, 1,
659  output2, 1);
660  }
661  }
662 
663  virtual void operator()(
664  int dir,
665  const Array<OneD, const NekDouble> &input,
666  Array<OneD, NekDouble> &output,
667  Array<OneD, NekDouble> &wsp)
668  {
669  const int nqtot = m_nquad0 * m_nquad1;
670  const int nqcol = nqtot*m_numElmt;
671 
672  ASSERTL1(wsp.num_elements() == m_wspSize,
673  "Incorrect workspace size");
674  ASSERTL1(input.num_elements() >= nqcol,
675  "Incorrect input size");
676 
677  Array<OneD, NekDouble> diff0(nqcol, wsp );
678  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
679 
680  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
681  m_nquad0, 1.0, m_Deriv0, m_nquad0,
682  input.get(), m_nquad0, 0.0,
683  diff0.get(), m_nquad0);
684 
685  int cnt = 0;
686  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
687  {
688  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
689  input.get() + cnt, m_nquad0,
690  m_Deriv1, m_nquad1, 0.0,
691  diff1.get() + cnt, m_nquad0);
692  }
693 
694  Vmath::Vmul (nqcol, m_derivFac[2*dir] , 1, diff0, 1, output, 1);
695  Vmath::Vvtvp (nqcol, m_derivFac[2*dir+1], 1, diff1, 1, output, 1,
696  output, 1);
697  }
698 
699  protected:
700  int m_coordim;
701  const int m_nquad0;
702  const int m_nquad1;
703  Array<TwoD, const NekDouble> m_derivFac;
704  NekDouble *m_Deriv0;
705  NekDouble *m_Deriv1;
706 
707  private:
708  PhysDeriv_SumFac_Quad(
709  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
710  CoalescedGeomDataSharedPtr pGeomData)
711  : Operator (pCollExp, pGeomData),
712  m_nquad0 (m_stdExp->GetNumPoints(0)),
713  m_nquad1 (m_stdExp->GetNumPoints(1))
714  {
715  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
716  m_coordim = m_stdExp->GetCoordim();
717 
718  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
719 
720  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
721  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
722  m_wspSize = 2 * m_nquad0*m_nquad1*m_numElmt;
723  }
724 };
725 
726 /// Factory initialisation for the PhysDeriv_SumFac_Quad operators
727 OperatorKey PhysDeriv_SumFac_Quad::m_type = GetOperatorFactory().
728  RegisterCreatorFunction(
729  OperatorKey(eQuadrilateral, ePhysDeriv, eSumFac, false),
730  PhysDeriv_SumFac_Quad::create, "PhysDeriv_SumFac_Quad");
731 
732 
733 /**
734  * @brief Phys deriv operator using sum-factorisation (Tri)
735  */
736 class PhysDeriv_SumFac_Tri : public Operator
737 {
738  public:
739  OPERATOR_CREATE(PhysDeriv_SumFac_Tri)
740 
741  virtual ~PhysDeriv_SumFac_Tri()
742  {
743  }
744 
745  virtual void operator()(
746  const Array<OneD, const NekDouble> &input,
747  Array<OneD, NekDouble> &output0,
748  Array<OneD, NekDouble> &output1,
749  Array<OneD, NekDouble> &output2,
750  Array<OneD, NekDouble> &wsp)
751  {
752  const int nqtot = m_nquad0 * m_nquad1;
753  const int nqcol = nqtot*m_numElmt;
754 
755  ASSERTL1(wsp.num_elements() == m_wspSize,
756  "Incorrect workspace size");
757  ASSERTL1(input.num_elements() >= nqcol,
758  "Incorrect input size");
759 
760  Array<OneD, NekDouble> diff0(nqcol, wsp );
761  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
762 
763  // Tensor Product Derivative
764  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
765  m_nquad0, 1.0, m_Deriv0, m_nquad0,
766  input.get(), m_nquad0, 0.0,
767  diff0.get(), m_nquad0);
768 
769  int cnt = 0;
770  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
771  {
772  // scale diff0 by geometric factor: 2/(1-z1)
773  Vmath::Vmul(nqtot,&m_fac1[0],1,diff0.get()+cnt,1,
774  diff0.get()+cnt,1);
775 
776  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
777  input.get() + cnt, m_nquad0,
778  m_Deriv1, m_nquad1, 0.0,
779  diff1.get() + cnt, m_nquad0);
780 
781  // add to diff1 by diff0 scaled by: (1_z0)/(1-z1)
782  Vmath::Vvtvp(nqtot,m_fac0.get(),1,diff0.get()+cnt,1,
783  diff1.get()+cnt,1,diff1.get()+cnt,1);
784  }
785 
786 
787  Vmath::Vmul (nqcol, m_derivFac[0], 1, diff0, 1, output0, 1);
788  Vmath::Vvtvp (nqcol, m_derivFac[1], 1, diff1, 1,
789  output0, 1, output0, 1);
790  Vmath::Vmul (nqcol, m_derivFac[2], 1, diff0, 1, output1, 1);
791  Vmath::Vvtvp (nqcol, m_derivFac[3], 1, diff1, 1,
792  output1, 1, output1, 1);
793 
794  if (m_coordim == 3)
795  {
796  Vmath::Vmul (nqcol, m_derivFac[4], 1, diff0, 1,
797  output2, 1);
798  Vmath::Vvtvp (nqcol, m_derivFac[5], 1, diff1, 1,
799  output2, 1, output2, 1);
800  }
801  }
802 
803  virtual void operator()(
804  int dir,
805  const Array<OneD, const NekDouble> &input,
806  Array<OneD, NekDouble> &output,
807  Array<OneD, NekDouble> &wsp)
808  {
809  const int nqtot = m_nquad0 * m_nquad1;
810  const int nqcol = nqtot*m_numElmt;
811 
812  ASSERTL1(wsp.num_elements() == m_wspSize,
813  "Incorrect workspace size");
814  ASSERTL1(input.num_elements() >= nqcol,
815  "Incorrect input size");
816 
817  Array<OneD, NekDouble> diff0(nqcol, wsp );
818  Array<OneD, NekDouble> diff1(nqcol, wsp + nqcol);
819 
820  // Tensor Product Derivative
821  Blas::Dgemm('N', 'N', m_nquad0, m_nquad1*m_numElmt,
822  m_nquad0, 1.0, m_Deriv0, m_nquad0,
823  input.get(), m_nquad0, 0.0,
824  diff0.get(), m_nquad0);
825 
826  int cnt = 0;
827  for (int i = 0; i < m_numElmt; ++i, cnt += nqtot)
828  {
829  // scale diff0 by geometric factor: 2/(1-z1)
830  Vmath::Vmul(nqtot,&m_fac1[0],1,diff0.get()+cnt,1,
831  diff0.get()+cnt,1);
832 
833  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1, 1.0,
834  input.get() + cnt, m_nquad0,
835  m_Deriv1, m_nquad1, 0.0,
836  diff1.get() + cnt, m_nquad0);
837 
838  // add to diff1 by diff0 scaled by: (1_z0)/(1-z1)
839  Vmath::Vvtvp(nqtot,m_fac0.get(),1,diff0.get()+cnt,1,
840  diff1.get()+cnt,1,diff1.get()+cnt,1);
841  }
842 
843 
844  Vmath::Vmul (nqcol, m_derivFac[2*dir] , 1, diff0, 1, output, 1);
845  Vmath::Vvtvp (nqcol, m_derivFac[2*dir+1], 1, diff1, 1, output, 1,
846  output, 1);
847  }
848 
849  protected:
850  int m_coordim;
851  const int m_nquad0;
852  const int m_nquad1;
853  Array<TwoD, const NekDouble> m_derivFac;
854  NekDouble *m_Deriv0;
855  NekDouble *m_Deriv1;
856  Array<OneD, NekDouble> m_fac0;
857  Array<OneD, NekDouble> m_fac1;
858 
859  private:
860  PhysDeriv_SumFac_Tri(
861  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
862  CoalescedGeomDataSharedPtr pGeomData)
863  : Operator (pCollExp, pGeomData),
864  m_nquad0 (m_stdExp->GetNumPoints(0)),
865  m_nquad1 (m_stdExp->GetNumPoints(1))
866  {
867  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
868  m_coordim = m_stdExp->GetCoordim();
869 
870  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
871 
872  const Array<OneD, const NekDouble>& z0
873  = m_stdExp->GetBasis(0)->GetZ();
874  const Array<OneD, const NekDouble>& z1
875  = m_stdExp->GetBasis(1)->GetZ();
876  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
877  // set up geometric factor: 0.5*(1+z0)
878  for (int i = 0; i < m_nquad0; ++i)
879  {
880  for(int j = 0; j < m_nquad1; ++j)
881  {
882  m_fac0[i+j*m_nquad0] = 0.5*(1+z0[i]);
883  }
884  }
885 
886  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1);
887  // set up geometric factor: 2/(1-z1)
888  for (int i = 0; i < m_nquad0; ++i)
889  {
890  for(int j = 0; j < m_nquad1; ++j)
891  {
892  m_fac1[i+j*m_nquad0] = 2.0/(1-z1[j]);
893  }
894  }
895 
896 
897  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
898  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
899  m_wspSize = 2 * m_nquad0*m_nquad1*m_numElmt;
900  }
901 };
902 
903 /// Factory initialisation for the PhysDeriv_SumFac_Tri operators
904 OperatorKey PhysDeriv_SumFac_Tri::m_typeArr[] =
905 {
906  GetOperatorFactory().RegisterCreatorFunction(
907  OperatorKey(eTriangle, ePhysDeriv, eSumFac,false),
908  PhysDeriv_SumFac_Tri::create, "PhysDeriv_SumFac_Tri"),
909  GetOperatorFactory().RegisterCreatorFunction(
910  OperatorKey(eTriangle, ePhysDeriv, eSumFac,true),
911  PhysDeriv_SumFac_Tri::create, "PhysDeriv_SumFac_NodalTri")
912 };
913 
914 
915 /**
916  * @brief Phys deriv operator using sum-factorisation (Hex)
917  */
918 class PhysDeriv_SumFac_Hex : public Operator
919 {
920  public:
921  OPERATOR_CREATE(PhysDeriv_SumFac_Hex)
922 
923  virtual ~PhysDeriv_SumFac_Hex()
924  {
925  }
926 
927  virtual void operator()(
928  const Array<OneD, const NekDouble> &input,
929  Array<OneD, NekDouble> &output0,
930  Array<OneD, NekDouble> &output1,
931  Array<OneD, NekDouble> &output2,
932  Array<OneD, NekDouble> &wsp)
933  {
934  int nPhys = m_stdExp->GetTotPoints();
935  int ntot = m_numElmt*nPhys;
936  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
937  Array<OneD, Array<OneD, NekDouble> > Diff(3);
938  Array<OneD, Array<OneD, NekDouble> > out(3);
939  out[0] = output0; out[1] = output1; out[2] = output2;
940 
941  for(int i = 0; i < m_dim; ++i)
942  {
943  Diff[i] = wsp + i*ntot;
944  }
945 
946  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
947  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
948  m_nquad0,0.0,&Diff[0][0],m_nquad0);
949 
950  for(int i = 0; i < m_numElmt; ++i)
951  {
952  for (int j = 0; j < m_nquad2; ++j)
953  {
954  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
955  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
956  m_nquad0, m_Deriv1, m_nquad1, 0.0,
957  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
958  m_nquad0);
959  }
960 
961  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
962  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
963  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
964  m_nquad0*m_nquad1);
965  }
966 
967  // calculate full derivative
968  for(int i = 0; i < m_coordim; ++i)
969  {
970  Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
971  for(int j = 1; j < m_dim; ++j)
972  {
973  Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
974  Diff[j], 1,
975  out[i], 1,
976  out[i], 1);
977  }
978  }
979  }
980 
981  virtual void operator()(
982  int dir,
983  const Array<OneD, const NekDouble> &input,
984  Array<OneD, NekDouble> &output,
985  Array<OneD, NekDouble> &wsp)
986  {
987  int nPhys = m_stdExp->GetTotPoints();
988  int ntot = m_numElmt*nPhys;
989  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
990  Array<OneD, Array<OneD, NekDouble> > Diff(3);
991 
992  for(int i = 0; i < m_dim; ++i)
993  {
994  Diff[i] = wsp + i*ntot;
995  }
996 
997  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
998  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
999  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1000 
1001  for(int i = 0; i < m_numElmt; ++i)
1002  {
1003  for (int j = 0; j < m_nquad2; ++j)
1004  {
1005  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1006  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1007  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1008  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1009  m_nquad0);
1010  }
1011 
1012  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1013  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1014  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1015  m_nquad0*m_nquad1);
1016  }
1017 
1018  // calculate full derivative
1019  Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
1020  for(int j = 1; j < m_dim; ++j)
1021  {
1022  Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
1023  Diff[j], 1,
1024  output, 1,
1025  output, 1);
1026  }
1027  }
1028 
1029  protected:
1030  Array<TwoD, const NekDouble> m_derivFac;
1031  int m_dim;
1032  int m_coordim;
1033  const int m_nquad0;
1034  const int m_nquad1;
1035  const int m_nquad2;
1036  NekDouble *m_Deriv0;
1037  NekDouble *m_Deriv1;
1038  NekDouble *m_Deriv2;
1039 
1040  private:
1041  PhysDeriv_SumFac_Hex(
1042  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1043  CoalescedGeomDataSharedPtr pGeomData)
1044  : Operator(pCollExp, pGeomData),
1045  m_nquad0 (m_stdExp->GetNumPoints(0)),
1046  m_nquad1 (m_stdExp->GetNumPoints(1)),
1047  m_nquad2 (m_stdExp->GetNumPoints(2))
1048  {
1049  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1050 
1051  m_dim = PtsKey.size();
1052  m_coordim = m_stdExp->GetCoordim();
1053 
1054  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1055 
1056  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1057  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1058  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1059 
1060  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1061  }
1062 };
1063 
1064 /// Factory initialisation for the PhysDeriv_SumFac_Hex operators
1065 OperatorKey PhysDeriv_SumFac_Hex::m_typeArr[] =
1066 {
1067  GetOperatorFactory().RegisterCreatorFunction(
1068  OperatorKey(eHexahedron, ePhysDeriv, eSumFac, false),
1069  PhysDeriv_SumFac_Hex::create, "PhysDeriv_SumFac_Hex")
1070 };
1071 
1072 
1073 /**
1074  * @brief Phys deriv operator using sum-factorisation (Tet)
1075  */
1076 class PhysDeriv_SumFac_Tet : public Operator
1077 {
1078  public:
1079  OPERATOR_CREATE(PhysDeriv_SumFac_Tet)
1080 
1081  virtual ~PhysDeriv_SumFac_Tet()
1082  {
1083  }
1084 
1085  virtual void operator()(
1086  const Array<OneD, const NekDouble> &input,
1087  Array<OneD, NekDouble> &output0,
1088  Array<OneD, NekDouble> &output1,
1089  Array<OneD, NekDouble> &output2,
1090  Array<OneD, NekDouble> &wsp)
1091  {
1092  int nPhys = m_stdExp->GetTotPoints();
1093  int ntot = m_numElmt*nPhys;
1094  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1095  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1096  Array<OneD, Array<OneD, NekDouble> > out(3);
1097  out[0] = output0; out[1] = output1; out[2] = output2;
1098 
1099  for(int i = 0; i < m_dim; ++i)
1100  {
1101  Diff[i] = wsp + i*ntot;
1102  }
1103 
1104  // dEta0
1105  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1106  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1107  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1108 
1109  // dEta2
1110  for(int i = 0; i < m_numElmt; ++i)
1111  {
1112  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1113  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1114  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1115  m_nquad0*m_nquad1);
1116  }
1117 
1118  for(int i = 0; i < m_numElmt; ++i)
1119  {
1120 
1121  // dEta1
1122  for (int j = 0; j < m_nquad2; ++j)
1123  {
1124  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1125  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1126  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1127  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1128  m_nquad0);
1129  }
1130 
1131  // dxi2 = (1 + eta_1)/(1 -eta_2)*dEta1 + dEta2
1132  Vmath::Vvtvp(nPhys, m_fac3.get(), 1,
1133  Diff[1].get() + i*nPhys, 1,
1134  Diff[2].get() + i*nPhys, 1,
1135  Diff[2].get() + i*nPhys, 1);
1136 
1137  // dxi1 = 2/(1 - eta_2) dEta1
1138  Vmath::Vmul(nPhys, m_fac2.get(), 1,
1139  Diff[1].get() + i*nPhys, 1,
1140  Diff[1].get() + i*nPhys, 1);
1141 
1142  // dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi1
1143  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1144  Diff[0].get() + i*nPhys, 1,
1145  Diff[1].get() + i*nPhys, 1,
1146  Diff[1].get() + i*nPhys, 1);
1147 
1148  // dxi2 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi2
1149  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1150  Diff[0].get() + i*nPhys, 1,
1151  Diff[2].get() + i*nPhys, 1,
1152  Diff[2].get() + i*nPhys, 1);
1153 
1154  // dxi0 = 4.0/((1-eta_1)(1-eta_2)) dEta0
1155  Vmath::Vmul(nPhys, m_fac0.get(), 1,
1156  Diff[0].get() + i*nPhys, 1,
1157  Diff[0].get() + i*nPhys, 1);
1158 
1159  }
1160 
1161  // calculate full derivative
1162  for(int i = 0; i < m_coordim; ++i)
1163  {
1164  Vmath::Vmul(ntot,m_derivFac[i*m_dim],1,Diff[0],1,out[i],1);
1165  for(int j = 1; j < m_dim; ++j)
1166  {
1167  Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+j], 1,
1168  Diff[j], 1, out[i], 1, out[i], 1);
1169  }
1170  }
1171  }
1172 
1173  virtual void operator()(
1174  int dir,
1175  const Array<OneD, const NekDouble> &input,
1176  Array<OneD, NekDouble> &output,
1177  Array<OneD, NekDouble> &wsp)
1178  {
1179  int nPhys = m_stdExp->GetTotPoints();
1180  int ntot = m_numElmt*nPhys;
1181  Array<OneD, NekDouble> tmp0,tmp1,tmp2;
1182  Array<OneD, Array<OneD, NekDouble> > Diff(3);
1183 
1184  for(int i = 0; i < m_dim; ++i)
1185  {
1186  Diff[i] = wsp + i*ntot;
1187  }
1188 
1189  // dEta0
1190  Blas::Dgemm('N','N', m_nquad0,m_nquad1*m_nquad2*m_numElmt,
1191  m_nquad0,1.0, m_Deriv0,m_nquad0,&input[0],
1192  m_nquad0,0.0,&Diff[0][0],m_nquad0);
1193 
1194  // dEta2
1195  for(int i = 0; i < m_numElmt; ++i)
1196  {
1197  Blas::Dgemm('N','T',m_nquad0*m_nquad1,m_nquad2,m_nquad2,
1198  1.0, &input[i*nPhys],m_nquad0*m_nquad1,
1199  m_Deriv2,m_nquad2, 0.0,&Diff[2][i*nPhys],
1200  m_nquad0*m_nquad1);
1201  }
1202 
1203  for(int i = 0; i < m_numElmt; ++i)
1204  {
1205 
1206  // dEta1
1207  for (int j = 0; j < m_nquad2; ++j)
1208  {
1209  Blas::Dgemm('N', 'T', m_nquad0, m_nquad1, m_nquad1,
1210  1.0, &input[i*nPhys+j*m_nquad0*m_nquad1],
1211  m_nquad0, m_Deriv1, m_nquad1, 0.0,
1212  &Diff[1][i*nPhys+j*m_nquad0*m_nquad1],
1213  m_nquad0);
1214  }
1215 
1216  // dxi2 = (1 + eta_1)/(1 -eta_2)*dEta1 + dEta2
1217  Vmath::Vvtvp(nPhys, m_fac3.get(), 1,
1218  Diff[1].get() + i*nPhys, 1,
1219  Diff[2].get() + i*nPhys, 1,
1220  Diff[2].get() + i*nPhys, 1);
1221 
1222  // dxi1 = 2/(1 - eta_2) dEta1
1223  Vmath::Vmul(nPhys, m_fac2.get(), 1,
1224  Diff[1].get() + i*nPhys, 1,
1225  Diff[1].get() + i*nPhys, 1);
1226 
1227  // dxi1 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi1
1228  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1229  Diff[0].get() + i*nPhys, 1,
1230  Diff[1].get() + i*nPhys, 1,
1231  Diff[1].get() + i*nPhys, 1);
1232 
1233  // dxi2 = 2.0(1+eta_0)/((1-eta_1)(1-eta_2)) dEta0 + dxi2
1234  Vmath::Vvtvp(nPhys, m_fac1.get(), 1,
1235  Diff[0].get() + i*nPhys, 1,
1236  Diff[2].get() + i*nPhys, 1,
1237  Diff[2].get() + i*nPhys, 1);
1238 
1239  // dxi0 = 4.0/((1-eta_1)(1-eta_2)) dEta0
1240  Vmath::Vmul(nPhys, m_fac0.get(), 1,
1241  Diff[0].get() + i*nPhys, 1,
1242  Diff[0].get() + i*nPhys, 1);
1243 
1244  }
1245 
1246  // calculate full derivative
1247  Vmath::Vmul(ntot,m_derivFac[dir*m_dim],1,Diff[0],1,output,1);
1248  for(int j = 1; j < m_dim; ++j)
1249  {
1250  Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+j], 1,
1251  Diff[j], 1, output, 1, output, 1);
1252  }
1253  }
1254 
1255  protected:
1256  Array<TwoD, const NekDouble> m_derivFac;
1257  int m_dim;
1258  int m_coordim;
1259  const int m_nquad0;
1260  const int m_nquad1;
1261  const int m_nquad2;
1262  NekDouble *m_Deriv0;
1263  NekDouble *m_Deriv1;
1264  NekDouble *m_Deriv2;
1265  Array<OneD, NekDouble> m_fac0;
1266  Array<OneD, NekDouble> m_fac1;
1267  Array<OneD, NekDouble> m_fac2;
1268  Array<OneD, NekDouble> m_fac3;
1269 
1270  private:
1271  PhysDeriv_SumFac_Tet(
1272  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1273  CoalescedGeomDataSharedPtr pGeomData)
1274  : Operator(pCollExp, pGeomData),
1275  m_nquad0 (m_stdExp->GetNumPoints(0)),
1276  m_nquad1 (m_stdExp->GetNumPoints(1)),
1277  m_nquad2 (m_stdExp->GetNumPoints(2))
1278  {
1279  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1280 
1281  m_dim = PtsKey.size();
1282  m_coordim = m_stdExp->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 < m_dim; ++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*m_dim],1,Diff[0],1,out[i],1);
1405  for(int j = 1; j < m_dim; ++j)
1406  {
1407  Vmath::Vvtvp (ntot, m_derivFac[i*m_dim+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 < m_dim; ++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*m_dim],1,Diff[0],1,output,1);
1466  for(int j = 1; j < m_dim; ++j)
1467  {
1468  Vmath::Vvtvp (ntot, m_derivFac[dir*m_dim+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_dim;
1476  int m_coordim;
1477  const int m_nquad0;
1478  const int m_nquad1;
1479  const int m_nquad2;
1480  NekDouble *m_Deriv0;
1481  NekDouble *m_Deriv1;
1482  NekDouble *m_Deriv2;
1483  Array<OneD, NekDouble> m_fac0;
1484  Array<OneD, NekDouble> m_fac1;
1485 
1486  private:
1487  PhysDeriv_SumFac_Prism(
1488  vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1489  CoalescedGeomDataSharedPtr pGeomData)
1490  : Operator(pCollExp, pGeomData),
1491  m_nquad0 (m_stdExp->GetNumPoints(0)),
1492  m_nquad1 (m_stdExp->GetNumPoints(1)),
1493  m_nquad2 (m_stdExp->GetNumPoints(2))
1494  {
1495  LibUtilities::PointsKeyVector PtsKey = m_stdExp->GetPointsKeys();
1496 
1497  m_dim = PtsKey.size();
1498  m_coordim = m_stdExp->GetCoordim();
1499 
1500  m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1501 
1502  const Array<OneD, const NekDouble>& z0
1503  = m_stdExp->GetBasis(0)->GetZ();
1504  const Array<OneD, const NekDouble>& z2
1505  = m_stdExp->GetBasis(2)->GetZ();
1506  m_fac0 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1507  m_fac1 = Array<OneD, NekDouble>(m_nquad0*m_nquad1*m_nquad2);
1508  for (int i = 0; i < m_nquad0; ++i)
1509  {
1510  for(int j = 0; j < m_nquad1; ++j)
1511  {
1512  for(int k = 0; k < m_nquad2; ++k)
1513  {
1514  m_fac0[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
1515  2.0/(1-z2[k]);
1516  m_fac1[i+j*m_nquad0 + k*m_nquad0*m_nquad1] =
1517  0.5*(1+z0[i]);
1518  }
1519  }
1520  }
1521 
1522 
1523 
1524  m_Deriv0 = &((m_stdExp->GetBasis(0)->GetD())->GetPtr())[0];
1525  m_Deriv1 = &((m_stdExp->GetBasis(1)->GetD())->GetPtr())[0];
1526  m_Deriv2 = &((m_stdExp->GetBasis(2)->GetD())->GetPtr())[0];
1527 
1528  m_wspSize = 3*m_nquad0*m_nquad1*m_nquad2*m_numElmt;
1529  }
1530 };
1531 
1532 /// Factory initialisation for the PhysDeriv_SumFac_Prism operators
1533 OperatorKey PhysDeriv_SumFac_Prism::m_typeArr[] = {
1534  GetOperatorFactory().RegisterCreatorFunction(
1535  OperatorKey(ePrism, ePhysDeriv, eSumFac, false),
1536  PhysDeriv_SumFac_Prism::create, "PhysDeriv_SumFac_Prism")
1537 };
1538 
1539 
1540 }
1541 }
Nektar::Collections::PhysDeriv_SumFac_Prism
Phys deriv operator using sum-factorisation (Prism)
Definition: PhysDeriv.cpp:1338
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
Nektar::Collections::PhysDeriv_StdMat::m_dim
int m_dim
Definition: PhysDeriv.cpp:148
Nektar::LibUtilities::PointsKeyVector
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:242
Nektar::Collections::PhysDeriv_SumFac_Quad::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:702
Nektar::Collections::PhysDeriv_SumFac_Quad::PhysDeriv_SumFac_Quad
PhysDeriv_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:708
Nektar::Collections::PhysDeriv_SumFac_Quad::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:700
Nektar::Collections::PhysDeriv_SumFac_Hex::PhysDeriv_SumFac_Hex
PhysDeriv_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1041
Nektar::Collections::PhysDeriv_SumFac_Tet::PhysDeriv_SumFac_Tet
PhysDeriv_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1271
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_Prism::m_dim
int m_dim
Definition: PhysDeriv.cpp:1475
Nektar::Collections::OperatorKey
boost::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:159
Nektar::Collections::PhysDeriv_IterPerExp::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:317
Nektar::Collections::PhysDeriv_SumFac_Seg::PhysDeriv_SumFac_Seg
PhysDeriv_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:579
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:549
Nektar::MemoryManager
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
Definition: NekMemoryManager.hpp:102
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:745
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:506
Nektar::Collections::PhysDeriv_NoCollection::m_expList
vector< StdRegions::StdExpansionSharedPtr > m_expList
Definition: PhysDeriv.cpp:455
Operator.h
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: PhysDeriv.cpp:1267
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:442
Nektar::Collections::PhysDeriv_SumFac_Hex::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1030
Nektar::Collections::PhysDeriv_SumFac_Prism::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1482
Nektar::Collections::PhysDeriv_SumFac_Quad::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:703
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:1035
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:803
Collection.h
Nektar::Collections::PhysDeriv_SumFac_Tet::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1260
Nektar::Collections::PhysDeriv_SumFac_Prism::PhysDeriv_SumFac_Prism
PhysDeriv_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:1487
Nektar::Collections::PhysDeriv_SumFac_Tri::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:856
Nektar::Collections::PhysDeriv_SumFac_Tri
Phys deriv operator using sum-factorisation (Tri)
Definition: PhysDeriv.cpp:736
Nektar::Collections::Operator
Base class for operators on a collection of elements.
Definition: Operator.h:108
Nektar::Collections::PhysDeriv_SumFac_Hex::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1032
Nektar::Collections::PhysDeriv_SumFac_Prism::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1478
Nektar::Collections::PhysDeriv_SumFac_Seg::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:574
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:616
Nektar::Collections::PhysDeriv_SumFac_Prism::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:1483
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1038
Nektar::Collections::PhysDeriv_SumFac_Tri::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:853
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:1480
Nektar::Collections::PhysDeriv_SumFac_Seg::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:576
Nektar::Collections::PhysDeriv_SumFac_Quad::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:701
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:927
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:1085
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:1261
Nektar::Collections::PhysDeriv_SumFac_Tet::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1259
Nektar::Collections::PhysDeriv_SumFac_Hex::m_dim
int m_dim
Definition: PhysDeriv.cpp:1031
Nektar::Collections::PhysDeriv_SumFac_Seg::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:575
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:981
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:1484
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:1173
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:1476
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:663
Nektar::Collections::PhysDeriv_SumFac_Tri::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:854
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:44
Nektar::Collections::GetOperatorFactory
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:110
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1036
Nektar::Collections::PhysDeriv_SumFac_Hex::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1037
Nektar::Collections::PhysDeriv_SumFac_Tet
Phys deriv operator using sum-factorisation (Tet)
Definition: PhysDeriv.cpp:1076
Nektar::Collections::PhysDeriv_SumFac_Hex::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:1034
Nektar::Collections::PhysDeriv_SumFac_Prism::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:1481
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:515
Nektar::Collections::PhysDeriv_SumFac_Hex::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1033
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: PhysDeriv.cpp:1266
Nektar::Collections::PhysDeriv_SumFac_Quad::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:704
Nektar::Collections::PhysDeriv_SumFac_Tet::m_dim
int m_dim
Definition: PhysDeriv.cpp:1257
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: PhysDeriv.cpp:1265
Nektar::Collections::PhysDeriv_SumFac_Tet::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:1258
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:437
Nektar::Collections::PhysDeriv_SumFac_Tri::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:850
OPERATOR_CREATE
#define OPERATOR_CREATE(cname)
Definition: Operator.h:44
Nektar::Collections::PhysDeriv_SumFac_Tet::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: PhysDeriv.cpp:1256
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:855
Nektar::Collections::PhysDeriv_SumFac_Tet::m_fac3
Array< OneD, NekDouble > m_fac3
Definition: PhysDeriv.cpp:1268
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:857
Nektar::Collections::CoalescedGeomDataSharedPtr
boost::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
Definition: CoalescedGeomData.h:75
Nektar::Collections::PhysDeriv_SumFac_Tri::PhysDeriv_SumFac_Tri
PhysDeriv_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:860
Nektar::Collections::PhysDeriv_SumFac_Hex
Phys deriv operator using sum-factorisation (Hex)
Definition: PhysDeriv.cpp:918
Nektar::Collections::PhysDeriv_SumFac_Prism::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:1477
Nektar::Collections::PhysDeriv_StdMat::PhysDeriv_StdMat
PhysDeriv_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:152
Nektar::Collections::PhysDeriv_SumFac_Tet::m_Deriv0
NekDouble * m_Deriv0
Definition: PhysDeriv.cpp:1262
Nektar::Collections::PhysDeriv_SumFac_Tet::m_Deriv2
NekDouble * m_Deriv2
Definition: PhysDeriv.cpp:1264
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:1263
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:373
Nektar::Collections::PhysDeriv_SumFac_Tri::m_nquad0
const int m_nquad0
Definition: PhysDeriv.cpp:851
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
Nektar::Collections::PhysDeriv_SumFac_Quad
Phys deriv operator using sum-factorisation (Quad)
Definition: PhysDeriv.cpp:607
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:1061
Nektar::Collections::PhysDeriv_SumFac_Tri::m_nquad1
const int m_nquad1
Definition: PhysDeriv.cpp:852
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:1479
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:183
Nektar::Collections::PhysDeriv_SumFac_Seg::m_coordim
int m_coordim
Definition: PhysDeriv.cpp:573
Nektar::Collections::PhysDeriv_SumFac_Quad::m_Deriv1
NekDouble * m_Deriv1
Definition: PhysDeriv.cpp:705
Nektar::Collections::PhysDeriv_NoCollection::PhysDeriv_NoCollection
PhysDeriv_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:458
Nektar::Collections::PhysDeriv_IterPerExp::PhysDeriv_IterPerExp
PhysDeriv_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData)
Definition: PhysDeriv.cpp:322
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215
Generated on Sat Nov 25 2017 16:18:15 for Nektar++ by   doxygen 1.8.8